Is there a way to calculate Fst or other metrics alike for snpclone object in poppr?

[Yu NING]

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

[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. 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...@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 to poppr+un...@googlegroups.com.
To post to this group, send email to po...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/poppr/56affc46-45a0-49ab-9511-f38f29b908f3%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[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写道：

[Yu NING]

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写道：

[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

-----

Zhian N. Kamvar, Ph. D.

Postdoctoral Researcher (Everhart Lab)

Department of Plant Pathology

University of Nebraska-Lincoln

ORCID: 0000-0003-1458-7108

To view this discussion on the web visit https://groups.google.com/d/msgid/poppr/135db219-5a16-45ed-b32f-487adf73ed41%40googlegroups.com.

[Yu NING]

Oops, I pasted the wrong line. Actually these codes work:

popslot<-pop(Bly_snpclo)

Blygl_dis<- bitwise.dist(Bly_gl, percent = FALSE)

Bly_amova<- pegas::amova(Blygl_dis ~ popslot, data = strata(Bly_gl), is.squared = TRUE)

Thanks for your advice. I'll go through the strata() documentation. 

在 2018年2月12日星期一 UTC+8下午11:31:06，Zhian Kamvar写道：

[Yu NING]

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,

[Zhian Kamvar]

Hello all,

I punted this thread over to R-SIG-GENTICS: https://stat.ethz.ch/pipermail/r-sig-genetics/2018-February/000181.html

I would recommend signing up for this list. It's quite low volume and allows you to ask general population genetics questions to a wider community. 

Best,

Zhian

-----

Zhian N. Kamvar, Ph. D.

Postdoctoral Researcher (Everhart Lab)

Department of Plant Pathology

University of Nebraska-Lincoln

ORCID: 0000-0003-1458-7108

To view this discussion on the web visit https://groups.google.com/d/msgid/poppr/0d0f8a6f-807f-47c2-ab72-ef43457da330%40googlegroups.com.

[Yu NING]

Big thanks for all the help, Zhian. I am joining the list ^_^

在 2018年2月15日星期四 UTC+8下午11:41:44，Zhian Kamvar写道：