---
title: "Parent-wise and standard cross-validation"
author: "Marnin Wolfe"
date: "2021-May-14"
output: 
  workflowr::wflow_html:
    toc: true
editor_options:
  chunk_output_type: inline
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = F, 
                      eval = FALSE, # <- NOTE THAT EVAL SET TO FALSE!
                      tidy='styler', tidy.opts=list(strict=FALSE,width.cutoff=100), highlight=TRUE)
```

# Previous step

4.  [Preprocess data files](04-PreprocessDataFiles.html): 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.**

```{bash parentwiseCV R env, eval=F}
# 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
```
```{r parentwiseCV input}
# 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) 
```
```{r parentwiseCV analysis chunks}
# 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"

# 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"

```
```{r cbsulm17 error on reduced_set}
# cbsulm17 - 112 cores, 512 GB RAM - 2021 Aug 10 - 1:10pm
# FAILED WITH ERRORS 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
# > 

```
```{r cbsulm30 debug the error}
# 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)

cvPredMeans<-genomicMateSelectR:::predictCrossMeans(modelType="DirDom",ncores=20,
                                                    snpeffs=markEffs,
                                                    parentfolds=parentfolds,
                                                    doseMat=dosages)
print(paste0("Took  ",round((proc.time()[3] - starttime)/60/60,5)," hrs"))


```


```{r run parentwiseCV for each snpset}
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.** 

```{bash standardCV R env, eval=F}
# 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
```
```{r standardCV inputs}
# 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) 
```
```{r standardCV analysis chunks}
# 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"

```
```{r run standardCV for each snpset}
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

6.  [Genomic predictions](06-GenomicPredictions.html):   

- 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](07-Results.html#Prediction_accuracy_estimates)



