[R-sig-genetics] [poppr] Is there a way to calculate Fst or other metrics alike for snpclone object in poppr?
Emmanuel Paradis
Emmanuel.Paradis at ird.fr
Mon Feb 19 12:11:52 CET 2018
Hi,
The improvement mentioned by Zhian is available in the version of pegas
that is currently on GitHub (which will be the next version on CRAN).
Best,
Emmanuel
Le 15/02/2018 à 16:34, Zhian Kamvar a écrit :
> Hi,
>
> This is more of a question for the R-SIG-GENETICS group as it no longer
> deals with poppr any more, so I'm continuing it here. More knowledgable
> folk, feel free to correct me.
>
>> (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?
> sigma2 represents the variance, so you should use that to estimate
> relative contribution. The next version of pegas will present the phi
> values for you.
>>
>> (2) What is the meaning of "sigma2" and "a"?
> See above.
>>
>> (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 personally put very little trust in p-values for deciding on whether
> or not there is an effect. The choice of 0.05 as the cutoff is
> well-known to be arbitrary, so I think you should use your judgement of
> what you know about the natural history of these populations and your
> sample size to aid your conclusions.
>
> Hope that helps,
> Zhian
>
> -----
> Zhian N. Kamvar, Ph. D.
> Postdoctoral Researcher (Everhart Lab)
> Department of Plant Pathology
> University of Nebraska-Lincoln
> ORCID: 0000-0003-1458-7108
>
>
>
>
>> On Feb 13, 2018, at 23:00 , Yu NING <ningyu_sino at qq.com
>> <mailto: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... at 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:
>>> |
>>> > Bly_gl
>>> /// 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
>>>
>>>
>>>
>>>
>>>> On Jan 23, 2018, at 08:51 , Yu NING <ningy... at qq.com
>>>> <http://qq.com/>> wrote:
>>>>
>>>> Hi~
>>>>
>>>> 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 ? Moreover, If I
>>>> want to export my snpclone object which is MLG filtered,
>>>> how should I code? Thanks very much
>>>>
>>>>
>>>>
>>>> --
>>>> You received this message because you are subscribed to
>>>> the Google Groups "poppr" group.
>>>> To unsubscribe from this group and stop receiving emails
>>>> from it, send an email topoppr+un... at googlegroups.
>>>> <http://googlegroups.com/>com <http://googlegroups.com/>.
>>>> To post to this group, send email
>>>> topo... at googlegroups.com <http://googlegroups.com/>.
>>>> To view this discussion on the web
>>>> visithttps://groups.google.com/d/msgid/poppr/56affc46-45a0-49ab-9511-f38f29b908f3%40googlegroups.com
>>>> <https://groups.google.com/d/msgid/poppr/56affc46-45a0-49ab-9511-f38f29b908f3%40googlegroups.com?utm_medium=email&utm_source=footer>.
>>>> For more options,
>>>> visithttps://groups.google.com/d/optout
>>>> <https://groups.google.com/d/optout>.
>>>
>>>
>>> --
>>> You received this message because you are subscribed to the
>>> Google Groups "poppr" group.
>>> To unsubscribe from this group and stop receiving emails from it,
>>> send an email topoppr+un... at googlegroups.com
>>> <http://googlegroups.com/>.
>>> To post to this group, send email topo... at googlegroups.com
>>> <http://googlegroups.com/>.
>>> To view this discussion on the web
>>> visithttps://groups.google.com/d/msgid/poppr/135db219-5a16-45ed-b32f-487adf73ed41%40googlegroups.com
>>> <https://groups.google.com/d/msgid/poppr/135db219-5a16-45ed-b32f-487adf73ed41%40googlegroups.com?utm_medium=email&utm_source=footer>.
>>> For more options, visithttps://groups.google.com/d/optout
>>> <https://groups.google.com/d/optout>.
>>
>>
>> --
>> You received this message because you are subscribed to the Google
>> Groups "poppr" group.
>> To unsubscribe from this group and stop receiving emails from it, send
>> an email topoppr+unsubscribe at googlegroups.com
>> <mailto:poppr+unsubscribe at googlegroups.com>.
>> To post to this group, send email topoppr at googlegroups.com
>> <mailto:poppr at googlegroups.com>.
>> To view this discussion on the web
>> visithttps://groups.google.com/d/msgid/poppr/0d0f8a6f-807f-47c2-ab72-ef43457da330%40googlegroups.com
>> <https://groups.google.com/d/msgid/poppr/0d0f8a6f-807f-47c2-ab72-ef43457da330%40googlegroups.com?utm_medium=email&utm_source=footer>.
>> For more options, visithttps://groups.google.com/d/optout.
>
>
>
>
>
> _______________________________________________
> R-sig-genetics mailing list
> R-sig-genetics at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-genetics
>
More information about the R-sig-genetics
mailing list