HWE calculation question
Hi there, you can use seppop to split the data per pop and then use lapply to apply hw.test to each pop. Here's an example using sim2pop: library(adegenet) library(pegas) data(sim2pop)
lapply(seppop(sim2pop), hw.test)
$P01
chi^2 df Pr(chi^2 >) Pr.exact
L01 41.839250 45 0.6066291534 0.261
L02 26.934437 36 0.8629457817 0.523
L03 38.244100 45 0.7517467041 0.686
L04 29.471569 21 0.1031324452 0.075
L05 9.617512 21 0.9834635531 0.942
L06 23.281206 21 0.3291702652 0.241
L07 9.626102 10 0.4738857234 0.626
L08 33.603928 36 0.5830833794 0.578
L09 13.140484 15 0.5914478782 0.448
L10 22.347974 36 0.9634809208 0.971
L11 20.932330 21 0.4630835512 0.249
L12 32.740460 36 0.6244206803 0.446
L13 22.833781 36 0.9567464865 0.909
L14 85.864340 45 0.0002328741 0.002
L15 63.785652 66 0.5543667907 0.831
L16 81.847750 66 0.0902475741 0.220
L17 27.407081 36 0.8476500215 0.530
L18 52.854750 55 0.5570270425 0.402
L19 29.004270 36 0.7895488190 0.245
L20 39.151458 45 0.7171143327 0.652
$P02
chi^2 df Pr(chi^2 >) Pr.exact
L01 37.957211 36 0.38021014 0.450
L02 10.264886 21 0.97526360 0.918
L03 20.045748 15 0.17018475 0.139
L04 25.179973 55 0.99981404 0.954
L05 16.805556 10 0.07877880 0.079
L06 19.361612 15 0.19781251 0.065
L07 1.673554 3 0.64282933 0.839
L08 14.236059 15 0.50770646 0.346
L09 21.052083 28 0.82317084 0.804
L10 29.842469 28 0.37078759 0.120
L11 23.861769 15 0.06747837 0.010
L12 37.108844 36 0.41768024 0.301
L13 12.153169 21 0.93539018 0.959
L14 37.457849 36 0.40208548 0.300
L15 20.244595 36 0.98409332 0.993
L16 25.891684 36 0.89341734 0.663
L17 46.358418 45 0.41607456 0.219
L18 46.151990 45 0.42441641 0.346
L19 26.690888 28 0.53512379 0.433
L20 42.357645 36 0.21570555 0.137
Best
Thibaut
--
Dr Thibaut Jombart
Lecturer, Department of Infectious Disease Epidemiology, Imperial College
London
Head of RECON: repidemicsconsortium.org
sites.google.com/site/thibautjombart/
github.com/thibautjombart
Twitter: @TeebzR <http://twitter.com/TeebzR>
+44(0)20 7594 3658
On 3 March 2017 at 08:57, Emmanuel Paradis <Emmanuel.Paradis at ird.fr> wrote:
Hi Felipe,
Your question is simple but this is not so easy. First you need to convert
your genind object into a loci object with as.loci. Then you can subset
this object for each population and call hw.test() on each subset, but a
problem is that in the "loci" class, each locus is stored as a factor, and
I guess since you have many populations, not all alleles are observed in
all populations. To solve this, you need to drop the extra levels for each
subset. A way to do this is (using the jaguar data in pegas):
library(pegas)
data(jaguar)
for (i in levels(jaguar$population)) {
dat <- jaguar[jaguar$population == i, ]
for (j in attr(dat, "locicol"))
dat[, j] <- factor(dat[, j])
print(hw.test(dat))
}
Note that instead of printing the output of hw.test() you can store it in
a list since this is a data frame.
HTH
Best,
Emmanuel
Le 02/03/2017 ? 17:46, Felipe Hern?ndez a ?crit :
Hi everyone, I have used the pegas package to estimate HWE using a genind object. I got the results from chi-squared and exact test based on Monte Carlo for my whole set of loci (52 loci) across all my 29 populations. I wonder if there is any option to calculate HWE and get exact test Monte Carlo results, but for each population separately. Sorry if the question is so basic, but I would appreciate any helpful advice! Best, Felipe
_______________________________________________ R-sig-genetics mailing list R-sig-genetics at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-genetics