Last updated: 2021-08-26

Checks: 7 0

Knit directory: IITA_2021GS/

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.

R Markdown file: up-to-date

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.

Environment: empty

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.

Seed: set.seed(20210504)

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.

Session information: recorded

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

Cache: none

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

File paths: relative

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

Repository version: c2c7dae

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 c2c7dae. 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:    code/.DS_Store
    Ignored:    data/.DS_Store

Untracked files:
    Untracked:  data/DatabaseDownload_2021Aug08/
    Untracked:  data/DatabaseDownload_2021May04/
    Untracked:  data/GBSdataMasterList_31818.csv
    Untracked:  data/IITA_GBStoPhenoMaster_33018.csv
    Untracked:  data/NRCRI_GBStoPhenoMaster_40318.csv
    Untracked:  data/PedigreeGeneticGainCycleTime_aafolabi_01122020.xls
    Untracked:  data/Report-DCas21-6038/
    Untracked:  data/blups_forGP.rds
    Untracked:  data/chr1_RefPanelAndGSprogeny_ReadyForGP_72719.fam
    Untracked:  data/dosages_IITA_2021Aug09.rds
    Untracked:  data/haps_IITA_2021Aug09.rds
    Untracked:  data/recombFreqMat_1minus2c_2021Aug02.qs
    Untracked:  output/

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
html 19c3a38 wolfemd 2021-08-19 Build site.
html e029efc wolfemd 2021-08-12 Build site.
Rmd efebeab wolfemd 2021-08-12 Cross-validation and genomic mate predictions complete. All results updated.
html 1c03315 wolfemd 2021-08-11 Build site.
Rmd e4df79f wolfemd 2021-08-11 Completed IITA_2021GS pipeline including imputation and genomic prediction. Last bit of cross-validation and cross-prediction finishes in 24 hrs.
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.

Parent-wise cross-validation

Assess the accuracy of predicted previously unobserved crosses.

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

# 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_forGP.rds")) %>% 
  dplyr::select(-varcomp)

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

# RECOMBINATION FREQUENCY MATRIX (UNFILTERED)
recombFreqMat<-qread(file=here::here("data",
                                     "recombFreqMat_1minus2c_2021Aug02.qs"))

# HAPLOTYPE MATRIX (UNFILTERED)
## 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_2021Aug09.rds"))
parents<-union(ped$sireID,ped$damID) 
parenthaps<-sort(c(paste0(parents,"_HapA"),
                   paste0(parents,"_HapB")))
haploMat<-haploMat[parenthaps,]

# SNP SETS TO ANALYZE
snpsets<-readRDS(file = here::here("data","snpsets.rds"))

# 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) 
# cbsulm14 - 112 cores, 512 GB RAM - 2021 Aug 10 - 1:10pm
snpSet<-"full_set"
grms<-list(A=readRDS(here::here("output","kinship_A_IITA_2021Aug09.rds")), 
           D=readRDS(here::here("output","kinship_Dgeno_IITA_2021Aug09.rds")))
# [1] "Marker-effects Computed. Took  3.08366 hrs"
# [1] "Done predicting fam vars. Took 51.41 mins for 209 crosses"
# [1] "Done predicting fam vars. Took 20.63 mins for 209 crosses"
# [1] "Done predicting fam vars. Took 39.99 mins for 155 crosses"
# [1] "Done predicting fam vars. Took 15.56 mins for 155 crosses"
# [1] "Done predicting fam vars. Took 57.26 mins for 226 crosses"
# [1] "Done predicting fam vars. Took 21.83 mins for 226 crosses"
# [1] "Done predicting fam vars. Took 53.94 mins for 223 crosses"
# [1] "Done predicting fam vars. Took 20.51 mins for 223 crosses"
# [1] "Done predicting fam vars. Took 53.21 mins for 204 crosses"
# [1] "Done predicting fam vars. Took 20.25 mins for 204 crosses"
# [1] "Done predicting fam vars. Took 48.73 mins for 188 crosses"
# [1] "Done predicting fam vars. Took 18.62 mins for 188 crosses"
# [1] "Done predicting fam vars. Took 55.69 mins for 228 crosses"
# [1] "Done predicting fam vars. Took 21.25 mins for 228 crosses"
# [1] "Done predicting fam vars. Took 56.76 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 22.43 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 47.4 mins for 194 crosses"
# [1] "Done predicting fam vars. Took 18.68 mins for 194 crosses"
# [1] "Done predicting fam vars. Took 42.87 mins for 175 crosses"
# [1] "Done predicting fam vars. Took 16.79 mins for 175 crosses"
# [1] "Done predicting fam vars. Took 54.11 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 21.92 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 56.88 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 22.41 mins for 235 crosses"
# ....
# [1] "Accuracies predicted. Took  32.68635 hrs total.Goodbye!"
# [1] "Time elapsed: 1961.205 mins"

