hw.test returning all NA

[dani-davenport]

Dear All,

I am using adegenet on AFLP-like data: Looks like...

Now trying to use hw.test and all values for more then 15,000 loci are NA. Initially I thought maybe some rows/cols of the genind@tab object may only have all 0, and was throwing off the calculations for some reason, but then shouldn't only that loci come out at NA? I tried using seppop(gemindi.obj.obj) %>% lapply(hw.test, B = 0) and was returned the same, all NA.

Am i missing something here?

[thibautjombart]

Yes, the behaviour of tab is no longer what you might have been used to with na.replace. tab does not replace missing values in a genind object, but extracts the data table from it. Check the input of hw.test, it is most likely a genind or a loci object, not a data.frame as returned by tab.

[zkamvar]

I think the problem is a bit more basic:

If you have dominant data, you cannot calculate Hardy-Weinberg Equilibrium. Case in point:

## Created with package reprex version 0.1.1
library("adegenet")
#> Loading required package: ade4
#> 
#>    /// adegenet 2.0.2 is loaded ////////////
#> 
#>    > overview: '?adegenet'
#>    > tutorials/doc/questions: 'adegenetWeb()' 
#>    > bug reports/feature requests: adegenetIssues()
library("pegas")
#> Loading required package: ape
#> 
#> Attaching package: 'pegas'
#> The following object is masked from 'package:ape':
#> 
#>     mst
#> The following object is masked from 'package:ade4':
#> 
#>     amova
dat <- read.table(system.file("files/AFLP.txt", package = "adegenet"), header = TRUE)
obj <- df2genind(dat, ploidy = 2, type = "PA", ncode = 1)
obj
#> /// GENIND OBJECT /////////
#> 
#>  // 7 individuals; 4 loci; 4 alleles; size: 4.2 Kb
#> 
#>  // Basic content
#>    @tab:  7 x 4 matrix of allele counts
#>    @loc.n.all: number of alleles per locus (range: 4-4)
#>    @ploidy: ploidy of each individual  (range: 2-2)
#>    @type:  PA
#>    @call: df2genind(X = dat, ncode = 1, ploidy = 2, type = "PA")
#> 
#>  // Optional content
#>    - empty -
obj <- obj[-4, ]  # removing missing data
pop(obj) <- rep(letters[1:2], 3)  # giving an arbitrary population
hw.test(obj, B = 0)
#>      chi^2 df Pr(chi^2 >)
#> loc1    NA NA          NA
#> loc2    NA NA          NA
#> loc3    NA NA          NA
#> loc4    NA NA          NA

Session info

