A population is another word for a group of individuals.
For population studies we often want to analyse groups of sequence data.
As we saw in the previous post (Introduction to sequence analysis in R) we can read in sequence data, perform analysis and view the results all in Rstudio.
This post introduces a workflow for running population analysis on sequence data with defined groups.
An Example with real data
We will use the small example dataset from the previous post, available here.
The example fasta file (example.fasta) contains a few COX-1 DNA sequences of fresh water snails collected from different localities.
>HG932219.1
GGCTTATCTTACTAATTCGACTAGAGTTAGGAACTTCCTTAGTACTTATTGATGAACATTTTTATAATGTTATTGTACAGCTCATGCTTTTGTTATGATTTTTTTTATAGGATACCTATAAATTGGTGGGTTTGGGAATTGAATAGTTCCGTTATTAATTGGTGCACCTGATATAAGTTTTCCTCGTATAAATAATATAAGATTTTGATTGCTTCCCCCTTCTTTTATTTTGTTACTTTGTTCTAGAATAGTAGAAGGCGGAGTTGGAACTGGCTGAACAGTCTATCCTCCTCTTAGAGGGCCAATTGCCCACGGAGGTTCTTCAGTAGACTTAGCTATTTTCTCTTTACATTTAGCAGGCCTATCAAGAATTTTAGGAGCTATTAATTTTATTACTACAATTTTAACATACGATCTCCTGGTATTAGATAGAACGAATAAGTTTATTTGTTTGGTCTGTATTGATTACTGCATTTTTATTATTATATCCTTACCTGTTGCAGGGGATTACTATATTATTAACTGATCGAA
>HG932218.1
GGCTTATCTTACTAATTCGACTAGAGTTAGGAACTTCCTTAGTACTTATTGATGAACATTTTTATAATGTTATTGTACAGCTCATGCTTTTGTTATGATTTTTTTTATAGGATACCTATAAATTGGTGGGTTTGGGAATTGAATAGTTCCGTTATTAATTGGTGCACCTGATATAAGTTTTCCTCGTATAAATAATATAAGATTTTGATTGCTTCCCCCTTCTTTTATTTTGTTACTTTGTTCTAGAATAGTAGAAGGCGGAGTTGGAACTGGCTGAACAGTCTATCCTCCTCTTAGAGGGCCAATTGCCCACGGAGGTTCTTCAGTAGACTTAGCTATTTTCTCTTTACATTTAGCAGGCCTATCAAGAATTTTAGGAGCTATTAATTTTATTACTACAATTTTAACATACGATCTCCTGGTATTAGATAGAACGAATAAGTTTATTTGTTTGGTCTGTATTGATTACTGCATTTTTATTATTATATCCTTACCTGTTGGAGGGGATTACTATATTATTAACTGATCGAA
>KP242598.1
GGCTTATCTTACTAATTCGACTAGAGTTAGGAACTTCCTTAGTACTTATTGATGAACATTTTTATAATGTTATTGTACAGCTCATGCTTTTGTTATGATTTTTTTTATAGGATACCTATAAATTGGTGGGTTTGGGAATTGAATAGTTCCGTTACTAATTGGTGCACCTGATATAAGTTTTCCTCGCATAAATAATATAAGATTTTGATTGCTTCCCCCTTCTTTTATTTTGTTACTTTGTTCTAGAATAGTAGAAGGCGGAGTTGGAACTGGTTGAACAGTCTATCCTCCTCTTAGAGGGCCAATTGCCCACGGAGGTTCTTCAGTAGACTTAGCTATTTTTTCTTTACATTTAGCAGGCCTATCAAGAATTTTAGGAGCTATTAATTTTATTACTACAATTTTAACATACGATCTCCTGGTATTAGATAGAACGAATAAGTTTATTTGTTTGGTCTGTATTGATTACTGCATTTTTATTATTATATCCTTACCTGTTGGAGGGGATTACTATATTATTAACTGATCGAA
>KP242599.1
GGCTTATCTTACTAATTCGACTAGAGTTAGGAACTTCCTTAGTACTTATTGATGAACATTTTTATAATGTTATTGTACAGCTCATGCTTTTGTTATGATTTTTTTTATAGGATACCTATAAATTGGTGGGTTTGGGAATTGAATAGTTCCGTTACTAATTGGTGCACCTGATATAAGTTTTCCTCGCATAAATAATATAAGATTTTGATTGCTTCCCCCTTCTTTTATTTTGTTACTTTGTTCTAGAATAGTAGAAGGCGGAGTTGGAACTGGTTGAACAGTCTATCCTCCTCTTAGAGGGCCAATTGCCCACGGAGGTTCTTCAGTAGACTTAGCTATTTTTTCTTTACATTTAGCAGGCCTATCAAGAATTTTAGGAGCTATTAATTTTATTACTACAATTTTAACATACGATCTCCTGGTATTAGATAGAACGAATAAGTTTATTTGTTTGGTCTGTATTGATTACTGCATTTTTATTATTATATCCTTACCTGTTGGAGGGGATTACTATATTATTAACTGATCGAA
>KT337568.1
GGCTTATCTTACTAATTCGACTAGAGTTAGGAACTTCCTTAGTACTTATTGATGAACATTTTTATAATGTTATTGTACAGCTCATGCTTTTGTTATGATTTTTTTTATAGGATACCTATAAATTGGTGGGTTTGGGAATTGAATAGTTCCGTTATTAATTGGTGCACCTGATATAAGTTTTCCTCGCATAAATAATATAAGATTTTGATTGCTTCCCCCTTCTTTTATTTTGTTACTTTGTTCTAGAATAGTAGAAGGCGGAGTTGGAACTGGCTGAACAGTCTATCCTCCTCTTAGAGGGCCAATTGCCCACGGAGGTTCTTCAGTAGACTTAGCTATTTTTTCTTTACATTTAGCAGGCCTATCAAGAATTTTAGGAGCTATTAATTTTATTACTACAATTTTAACATACGATCTCCTGGTATTAGATAGATCGAATAAGTTTATTTGTTTGGTCTGTATTGATTACTGCATTTTTATTATTATATCCTTACCTGTTGGAGGGGATTACTATATTATTAACTGATCGAA
>KT337565.1
GGCTTATCTTACTAATTCGACTAGAGTTAGGAACTTCCTTAGTACTTATTGATGAACATTTTTATAATGTTATTGTACAGCTCATGCTTTTGTTATGATTTTTTTTATAGGATACCTATAAATTGGTGGGTTTGGGAATTGAATAGTTCCGTTATTAATTGGTGCACCTGATATAAGTTTTCCTCGCATAAATAATATAAGATTTTGATTGCTTCCCCCTTCTTTTATTTTGTTACTTTGTTCTAGAATAGTAGAAGGCGGAGTTGGAACTGGCTGAACAGTCTATCCTCCTCTTAGAGGGCCAATTGCCCACGGAGGTTCTTCAGTAGACTTAGCTATTTTTTCTTTACATTTAGCAGGCCTATCAAGAATTTTAGGAGCTATTAATTTTATTACTACAATTTTAACATACGATCTCCTGGTATTAGATAGATCGAATAAGTTTATTTGTTTGGTCTGTATTGATTACTGCATTTTTATTATTATATCCTTACCTGTTGGAGGGGATTACTATATTATTAACTGATCGAAThe associated metadata example_metadata.csv contains information for each sequence, telling us the sequence ID and where each snail was collected from.
For example, according to our metadata, the sequence with accession number KT337565.1 is from a snail collected in the UK which has been assigned to Group 1.
Accession Locality Group
1 KT337565.1 UK 1
2 KT337568.1 UK 1
3 KP242599.1 France 2
4 KP242598.1 France 2
5 HG932219.1 Italy 3
6 HG932218.1 Italy 3Note that the same Accession information is present in both files, this is important and we will use it to link the metadata to our sequence data.
Calculating statistics for geographic populations
To start with we will calculate the nucleotide diversity and haplotype diversity for each group in our example dataset.
There are different ways to define populations.
Groups in this example represent geographic populations.
Below is an example workflow for performing some simple population analyses in R.
1.Read in your data
First we need to read in our sequence alignment and our metadata.
alignment<-read.dna("example.fasta", format = "fasta")
metadata<-read.csv("example_metadata.csv", header = TRUE)2. Sort your data
It’s always a good idea to sort your data, we can use R’s function order() and specify to order the metadata by the ‘Group’ column in our csv file.
metadata_sorted <- metadata[order(metadata$Group),]3. Define your function
Once you have decided what kind of analysis you want to run on each group, first create a function which performs this task on one sequence at a time.
Our pop_stats function below calculates nucleotide diversity and haplotype diversity using the pegas libraries nuc.div() and hap.div().
The function takes a DNA sequence and a group ID and returns the stats in a data frame.
pop_stats<-function(seqs, group){
n<-nuc.div(seqs, variance = FALSE, pairwise.deletion = FALSE)
h<-hap.div(seqs)
stats_out<-data.frame(Group=paste(group),Pi=n, HapDiv=h)
return(stats_out)
}4. Split your data by group
Before we can run our pop_stats function we need to split our metadata by the ‘Group’ column using the split()function.
This creates a list, with an entry for each group, which we can iterate through later:
split_groups<- split(metadata_sorted, metadata_sorted$Group)5. Create a function to call pop_stats
Next we create a group_analysis function that will extract the sequences for each group and call our pop_stats function.
Here we use the Accession ID’s to extract the sequences for each group by their rownames.
These sequences are then passed to the pop_stats function we created earlier, along with their group ID.
group_analysis <- function(each_group) {
sequences <- alignment[rownames(alignment) %in% each_group$Accession, ]
pop_stats(sequences, paste(each_group$Group[1]))
}6. Use lapply and run the group_analysis
Finally we can use lapply() to run the group_analysis function on our list of split_groups:
group_summary_stats<-lapply(split_groups, group_analysis)7.rbind the final result
The final result ( group_summary_stats) is a list, to make the output more reader friendly we can use do.call and rbind to bind the list into a dataframe.
group_results <- do.call(rbind, group_summary_stats)8. Take a look at the output
Now when we call group_results, we can see our population summary stats are in a much more friendly table format.
In this format we can easily save the results into a csv file for later:
write.csv(group_results, file = "group_stats.csv", row.names = FALSE)