# cbsulm30 - 112 cores, 512 GB RAM - 2021 Aug 10 - 2:40pm
snpSet<-"medium_set"
grms<-list(A=readRDS(here::here("output","kinship_A_MediumSNPset_IITA_2021Aug09.rds")), 
           D=readRDS(here::here("output","kinship_Dgeno_MediumSNPset_IITA_2021Aug09.rds")))
# [1] "Marker-effects Computed. Took  1.38038 hrs"
# [1] "Predicting cross variances and covariances"
# Joining, by = c("Repeat", "Fold")
# [1] "Done predicting fam vars. Took 10 mins for 209 crosses"
# [1] "Done predicting fam vars. Took 4.03 mins for 209 crosses"
# [1] "Done predicting fam vars. Took 7.57 mins for 155 crosses"
# [1] "Done predicting fam vars. Took 3.19 mins for 155 crosses"
# [1] "Done predicting fam vars. Took 11.54 mins for 226 crosses"
# [1] "Done predicting fam vars. Took 4.36 mins for 226 crosses"
# [1] "Done predicting fam vars. Took 10.82 mins for 223 crosses"
# [1] "Done predicting fam vars. Took 4.11 mins for 223 crosses"
# [1] "Done predicting fam vars. Took 10.27 mins for 204 crosses"
# [1] "Done predicting fam vars. Took 3.93 mins for 204 crosses"
# [1] "Done predicting fam vars. Took 9.74 mins for 188 crosses"
# [1] "Done predicting fam vars. Took 3.7 mins for 188 crosses"
# [1] "Done predicting fam vars. Took 11 mins for 228 crosses"
# [1] "Done predicting fam vars. Took 4.19 mins for 228 crosses"
# [1] "Done predicting fam vars. Took 11.14 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 4.25 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 9.32 mins for 194 crosses"
# [1] "Done predicting fam vars. Took 3.77 mins for 194 crosses"
# [1] "Done predicting fam vars. Took 8.57 mins for 175 crosses"
# [1] "Done predicting fam vars. Took 3.43 mins for 175 crosses"
# [1] "Accuracies predicted. Took  7.3515 hrs total.Goodbye!"
# > saveRDS(parentWiseCV,
# +         file = here::here("output",
# +                           paste0("parentWiseCV_",snpSet,"_CrossPredAccuracy.rds")))
# > endtime<-proc.time()[3]; print(paste0("Time elapsed: ",
# +                                       round((endtime-starttime)/60,3)," mins"))
# [1] "Time elapsed: 441.107 mins"



# cbsulm17 - 112 cores, 512 GB RAM - 2021 Aug 10 - 1:10pm
snpSet<-"reduced_set"
grms<-list(A=readRDS(here::here("output","kinship_A_ReducedSNPset_IITA_2021Aug09.rds")), 
           D=readRDS(here::here("output","kinship_Dgeno_ReducedSNPset_IITA_2021Aug09.rds")))
# FAILED WITH ERRORS BELOW

“reduced_set” failed with errors in the code chunk below:

snpSet<-"reduced_set"
grms<-list(A=readRDS(here::here("output","kinship_A_ReducedSNPset_IITA_2021Aug09.rds")), 
           D=readRDS(here::here("output","kinship_Dgeno_ReducedSNPset_IITA_2021Aug09.rds")))
