JiscaH/sequoia

Feature Request: Generation Time

Closed this issue · 5 comments

alexkrohn commented

I have a dataset of tortoises that I'm running through sequoia, then plotting with pedtools. I'm using the true birth years of the tortoises as a prior (birth years range from 1960-2020), but unfortunately the sibship analysis run by sequoia allows parent-offspring relations between tortoises that are only 1-3 years apart. In reality, we know it takes ~10 years for a tortoise to become reproductively active. When I change birth year to "generation" (e.g. 1, 2, 3, 4, etc.), I get many more spurious relationships (five tortoises become their own ancestors).

It seems like using the true birth years as a prior, then just manually changing the resulting pedigree from impossible parent-offspring relationships to sibling relationships resolves this issue. I wonder, however, how complicated it would be to include generation time or minimum age at reproduction, as another optional parameter.

Thanks for the great package. It's proving very useful in my work.

JiscaH commented

Hi Alex,

Actually, I've added this feature just a few weeks back, and it's in the version you can find here on GitHub! Please let me know if/how it works for you; since it's brand new it's not been extensively tested yet.

I'm still tripple checking for bugs, but plan to submit the new version to CRAN this/next week

edit: new parameter is called 'MinAgeParent', see ?MakeAgePrior

alexkrohn commented

What a wonderful coincidence! Thank you. I'll give it a shot.

alexkrohn commented

I encountered an error that I don't quite understand.

sibship <- sequoia(GenoM = geno.low.missing.inds, LifeHistData = life.history.table, Err = 0.05, Module="ped", quiet = FALSE, Plot = TRUE, args.AP = list(MinAgeParent = 10))
Warning:  There are 1058 monomorphic (fixed) SNPs, these will be excluded 

After exclusion, There are  87  individuals and  725  SNPs.

Ageprior: Flat 0/1, overlapping generations, MaxAgeParent = 62,62
Error in if ((Herm != "no" | any(LifeHistData$Sex == 4)) & length(DummyPrefix) ==  : 
  missing value where TRUE/FALSE needed

While the regular call to sequoia used to work, it also gives that same error. Any advice?

sibship <- sequoia(GenoM = geno.low.missing.inds, LifeHistData = life.history.table, Err = 0.05, Module="ped", quiet = FALSE, Plot = TRUE)

sessionInfo()
R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] googlesheets4_1.0.0 forcats_0.5.1       stringr_1.5.0       dplyr_1.1.0         purrr_1.0.0         readr_2.1.1        
 [7] tidyr_1.2.1         tibble_3.1.8        ggplot2_3.4.1       tidyverse_1.3.2     sequoia_2.5.0      

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.10       lubridate_1.8.0   prettyunits_1.1.1 ps_1.6.0          assertthat_0.2.1  rprojroot_2.0.2  
 [7] utf8_1.2.3        R6_2.5.1          cellranger_1.1.0  plyr_1.8.8        backports_1.4.1   reprex_2.0.1     
[13] httr_1.4.2        pillar_1.8.1      rlang_1.0.6       readxl_1.3.1      curl_4.3.2        rstudioapi_0.13  
[19] callr_3.7.0       desc_1.4.0        devtools_2.4.3    googledrive_2.0.0 munsell_0.5.0     broom_0.7.11     
[25] compiler_4.1.2    modelr_0.1.8      askpass_1.1       pkgconfig_2.0.3   pkgbuild_1.3.1    openssl_2.0.0    
[31] tidyselect_1.2.0  fansi_1.0.4       tzdb_0.2.0        crayon_1.5.2      dbplyr_2.1.1      withr_2.5.0      
[37] rappdirs_0.3.3    grid_4.1.2        jsonlite_1.8.4    gtable_0.3.1      lifecycle_1.0.3   DBI_1.1.3        
[43] magrittr_2.0.3    scales_1.2.1      stringi_1.7.12    cli_3.6.0         cachem_1.0.6      fs_1.6.1         
[49] remotes_2.4.2     testthat_3.1.1    xml2_1.3.3        ellipsis_0.3.2    generics_0.1.3    vctrs_0.5.2      
[55] tools_4.1.2       glue_1.6.2        hms_1.1.1         processx_3.5.2    pkgload_1.2.4     fastmap_1.1.0    
[61] colorspace_2.1-0  gargle_1.2.0      sessioninfo_1.2.2 rvest_1.0.2       memoise_2.0.1     haven_2.4.3      
[67] usethis_2.1.5
JiscaH commented

Thanks! That does sound like a bug that should be fairly easy to fix.

The alternative is to edit the age prior matrix 'by hand' in the same way that MinAgeParent does (should do) in the new version, see https://jiscah.github.io/articles/vignette_age/book/sec-customAP.html#other-tweaks

  • 👍1
JiscaH commented

Hi,
Turned out to be just a matter of a forgotten na.rm=TRUE . For now: if you replace all NA's in the Sex column of life.history.table by a 3 , it should work fine.

Cheers, Jisca

  • 👍1