``` r devtools::session_info() #> Session info -------------------------------------------------------------- #> setting value #> version R version 3.3.3 (2017-03-06) #> system x86_64, darwin13.4.0 #> ui X11 #> language (EN) #> collate en_US.UTF-8 #> tz America/Chicago #> date 2017-03-13 #> Packages ------------------------------------------------------------------ #> package * version date #> ade4 * 1.7-5 2016-12-13 #> adegenet * 2.0.2 2017-02-21 #> ape * 4.1 2017-02-14 #> assertthat 0.1 2013-12-06 #> backports 1.0.5 2017-01-18 #> boot 1.3-18 2016-02-23 #> cluster 2.0.5 2016-10-08 #> coda 0.19-1 2016-12-08 #> colorspace 1.3-2 2016-12-14 #> DBI 0.5-1 2016-09-10 #> deldir 0.1-12 2016-03-06 #> devtools 1.12.0 2016-06-24 #> digest 0.6.12 2017-01-27 #> dplyr 0.5.0 2016-06-24 #> evaluate 0.10 2016-10-11 #> formatR 1.4 2016-05-09 #> gdata 2.17.0 2015-07-04 #> ggplot2 2.2.1 2016-12-30 #> gmodels 2.16.2 2015-07-22 #> gtable 0.2.0 2016-02-26 #> gtools 3.5.0 2015-05-29 #> htmltools 0.3.5 2016-03-21 #> httpuv 1.3.3 2015-08-04 #> igraph 1.0.1 2015-06-26 #> knitr 1.15.13 2017-02-27 #> lattice 0.20-34 2016-09-06 #> lazyeval 0.2.0.9000 2016-07-01 #> LearnBayes 2.15 2014-05-29 #> magrittr 1.5 2014-11-22 #> MASS 7.3-45 2015-11-10 #> Matrix 1.2-8 2017-01-20 #> memoise 1.0.0 2016-01-29 #> mgcv 1.8-17 2017-02-08 #> mime 0.5 2016-07-07 #> munsell 0.4.3 2016-02-13 #> nlme 3.1-131 2017-02-06 #> pegas * 0.9 2016-04-16 #> permute 0.9-4 2016-09-09 #> plyr 1.8.4 2016-06-08 #> R6 2.2.0 2016-10-05 #> Rcpp 0.12.9 2017-01-14 #> reshape2 1.4.2 2016-10-22 #> rmarkdown 1.3.9004 2017-02-27 #> rprojroot 1.2 2017-01-16 #> scales 0.4.1 2016-11-09 #> seqinr 3.3-3 2016-10-13 #> shiny 1.0.0 2017-01-12 #> sp 1.2-4 2016-12-22 #> spdep 0.6-9 2017-01-09 #> stringi 1.1.2 2016-10-01 #> stringr 1.2.0 2017-02-18 #> tibble 1.2 2016-08-26 #> vegan 2.4-2 2017-01-17 #> withr 1.0.2 2016-06-20 #> xtable 1.8-2 2016-02-05 #> yaml 2.1.14 2016-11-12 #> source #> CRAN (R 3.3.2) #> Github (thibautjombart/adegenet@713b873) #> CRAN (R 3.3.2) #> CRAN (R 3.2.0) #> CRAN (R 3.3.2) #> CRAN (R 3.2.3) #> CRAN (R 3.3.0) #> cran (@0.19-1) #> CRAN (R 3.3.2) #> CRAN (R 3.3.0) #> CRAN (R 3.2.4) #> CRAN (R 3.3.0) #> CRAN (R 3.3.2) #> CRAN (R 3.3.0) #> cran (@0.10) #> CRAN (R 3.3.0) #> CRAN (R 3.2.0) #> CRAN (R 3.3.2) #> CRAN (R 3.2.0) #> CRAN (R 3.2.3) #> CRAN (R 3.2.0) #> CRAN (R 3.2.4) #> CRAN (R 3.2.0) #> CRAN (R 3.2.0) #> Github (yihui/knitr@121d9fd) #> CRAN (R 3.3.0) #> Github (hadley/lazyeval@c155c3d) #> CRAN (R 3.2.0) #> CRAN (R 3.2.0) #> CRAN (R 3.2.2) #> CRAN (R 3.3.2) #> CRAN (R 3.2.3) #> CRAN (R 3.3.2) #> cran (@0.5) #> CRAN (R 3.2.3) #> CRAN (R 3.3.2) #> CRAN (R 3.2.5) #> cran (@0.9-4) #> CRAN (R 3.3.0) #> cran (@2.2.0) #> CRAN (R 3.3.2) #> cran (@1.4.2) #> Github (rstudio/rmarkdown@e6cc75e) #> CRAN (R 3.3.2) #> CRAN (R 3.3.2) #> cran (@3.3-3) #> CRAN (R 3.3.2) #> cran (@1.2-4) #> cran (@0.6-9) #> CRAN (R 3.3.0) #> cran (@1.2.0) #> cran (@1.2) #> CRAN (R 3.3.2) #> cran (@1.0.2) #> CRAN (R 3.2.3) #> cran (@2.1.14) ```

[thibautjombart]

"AFLP". It did read "AFLP". My bad. Thanks @zkamvar ;)