# [1] "Time elapsed: 318.057 mins"
# Run `rlang::last_error()` to see where the error occurred.
# > saveRDS(parentWiseCV,
# +         file = here::here("output",
# +                           paste0("parentWiseCV_",snpSet,"_CrossPredAccuracy.rds")))
# Error in saveRDS(parentWiseCV, file = here::here("output", paste0("parentWiseCV_",  :
#   object 'parentWiseCV' not found
# > endtime<-proc.time()[3]; print(paste0("Time elapsed: ",
# +                                       round((endtime-starttime)/60,3)," mins"))
# [1] "Time elapsed: 318.057 mins"
# >
# > rlang::last_error()
# <error/dplyr:::mutate_error>
# Problem with `mutate()` column `predVars`.
# ℹ `predVars = future_map(...)`.
# ✖ error writing to connection
# Backtrace:
#   1. genomicMateSelectR::runParentWiseCrossVal(...)
#  25. base::.handleSimpleError(...)
#  26. dplyr:::h(simpleError(msg, call))
# Run `rlang::last_trace()` to see the full context.
# Backtrace:
#      █
#   1. ├─genomicMateSelectR::runParentWiseCrossVal(...)
#   2. │ └─genomicMateSelectR:::varPredAccuracy(...)
#   3. │   └─`%<>%`(...)
#   4. ├─dplyr::mutate(...)
#   5. ├─dplyr:::mutate.data.frame(...)
#   6. │ └─dplyr:::mutate_cols(.data, ..., caller_env = caller_env())
#   7. │   ├─base::withCallingHandlers(...)
#   8. │   └─mask$eval_all_mutate(quo)
#   9. ├─furrr::future_map(...)
#  10. │ └─furrr:::furrr_map_template(...)
#  11. │   └─furrr:::furrr_template(...)
#  12. │     └─future::future(...)
#  13. │       └─future:::makeFuture(...)
#  14. │         └─future:::strategy(..., envir = envir, workers = workers)
#  15. │           ├─future::run(future)
#  16. │           └─future:::run.ClusterFuture(future)
#  17. │             ├─base::suppressWarnings(...)
#  18. │             │ └─base::withCallingHandlers(...)
#  19. │             └─parallel::clusterCall(cl, fun = gassign, name, value)
#  20. │               └─parallel:::sendCall(cl[[i]], fun, list(...))
#  21. │                 └─parallel:::postNode(...)
#  22. │                   ├─parallel:::sendData(con, list(type = type, data = value, tag = tag))
#  23. │                   └─parallel:::sendData.SOCK0node(...)
#  24. │                     └─base::serialize(data, node$con, xdr = FALSE)
#  25. └─base::.handleSimpleError(...)
#  26.   └─dplyr:::h(simpleError(msg, call))
# <error/simpleError>
# error writing to connection
# >

In the following code chunk, I manually debug / run the component functions of runParentWiseCrossVal() one-by-one only to find no error. Perhaps a server memory issue? Either way, “reduced_set” cross-val. is also completed…

# TRY TO FIX FIND AND FIX THE ERROR THAT KILLED THE REDUCED_SET ANALYSIS
## EVEN THOUGH IT ISN'T A VIP RESULT
## THE ERROR THAT KILLED IT NEEDS FOUND!
snpSet<-"reduced_set"
grms<-list(A=readRDS(here::here("output","kinship_A_ReducedSNPset_IITA_2021Aug09.rds")), 
           D=readRDS(here::here("output","kinship_Dgeno_ReducedSNPset_IITA_2021Aug09.rds")))

snps2keep<-snpsets %>% 
  filter(Set==snpSet) %>% 
  select(snps2keep) %>% 
  unnest(snps2keep)

dosages<-dosages[,snps2keep$FULL_SNP_ID]
haploMat<-haploMat[,snps2keep$FULL_SNP_ID]
recombFreqMat<-recombFreqMat[snps2keep$FULL_SNP_ID,snps2keep$FULL_SNP_ID]
rm(snpsets); gc()

parentfolds<-makeParentFolds(ped=ped,gid="GID",
                             nrepeats=5,
                             nfolds=5,
                             seed=121212)
markEffs<-genomicMateSelectR:::getMarkEffs(parentfolds,blups=blups,gid="GID",modelType="DirDom",
                                           grms=grms,dosages=dosages,
                                           ncores=20,nBLASthreads=5)
saveRDS(parentfolds,file=here::here("output","TMP_parentWiseCV_reduced_set_parentfolds.rds"))
saveRDS(markEffs,file=here::here("output","TMP_parentWiseCV_reduced_set_markerEffects.rds"))

starttime<-proc.time()[3]
cvPredVars<-genomicMateSelectR:::predictCrossVars(modelType="DirDom",snpeffs=markEffs,
                                                  parentfolds=parentfolds,
                                                  haploMat=haploMat,recombFreqMat=recombFreqMat,
                                                  ncores=20,nBLASthreads=5)
