Several R packages currently offer AMOVA (analysis of molecular variance) functionality. This post is a work in progress and I will update it as I make progress. My goal here was to obtain locus wise estimates of phi-statistic at two levels of hierarchy:
– Among populations
– Among Regions comprising of populations
Starting data in Structure format:
156 individuals, 14 populations, 2 regions, 5300 SNP sites
Hierarchy levels in text format:
ind pop reg
1 1 1 1
2 2 1 1
3 3 2 1
4 4 2 1
5 5 3 2
6 6 3 2
Let’s load the packages we need:
Import data into genind object:
mygenind <- import2genind('mydata.str', onerowperind=F, n.ind=156, n.loc=5300, col.lab=1, col.pop=2, ask=FALSE)
Convert to genclone object:
mygenclone <- as.genclone(mygenind)
Set hierarchy levels for the genclone object:
myhier <- read.table('myhier.txt', header=T)
sethierarchy(mygenclone) <- myhier
Separate loci using ‘seploc’ function:
mygenclone_seploc <- seploc(mygenclone)
Finally, let’s do AMOVA at two levels: regions and populations
amova_res <- lapply(mygenclone_seploc, poppr.amova, hier=~reg/pop, within=TRUE)
The last command currently does not work. To be updated soon with more info.
If you do not need to do hierarchical amova and simply need phi stats among all populations in your data sets, it is very easy to do. Package
mmod will estimate Meirman’s global phi-stat.
First, load package mmod:
Then do amova directly on the genind object:
amova_phistat_Meirmans <- Phi_st_Meirmans(mygenind)
This will generate global and individual locus phi estimates.
Acknowledgements: Thanks to Zhian Kamvar, an author of poppr package for help with the above code.