On Feb 13, 2018, at 23:00 , Yu NING <ningyu_sino at qq.com> wrote:
Hi, I have some difficulties in interpreting the result of AMOVA. The results show:
Bly_amova
Analysis of Molecular Variance
Call: pegas::amova(formula = Blygl_dis ~ popslot, data = strata(Bly_gl),
is.squared = TRUE)
SSD MSD df
popslot 40434.10 3110.316 13
Error 36754.33 1470.173 25
Total 77188.44 2031.275 38
Variance components:
sigma2 P.value
popslot 589.75 0.048
Error 1470.17
Variance coefficients:
a
2.781065
I wonder:
(1) Could I estimate the relative contribution of variances using SSD? e.g. 40434.10/77188.44=52.38%, It's the contribution of between groups variances?
(2) What is the meaning of "sigma2" and "a"?
(3) Is P<0.05 conventionally enough to conclude significant difference for molecular data ? My result is 0.048 but I doubt it's power according to field observation. Should I set more strict critiria? 0.01 for example
I've tried to searched for pegas tutorial or forum but failed. Any advice here will be highly appreciated. : )
? 2018?2?12???? UTC+8??11:31:06?Zhian Kamvar???
Hello,
I think I know why you were having this problem:
1. You didn't have anything defined in your strata [1].
2. You didn't have a variable called "pop" defined in your workspace
3. Pegas saw a variable called "pop", which was a function, and gave a warning ("elements in the rhs of the formula are not all factors")
However, given your code below, I'm not exactly sure why it works because "popslot" is not the same as "popslot_Kob". This tells me that you may run into problems in the future.
You can define your strata with a data frame where each row represents a sample and each column is a different population grouping. A simple way to define it for the population:
strata(Kob_snpclo) <- data.frame(population = pop(Kob_snpclo))
From there, you would use "population" on the right hand side of the formula in AMOVA.
Cheers,
Zhian
[1]: If you are unfamiliar with strata, check out the strata tutorial: https://github.com/thibautjombart/adegenet/wiki/Tutorials <https://github.com/thibautjombart/adegenet/wiki/Tutorials>
-----
Zhian N. Kamvar, Ph. D.
Postdoctoral Researcher (Everhart Lab)
Department of Plant Pathology
University of Nebraska-Lincoln
ORCID: 0000-0003-1458-7108
On Feb 12, 2018, at 04:47 , Yu NING <ningy...@ <>qq.com <http://qq.com/>> wrote:
Now I have found one way to work it out:
popslot_Kob<-pop(Kob_snpclo)
Blygl_dis<- bitwise.dist(Bly_gl, percent = FALSE)
Bly_amova<- pegas::amova(Blygl_dis ~ popslot, data = strata(Bly_gl), is.squared = TRUE)
Just add one line to extract the pop information. But I still wonder why error showed when I use:
Bly_amova<- pegas::amova(Blygl_dis ~ pop, data = strata(Bly_gl), is.squared = TRUE)
It's a little wired....
? 2018?2?11???? UTC+8??11:26:11?Yu NING???
Hello,
Thanks for your help. But other problems showed up:
Blygl_dis<- bitwise.dist(Bly_gl, percent = FALSE)
Bly_amova<- pegas::amova(Blygl_dis ~ pop, data = strata(Bly_gl), is.squared = TRUE)
Error in FUN(X[[i]], ...) : 'bin' must be numeric or a factor
In addition: Warning message:
In pegas::amova(Blygl_dis ~ pop, data = strata(Bly_gl), is.squared = TRUE) :
elements in the rhs of the formula are not all factors
Is there something wrong in my data organization? This is the data info:
/// GENLIGHT OBJECT /////////
// 39 genotypes, 7,227 binary SNPs, size: 969.5 Kb
11255 (3.99 %) missing data
// Basic content
@gen: list of 39 SNPbin
@ploidy: ploidy of each individual (range: 2-2)
// Optional content
@ind.names: 39 individual labels
@loc.names: 7227 locus labels
@chromosome: factor storing chromosomes of the SNPs
@position: integer storing positions of the SNPs
@pop: population of each individual (group size range: 2-3)
@other: a list containing: elements without names
? 2018?1?23???? UTC+8??11:38:39?Zhian Kamvar???
Hello,
I have browsed the manual but I couldn't find how to calculate metrics of differentiation in poppr , and poppr.amova() is also not applicable for snpclone object. Is there a way to solve it ?
Unfortunately, calculating AMOVA on genlight objects is not currently possible in poppr [1]. I would recommend to use the pegas implementation like so:
library("pegas")
myDist <- bitwise.dist(myGenlight, percent = FALSE)
res <- pegas::amova(myDist ~ pop/subpop, data = strata(myGenlight), is.squared = TRUE)
Replace "myGenlight" with the name of your data and "pop/subpop" with the hierarchy present in your strata slot.
Moreover, If I want to export my snpclone object which is MLG filtered, how should I code?
If you want to export the object for use in another Rscript, I would recommend using saveRDS() to save the object to a file and then readRDS() to read that file in a new session*. If you only need the MLG assignments, then he underlying genetic data does not change when you filter MLGs, so you can export the assignments as a CSV file along with your individual names and population data.
I hope that helps.
Zhian
[1]: https://github.com/grunwaldlab/poppr/issues/72 <https://github.com/grunwaldlab/poppr/issues/72>. This is not a simple problem and may still take a bit longer to implement as I am now only working on poppr in my spare time.
* If you defined your MLGs using a matrix, then you must also save that matrix and reload it with the same name.
-----
Zhian N. Kamvar, Ph. D.
Postdoctoral Researcher (Everhart Lab)
Department of Plant Pathology
University of Nebraska-Lincoln
ORCID: 0000-0003-1458-7108