# [1] "Done predicting fam vars. Took 4.56 mins for 209 crosses"
# [1] "Done predicting fam vars. Took 2.35 mins for 209 crosses"
# [1] "Done predicting fam vars. Took 3.57 mins for 155 crosses"
# [1] "Done predicting fam vars. Took 1.92 mins for 155 crosses"
# [1] "Done predicting fam vars. Took 4.94 mins for 226 crosses"
# [1] "Done predicting fam vars. Took 2.66 mins for 226 crosses"
# [1] "Done predicting fam vars. Took 4.89 mins for 223 crosses"
# [1] "Done predicting fam vars. Took 2.45 mins for 223 crosses"
# [1] "Done predicting fam vars. Took 4.75 mins for 204 crosses"
# [1] "Done predicting fam vars. Took 2.21 mins for 204 crosses"
# [1] "Done predicting fam vars. Took 3.95 mins for 188 crosses"
# [1] "Done predicting fam vars. Took 2.2 mins for 188 crosses"
# [1] "Done predicting fam vars. Took 4.88 mins for 228 crosses"
# [1] "Done predicting fam vars. Took 2.58 mins for 228 crosses"
# [1] "Done predicting fam vars. Took 5.01 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 2.74 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 4.31 mins for 194 crosses"
# [1] "Done predicting fam vars. Took 2.27 mins for 194 crosses"
# [1] "Done predicting fam vars. Took 3.98 mins for 175 crosses"
# [1] "Done predicting fam vars. Took 2.09 mins for 175 crosses"
# [1] "Done predicting fam vars. Took 4.57 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 2.63 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 4.83 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 2.62 mins for 235 crosses"
# [1] "Done predicting fam vars. Took 4.57 mins for 217 crosses"
# [1] "Done predicting fam vars. Took 2.49 mins for 217 crosses"
# [1] "Done predicting fam vars. Took 3.46 mins for 161 crosses"
# [1] "Done predicting fam vars. Took 1.94 mins for 161 crosses"
# [1] "Done predicting fam vars. Took 4.23 mins for 187 crosses"
# [1] "Done predicting fam vars. Took 2.32 mins for 187 crosses"
# [1] "Done predicting fam vars. Took 4.15 mins for 178 crosses"
# [1] "Done predicting fam vars. Took 2.07 mins for 178 crosses"
# [1] "Done predicting fam vars. Took 4.99 mins for 221 crosses"
# [1] "Done predicting fam vars. Took 2.54 mins for 221 crosses"
# [1] "Done predicting fam vars. Took 4.06 mins for 188 crosses"
# [1] "Done predicting fam vars. Took 2.28 mins for 188 crosses"
# [1] "Done predicting fam vars. Took 5.52 mins for 259 crosses"
# [1] "Done predicting fam vars. Took 2.86 mins for 259 crosses"
# [1] "Done predicting fam vars. Took 3.55 mins for 161 crosses"
# [1] "Done predicting fam vars. Took 1.91 mins for 161 crosses"
# [1] "Done predicting fam vars. Took 4.68 mins for 211 crosses"
# [1] "Done predicting fam vars. Took 2.48 mins for 211 crosses"
# [1] "Done predicting fam vars. Took 4.02 mins for 187 crosses"
# [1] "Done predicting fam vars. Took 2.21 mins for 187 crosses"
# [1] "Done predicting fam vars. Took 4.04 mins for 186 crosses"
# [1] "Done predicting fam vars. Took 2.11 mins for 186 crosses"
# [1] "Done predicting fam vars. Took 4.5 mins for 216 crosses"
# [1] "Done predicting fam vars. Took 2.45 mins for 216 crosses"

cvPredMeans<-genomicMateSelectR:::predictCrossMeans(modelType="DirDom",ncores=20,
                                                    snpeffs=markEffs,
                                                    parentfolds=parentfolds,
                                                    doseMat=dosages)

varPredAcc<-genomicMateSelectR:::varPredAccuracy(modelType = "DirDom",ncores=20,
                                                 crossValOut = cvPredVars,
                                                 snpeffs = markEffs,
                                                 ped = ped,selInd = TRUE,SIwts = SIwts)

## Mean prediction accuracies
meanPredAcc<-genomicMateSelectR:::meanPredAccuracy(modelType = "DirDom",
                                                   crossValOut = cvPredMeans,
                                                   snpeffs = markEffs,
                                                   ped = ped,selInd = TRUE,SIwts = SIwts)

accuracy_out<-list(meanPredAccuracy=meanPredAcc,
                   varPredAccuracy=varPredAcc)
