Last updated: 2021-08-02

Checks: 7 0

Knit directory: implementGMSinCassava/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20210504) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 511236f. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/.DS_Store
    Ignored:    analysis/accuracies.png
    Ignored:    analysis/fig2.png
    Ignored:    analysis/fig3.png
    Ignored:    analysis/fig4.png
    Ignored:    analysis/figure/
    Ignored:    code/.DS_Store
    Ignored:    data/.DS_Store

Untracked files:
    Untracked:  accuracies.png
    Untracked:  analysis/docs/
    Untracked:  analysis/speedUpPredCrossVar.Rmd
    Untracked:  data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719.bed
    Untracked:  data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719.bim
    Untracked:  data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719.fam
    Untracked:  data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719.hap.gz
    Untracked:  data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719.log
    Untracked:  data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719.nosex
    Untracked:  data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719.sample
    Untracked:  data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719.vcf.gz
    Untracked:  data/DatabaseDownload_2021May04/
    Untracked:  data/blups_forCrossVal.rds
    Untracked:  data/config.txt
    Untracked:  data/config_mw.txt
    Untracked:  data/dosages_IITA_2021May13.rds
    Untracked:  data/dosages_IITA_filtered_2021May13.rds
    Untracked:  data/genmap_2021May13.rds
    Untracked:  data/haps_IITA_2021May13.rds
    Untracked:  data/haps_IITA_filtered_2021May13.rds
    Untracked:  data/recombFreqMat_1minus2c_2021May13.rds
    Untracked:  fig2.png
    Untracked:  fig3.png
    Untracked:  figure/
    Untracked:  output/Cassava_HMII_V3_Marning_imputation_6-18-21.samples
    Untracked:  output/IITA_CleanedTrialData_2021May10.rds
    Untracked:  output/IITA_ExptDesignsDetected_2021May10.rds
    Untracked:  output/IITA_blupsForModelTraining_twostage_asreml_2021May10.rds
    Untracked:  output/IITA_trials_NOT_identifiable.csv
    Untracked:  output/alphaAssignOutput_iita_pedigree.txt.dams
    Untracked:  output/alphaAssignOutput_iita_pedigree.txt.dams.full
    Untracked:  output/alphaAssignOutput_iita_pedigree.txt.pedigree
    Untracked:  output/alphaAssignOutput_iita_pedigree.txt.pedigree.top
    Untracked:  output/alphaAssignOutput_iita_pedigree.txt.sires
    Untracked:  output/alphaAssignOutput_iita_pedigree.txt.sires.full
    Untracked:  output/crossValPredsA.rds
    Untracked:  output/crossValPredsAD.rds
    Untracked:  output/cvAD_1rep_markerEffects.rds
    Untracked:  output/cvAD_1rep_meanPredAccuracy.rds
    Untracked:  output/cvAD_1rep_parentfolds.rds
    Untracked:  output/cvAD_1rep_predAccuracy.rds
    Untracked:  output/cvAD_1rep_predMeans.rds
    Untracked:  output/cvAD_1rep_predVars.rds
    Untracked:  output/cvAD_1rep_varPredAccuracy.rds
    Untracked:  output/cvAD_5rep5fold_markerEffects.rds
    Untracked:  output/cvAD_5rep5fold_meanPredAccuracy.rds
    Untracked:  output/cvAD_5rep5fold_parentfolds.rds
    Untracked:  output/cvAD_5rep5fold_predMeans.rds
    Untracked:  output/cvAD_5rep5fold_predVars.rds
    Untracked:  output/cvAD_5rep5fold_varPredAccuracy.rds
    Untracked:  output/cvDirDom_5rep5fold_markerEffects.rds
    Untracked:  output/cvDirDom_5rep5fold_meanPredAccuracy.rds
    Untracked:  output/cvDirDom_5rep5fold_parentfolds.rds
    Untracked:  output/cvDirDom_5rep5fold_predMeans.rds
    Untracked:  output/cvDirDom_5rep5fold_predVars.rds
    Untracked:  output/cvDirDom_5rep5fold_varPredAccuracy.rds
    Untracked:  output/cvMeanPredAccuracyA.rds
    Untracked:  output/cvMeanPredAccuracyAD.rds
    Untracked:  output/cvPredMeansA.rds
    Untracked:  output/cvPredMeansAD.rds
    Untracked:  output/cvVarPredAccuracyA.rds
    Untracked:  output/cvVarPredAccuracyAD.rds
    Untracked:  output/genomicMatePredictions_top121parents_ModelAD.csv
    Untracked:  output/genomicMatePredictions_top121parents_ModelAD.rds
    Untracked:  output/genomicMatePredictions_top121parents_ModelDirDom.csv
    Untracked:  output/genomicMatePredictions_top121parents_ModelDirDom.rds
    Untracked:  output/genomicPredictions_ModelAD.csv
    Untracked:  output/genomicPredictions_ModelAD.rds
    Untracked:  output/genomicPredictions_ModelDirDom.csv
    Untracked:  output/genomicPredictions_ModelDirDom.rds
    Untracked:  output/kinship_A_IITA_2021May13.rds
    Untracked:  output/kinship_D_IITA_2021May13.rds
    Untracked:  output/kinship_domGenotypic_IITA_2021July5.rds
    Untracked:  output/markEffsTest.rds
    Untracked:  output/markerEffects.rds
    Untracked:  output/markerEffectsA.rds
    Untracked:  output/markerEffectsAD.rds
    Untracked:  output/maxNOHAV_byStudy.csv
    Untracked:  output/obsCrossMeansAndVars.rds
    Untracked:  output/parentfolds.rds
    Untracked:  output/ped2check_genome.rds
    Untracked:  output/ped2genos.txt
    Untracked:  output/pednames2keep.txt
    Untracked:  output/pednames_Prune100_25_pt25.log
    Untracked:  output/pednames_Prune100_25_pt25.nosex
    Untracked:  output/pednames_Prune100_25_pt25.prune.in
    Untracked:  output/pednames_Prune100_25_pt25.prune.out
    Untracked:  output/potential_dams.txt
    Untracked:  output/potential_sires.txt
    Untracked:  output/predVarTest.rds
    Untracked:  output/samples2keep_IITA_2021May13.txt
    Untracked:  output/samples2keep_IITA_MAFpt01_prune50_25_pt98.log
    Untracked:  output/samples2keep_IITA_MAFpt01_prune50_25_pt98.nosex
    Untracked:  output/samples2keep_IITA_MAFpt01_prune50_25_pt98.prune.in
    Untracked:  output/samples2keep_IITA_MAFpt01_prune50_25_pt98.prune.out
    Untracked:  output/samples2keep_notInPHGdb.txt
    Untracked:  output/stdcv_AD_predAccuracy.rds
    Untracked:  output/stdcv_DirDom_predAccuracy.rds
    Untracked:  output/test_cvAD_markerEffects.rds
    Untracked:  output/test_cvAD_meanPredAccuracy.rds
    Untracked:  output/test_cvAD_parentfolds.rds
    Untracked:  output/test_cvAD_predAccuracy.rds
    Untracked:  output/test_cvAD_predMeans.rds
    Untracked:  output/test_cvAD_predVars.rds
    Untracked:  output/test_cvAD_varPredAccuracy.rds
    Untracked:  output/top50crosses.csv
    Untracked:  output/verified_ped.txt

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/05-CrossValidation.Rmd) and HTML (docs/05-CrossValidation.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd 511236f wolfemd 2021-08-02 Simplify (remove debugging and prelim test sections) and organize for consistency.
Rmd 2d5a542 wolfemd 2021-08-02 Additional debugging of runCrossVal() completed with work shown. Full analysis also completed successfully.
Rmd 230cd14 wolfemd 2021-07-30 Develop upgraded runCrossVal() function to be integrated into code/gmsFunctions.R. Test is and then run full analysis. Debug work shown.
html 769fe81 wolfemd 2021-07-29 Build site.
Rmd f7af63a wolfemd 2021-07-26 Link to and add placeholders for extraction of a VCF from PHG DB and subsequent prepaation of inputs for subsequent analyses.
html 934141c wolfemd 2021-07-14 Build site.
html cc1eb4b wolfemd 2021-07-14 Build site.
Rmd 772750a wolfemd 2021-07-14 DirDom model and selection index calc fully integrated functions.
Rmd 97db806 wolfemd 2021-07-11 Work shown finding and fixing a bug where at least one getMarkEffs model failed. Problem was with use of plan(multicore) + OpenBLAS both using forking. Instead use plan(multisession).
Rmd 2bc9644 wolfemd 2021-07-09 Re-run cross-val with meanPredAccuracy SELIND handling fixed, but debug work not shown anymore.
Rmd 889d98a wolfemd 2021-07-09 test and fix bug in meanPredAccuracy() output when SIwts contain only subset of traits predicted.
Rmd 4308b87 wolfemd 2021-07-08 Full run 5-reps x 5-fold parent-wise cross-val both models DirDom and AD.
Rmd 7888dee wolfemd 2021-07-08 Work fully shown, testing and integrating DirDom model into crossval funcs. Now using R inside a singularity via rocker. Controlling OpenBLAS inside R session with RhpcBLASctl::blas_set_num_threads() and much more.
html 5e45aac wolfemd 2021-06-18 Build site.
Rmd fa20501 wolfemd 2021-06-18 Initial results are ready to publish and share with colleagues for
Rmd 12cc368 wolfemd 2021-06-18 runParentWiseCrossVal for 1 full rep, 5 folds. Found issue with CBSU R compilation but NOT with my code!
html e66bdad wolfemd 2021-06-10 Build site.
Rmd a8452ba wolfemd 2021-06-10 Initial build of the entire page upon completion of all
Rmd 6a5ef32 wolfemd 2021-06-09 meanPredAccuracy() now also included with function moved to “parentWiseCrossVal.R”. NOTE on previous commit: cross-validation functions are NOT in “predCrossVar.R”.
Rmd 63067f7 wolfemd 2021-06-07 Function varPredAccuracy() debugged / tested and moved to predCrossVar.R
Rmd 66c0bde wolfemd 2021-06-07 Remove old and unused code. STILL IN PROGRESS at the computeVarPredAccuracy step.
Rmd 3c085ee wolfemd 2021-06-07 Cross-validation code IN PROGRESS. Currently working on computeVarPredAccuracy.

Previous step

  1. Preprocess data files: Prepare haplotype and dosage matrices, pedigree and BLUPs, genetic map and recombination frequency matrix, for use in predictions.

Also:

  1. Extract and process PHG files: Extract a VCF file from the PHG *.db file produced by Evan Long. Subsequently, prepare haplotype and dosage matrices, genetic map and recombination frequency matrix, for use in predictions.

Automating cross-validation

In the manuscript, the cross-validation is documented many pages and scripts, documented here.

For ongoing GS, I have a function runCrossVal() that manages all inputs and outputs easy to work with pre-computed accuracies.

Goal here is to make a function: runParentWiseCrossVal(), or at least make progress towards developing one.

However, for computational reasons, I imagine it might still be best to separate the task into a few functions.

My goal is to simplify and integrate into the pipeline used for NextGen Cassava. In the paper, used multi-trait Bayesian ridge-regression (MtBRR) to obtain marker effects, and also stored posterior matrices on disk to later compute posterior mean variances. This was computationally expensive and different from my standard univariate REML approach. I think MtBRR and PMV are probably the least biased way to go… but…

For the sake of testing a simple integration into the in-use pipeline, I want to try univariate REML to get the marker effects, which I’ll subsequently use for the cross-validation.

Revised the functions in package:predCrossVar to increase the computational efficiency. Not yet included into the actual R package but instead sourced from code/predCrossVar.R. Additional speed increases were achieved by extra testing to optimize balance of OMP_NUM_THREADS setting (multi-core BLAS) and parallel processing of the crosses-being-predicted. Improvements will benefit users predicting with REML / Bayesian-VPM, but probably worse for Bayesian-PMV.

Set-up server computing env.

Use a a singularity image from the rocker project, as recommended by Qi Sun to get an OpenBLAS-linked R environment that packages can easily be installed on.

This first chunk is one-time only and doesn’t take long. Saves a 650Mb *.sif file to server’s /workdir/

# copy the project data
cd /home/jj332_cas/marnin/;
cp -R implementGMSinCassava /home/$USER;
# the project directory can be in my networked folder for 2 reasons:
# 1) singularity will automatically recognize and be able to access it
# 2) My analyses not read/write intensive; don't break server rules/etiquette 
# set up a working directory on the remote machine
mkdir /workdir/$USER
cd /workdir/$USER/; 

# pull a singularity image and save in the file rocker.sif
# next time you use the rocker.sif file to start the container
singularity pull rocker.sif docker://rocker/tidyverse:latest;

For analysis, operate each R session within a singularity Linux shell within a screen shell.

# 1) start a screen shell 
screen; # or screen -r if re-attaching...
# 2) start the singularity Linux shell inside that
#singularity shell /workdir/$USER/rocker.sif; 
singularity shell ~/rocker2.sif; 
# Project directory, so R will use as working dir.
cd /home/mw489/implementGMSinCassava/;
# 3) Start R
R

Parent-wise cross-validation

Fully-tested runParentWiseCrossVal() and component functions are in the code/parentWiseCrossVal.R script.

Below, source it and use it for a full cross-validation run.

# install.packages(c("RhpcBLASctl","here","rsample","sommer","psych","future.callr","furrr","lme4"))
# install.packages('future.callr')
require(tidyverse); require(magrittr); 
# 5 threads per Rsession for matrix math (openblas)
RhpcBLASctl::blas_set_num_threads(5)

# SOURCE CORE FUNCTIONS
source(here::here("code","parentWiseCrossVal.R"))
source(here::here("code","predCrossVar.R"))

# PEDIGREE
ped<-read.table(here::here("output","verified_ped.txt"),
                header = T, stringsAsFactors = F) %>% 
  rename(GID=FullSampleName,
         damID=DamID,
         sireID=SireID) %>% 
  dplyr::select(GID,sireID,damID)
# Keep only families with _at least_ 2 offspring
ped %<>%
  semi_join(ped %>% count(sireID,damID) %>% filter(n>1) %>% ungroup())

# BLUPs
blups<-readRDS(file=here::here("data","blups_forCrossVal.rds")) %>% 
  dplyr::select(-varcomp)

# GENOMIC RELATIONSHIP MATRICES (GRMS)
grms<-list(A=readRDS(file=here::here("output","kinship_A_IITA_2021May13.rds")),
           D=readRDS(file=here::here("output",
                                     "kinship_domGenotypic_IITA_2021July5.rds")))
## using A+domGenotypic (instead of domClassic used previously)
## will achieve appropriate dom effects for predicting family mean TGV
## but resulting add effects WILL NOT represent allele sub. effects and thus
## predictions won't equal GEBV, allele sub. effects will be post-computed
## as alpha = a + d(q-p)

# DOSAGE MATRIX
dosages<-readRDS(file=here::here("data",
                                 "dosages_IITA_filtered_2021May13.rds"))

# RECOMBINATION FREQUENCY MATRIX
recombFreqMat<-readRDS(file=here::here("data",
                                       "recombFreqMat_1minus2c_2021May13.rds"))
# HAPLOTYPE MATRIX
## keep only haplos for parents-in-the-pedigree
## those which will be used in prediction, saves memory
haploMat<-readRDS(file=here::here("data","haps_IITA_filtered_2021May13.rds"))
parents<-union(ped$sireID,ped$damID) 
parenthaps<-sort(c(paste0(parents,"_HapA"),
                   paste0(parents,"_HapB")))
haploMat<-haploMat[parenthaps,colnames(recombFreqMat)]

# SELECTION INDEX WEIGHTS
## from IYR+IK
## note that not ALL predicted traits are on index
SIwts<-c(logFYLD=20,
         HI=10,
         DM=15,
         MCMDS=-10,
         logRTNO=12,
         logDYLD=20,
         logTOPYLD=15,
         PLTHT=10) 

Full run - both models - parent-wise CV

Server 1: modelType=“AD”

cbsulm29 - 104 cores, 512 GB RAM

grmsAD<-list(A=readRDS(file=here::here("output","kinship_A_IITA_2021May13.rds")),
             D=readRDS(file=here::here("output",
                                       "kinship_D_IITA_2021May13.rds")))
rm(grms)
cvAD_5rep5fold<-runParentWiseCrossVal(nrepeats=5,nfolds=5,seed=84,
                                      modelType="AD",
                                      ncores=20,
                                      outName="output/cvAD_5rep5fold",
                                      ped=ped,
                                      blups=blups,
                                      dosages=dosages,
                                      haploMat=haploMat,
                                      grms=grmsAD,
                                      recombFreqMat = recombFreqMat,
                                      selInd = TRUE, SIwts = SIwts)
saveRDS(cvAD_5rep5fold,here::here("output","cvAD_5rep5fold_predAccuracy.rds"))
# [1] "Marker-effects Computed. Took  1.81086 hrs"
# [1] "Done predicting fam vars. Took 43.11 mins for 198 crosses"
# [1] "Done predicting fam vars. Took 47.04 mins for 216 crosses"
# .....
# [1] "Accuracies predicted. Took  19.68694 hrs total.\n Goodbye!"
# [1] "Accuracies predicted. Took  19.73242 hrs total.Goodbye!"
# > saveRDS(cvAD_5rep5fold,here::here("output","cvAD_5rep5fold_predAccuracy.rds"))

Server 2: modelType=“DirDom”

cbsulm17 - 112 cores, 512 GB RAM

cvDirDom_5rep5fold<-runParentWiseCrossVal(nrepeats=5,nfolds=5,seed=84,
                                          modelType="DirDom",
                                          ncores=20,nBLASthreads=5,
                                          outName="output/cvDirDom_5rep5fold",
                                          ped=ped,
                                          blups=blups,
                                          dosages=dosages,
                                          haploMat=haploMat,
                                          grms=grms,
                                          recombFreqMat = recombFreqMat,
                                          selInd = TRUE, SIwts = SIwts)
saveRDS(cvDirDom_5rep5fold,here::here("output","cvDirDom_5rep5fold_predAccuracy.rds"))
# [1] "Marker-effects Computed. Took  2.3851 hrs"
# [1] "Predicting cross variances and covariances"
# Joining, by = c("Repeat", "Fold")
# [1] "Done predicting fam vars. Took 59.08 mins for 198 crosses"
# [1] "Done predicting fam vars. Took 18.63 mins for 198 crosses"
# [1] "Done predicting fam vars. Took 64.82 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 20.41 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 46.42 mins for 156 crosses"
# [1] "Done predicting fam vars. Took 14.94 mins for 156 crosses"
# [1] "Done predicting fam vars. Took 63.45 mins for 210 crosses"
# [1] "Done predicting fam vars. Took 19.8 mins for 210 crosses"
# [1] "Done predicting fam vars. Took 50.62 mins for 171 crosses"
# [1] "Done predicting fam vars. Took 16.26 mins for 171 crosses"
# [1] "Done predicting fam vars. Took 49.87 mins for 163 crosses"
# [1] "Done predicting fam vars. Took 16.2 mins for 163 crosses"
# [1] "Done predicting fam vars. Took 73.37 mins for 253 crosses"
# [1] "Done predicting fam vars. Took 23.59 mins for 253 crosses"
# [1] "Done predicting fam vars. Took 56.32 mins for 190 crosses"
# [1] "Done predicting fam vars. Took 18.44 mins for 190 crosses"
# [1] "Done predicting fam vars. Took 47.33 mins for 161 crosses"
# [1] "Done predicting fam vars. Took 15.79 mins for 161 crosses"
# [1] "Done predicting fam vars. Took 59.18 mins for 189 crosses"
# [1] "Done predicting fam vars. Took 18.67 mins for 189 crosses"
# [1] "Done predicting fam vars. Took 64.72 mins for 205 crosses"
# [1] "Done predicting fam vars. Took 21.17 mins for 205 crosses"
# [1] "Done predicting fam vars. Took 63.97 mins for 213 crosses"
# [1] "Done predicting fam vars. Took 20.04 mins for 213 crosses"
# [1] "Done predicting fam vars. Took 53.03 mins for 180 crosses"
# [1] "Done predicting fam vars. Took 17.28 mins for 180 crosses"
# [1] "Done predicting fam vars. Took 58.67 mins for 199 crosses"
# [1] "Done predicting fam vars. Took 19.03 mins for 199 crosses"
# ....

# estimate 20 more hours, complete on July 12 very early AM?

# [1] "Accuracies predicted. Took  34.37369 hrs total.Goodbye!"
# Warning message:
# In for (ii in 1L:length(res)) { : closing unused connection 3 (localhost)
# > saveRDS(cvDirDom_5rep5fold,here::here("output","cvDirDom_5rep5fold_predAccuracy.rds"))

Standard clone-wise cross-validation

The new “parent-wise” cross-validation scheme assesses the accuracy of predicting the means and variances of crosses.

The “standard” cross-validation, used in all previous genomic selection analyses (e.g. IITA_2020GS CV, NRCRI_2021GS CV, TARI_2020GS CV), assesses the accuracy predicting the individual performance (breeding value or TGV).

[NEW]: Below, I upgrade the runCrossVal() function used previously for “standard” cross-validation. Include selection index (via selInd= and SIwts= arguments) and a modelType="DirDom" option.

# 1) start a screen shell 
screen; # or screen -r if re-attaching...
# 2) start the singularity Linux shell inside that
#singularity shell /workdir/$USER/rocker.sif; 
singularity shell ~/rocker2.sif; 
# Project directory, so R will use as working dir.
cd /home/mw489/implementGMSinCassava/;
# 3) Start R
R
require(tidyverse); require(magrittr); 

# SOURCE CORE FUNCTIONS
source(here::here("code","gmsFunctions.R"))
source(here::here("code","predCrossVar.R"))

# BLUPs
blups<-readRDS(file=here::here("data","blups_forCrossVal.rds")) %>% 
  dplyr::select(-varcomp) %>% 
  rename(TrainingData=blups) # for compatibility with downstream functions

# SELECTION INDEX WEIGHTS
## from IYR+IK
## note that not ALL predicted traits are on index
SIwts<-c(logFYLD=20,
         HI=10,
         DM=15,
         MCMDS=-10,
         logRTNO=12,
         logDYLD=20,
         logTOPYLD=15,
         PLTHT=10) 

Full run of both models - standard 5-fold CV

cbsulm15 - 112 cores, 512 GB RAM

# GENOMIC RELATIONSHIP MATRICES (GRMS)
grms_ad<-list(A=readRDS(file=here::here("output","kinship_A_IITA_2021May13.rds")),
              D=readRDS(file=here::here("output","kinship_D_IITA_2021May13.rds")))
## achieves the classic TGV=BV+TGV partition
stdcv_ad<-runCrossVal(blups=blups,
                      modelType="AD",
                      selInd=TRUE,SIwts=SIwts,
                      grms=grms_ad,dosages=NULL,
                      nrepeats=5,nfolds=5,
                      ncores=20,nBLASthreads=5,
                      gid="GID",seed=42)
saveRDS(stdcv_ad,here::here("output","stdcv_AD_predAccuracy.rds"))

stdcv_ad %>% 
  select(-splits) %>% 
  unnest(accuracyEstOut) %>% 
  select(-predVSobs)

cbsulm17 - 112 cores, 512 GB RAM

# GENOMIC RELATIONSHIP MATRICES (GRMS)
grms_dirdom<-list(A=readRDS(file=here::here("output","kinship_A_IITA_2021May13.rds")),
                  D=readRDS(file=here::here("output",
                                            "kinship_domGenotypic_IITA_2021July5.rds")))
## using A+domGenotypic (instead of domClassic used previously)
## will achieve appropriate dom effects for predicting family mean TGV
## but resulting add effects WILL NOT represent allele sub. effects and thus
## predictions won't equal GEBV, allele sub. effects will be post-computed
## as alpha = a + d(q-p)

# DOSAGE MATRIX
dosages<-readRDS(file=here::here("data",
                                 "dosages_IITA_filtered_2021May13.rds"))

stdcv_dirdom<-runCrossVal(blups=blups,
                          modelType="DirDom",
                          selInd=TRUE,SIwts=SIwts,
                          grms=grms_dirdom,dosages=dosages,
                          nrepeats=5,nfolds=5,
                          ncores=20,nBLASthreads=5,
                          gid="GID",seed=42)
saveRDS(stdcv_dirdom,here::here("output","stdcv_DirDom_predAccuracy.rds"))

[TO DO] PHG parent-wise cross-validation

Next step / Results

  1. Genomic predictions:
  • A. Standard genomic prediction of individual GEBV and GETGV for all selection candidates using all available data.
  • B. Predict cross means and variances for genomic mate selection

See Results: Home for plots and summary tables.


sessionInfo()