saveRDS(accuracy_out,
        file = here::here("output",
                          paste0("parentWiseCV_",snpSet,"_CrossPredAccuracy.rds")))
# In the end, the code completed successfully.
# There doesn't seem to be a problem with the "reduced_set" or the package...
# Assuming the "full_set" (still running) completes without problems
# It may have been a server error or memory issue.... 
snps2keep<-snpsets %>% 
  filter(Set==snpSet) %>% 
  select(snps2keep) %>% 
  unnest(snps2keep)

dosages<-dosages[,snps2keep$FULL_SNP_ID]
haploMat<-haploMat[,snps2keep$FULL_SNP_ID]
recombFreqMat<-recombFreqMat[snps2keep$FULL_SNP_ID,snps2keep$FULL_SNP_ID]
rm(snpsets); gc()

starttime<-proc.time()[3]
parentWiseCV<-runParentWiseCrossVal(nrepeats=5,nfolds=5,seed=121212,
                                    modelType="DirDom",
                                    ncores=20,nBLASthreads=5,
                                    outName=NULL,
                                    ped=ped,
                                    blups=blups,
                                    dosages=dosages,
                                    haploMat=haploMat,
                                    grms=grms,
                                    recombFreqMat = recombFreqMat,
                                    selInd = TRUE, SIwts = SIwts)
saveRDS(parentWiseCV,
        file = here::here("output",
                          paste0("parentWiseCV_",snpSet,"_CrossPredAccuracy.rds")))
endtime<-proc.time()[3]; print(paste0("Time elapsed: ",
                                      round((endtime-starttime)/60,3)," mins"))

Standard clone-wise cross-validation

Run k-fold cross-validation and assess the accuracy of predicted previously unobserved genotypes (individuals) based on the available training data.

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

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

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

# SNP SETS TO ANALYZE
snpsets<-readRDS(file = here::here("data","snpsets.rds"))

# 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) 
# cbsulm19 - 88 cores, 512 GB RAM - 2021 Aug 10 - 6:45am
snpSet<-"full_set"
grms<-list(A=readRDS(here::here("output","kinship_A_IITA_2021Aug09.rds")), 
           D=readRDS(here::here("output","kinship_Dgeno_IITA_2021Aug09.rds")))
# [1] "Time elapsed: 209.1 mins"

# cbsulm20 - 88 cores, 512 GB RAM - 2021 Aug 10 - 6:45am 
snpSet<-"reduced_set"
grms<-list(A=readRDS(here::here("output","kinship_A_ReducedSNPset_IITA_2021Aug09.rds")), 
           D=readRDS(here::here("output","kinship_Dgeno_ReducedSNPset_IITA_2021Aug09.rds")))
# [1] "Time elapsed: 157.748 mins"

# cbsulm19 - 88 cores, 512 GB RAM - 2021 Aug 10 - 2:45pm
snpSet<-"medium_set"
grms<-list(A=readRDS(here::here("output","kinship_A_MediumSNPset_IITA_2021Aug09.rds")), 
           D=readRDS(here::here("output","kinship_Dgeno_MediumSNPset_IITA_2021Aug09.rds")))
# [1] "Time elapsed: 169.183 mins"
snps2keep<-snpsets %>% 
  filter(Set==snpSet) %>% 
  select(snps2keep) %>% 
  unnest(snps2keep)

dosages<-dosages[,snps2keep$FULL_SNP_ID]
rm(snpsets); gc()

starttime<-proc.time()[3]
standardCV<-runCrossVal(blups=blups,
                        modelType="DirDom",
                        selInd=TRUE,SIwts=SIwts,
                        grms=grms,dosages=dosages,
                        nrepeats=5,nfolds=5,
                        ncores=17,nBLASthreads=5,
                        gid="GID",seed=424242)
saveRDS(standardCV,
        file = here::here("output",
                          paste0("standardCV_",snpSet,"_ClonePredAccuracy.rds")))
endtime<-proc.time()[3]; print(paste0("Time elapsed: ",
                                      round((endtime-starttime)/60,3)," mins"))

Next step / Results

  1. Genomic predictions:
  • First, predict of individual GEBV/GETGV for all selection candidates using all available data and return marker effects for use downstream.
  • Next, Select a top set of candidate parents, for whom we would like to predict cross performances.
  • Finally, predict all pairwise crosses of candidate parents and evaluate them for genomic mate selection.
  • Select the top crosses and plant a crossing nursery with the parents indicated.

Cross-validation results here


sessionInfo()