---
title: "Parent-wise cross-validation to check the accuracy of predicting cross (co)-variances"
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.

Also:

5.  [Extract and process PHG files](08-PHGfiles.html): 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](https://wolfemd.github.io/PredictOutbredCrossVar/).

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/**
```{bash, eval=F}
# 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.
```{bash set-up R environment, eval=F}
# 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. 

```{r install packages if needed}
# install.packages(c("RhpcBLASctl","here","rsample","sommer","psych","future.callr","furrr","lme4","qs"))
# install.packages('future.callr')
```

```{r parent-wise CV inputs}
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

```{r runParentWiseCrossVal - model AD}
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
```{r runParentWiseCrossVal - model DirDom}
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](https://wolfemd.github.io/IITA_2020GS/06-Results.html#cross-validation-accuracy), [NRCRI_2021GS CV](https://wolfemd.github.io/NRCRI_2021GS/05-Results.html#Prediction_accuracy), [TARI_2020GS CV](https://wolfemd.github.io/TARI_2020GS/05-Results.html#Prediction_accuracy)), 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.

```{bash R environment, eval=F}
# 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
```
```{r clone-wise CV inputs}
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
```{r runCrossVal - model AD}
# 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
```{r runCrossVal - model DirDom}
# 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"))

```

# Optimize marker density for speed

```{bash usual R environment, eval=F}
# 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
```

Cross-variance prediction is slow, but significant speed gains can be made by using fewer markers for the predictions. Examine the speed benefit vs. accuracy cost trade-off.

## LD-pruning to reduce SNP numbers

Two ways to accomplish reduced marker number: (1) random sample, (2) LD-pruning ([as in the preprocessing step _here_](https://wolfemd.github.io/implementGMSinCassava/04-PreprocessDataFiles.html#Variant_filters)). 

Try both below.

Use cross-validation to compare accuracy achieved.

**DirDom** model-only at this point.

Use more stringent versions of the `plink1.9 --indep-pairwise` used previously ([here](https://wolfemd.github.io/implementGMSinCassava/04-PreprocessDataFiles.html#Variant_filters))

```{plink, eval=F}
cd ~/implementGMSinCassava/
export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;
# plink --bfile data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719 \
#   --keep output/samples2keep_IITA_2021May13.txt \
#   --indep-pairwise 250 'kb' 125 0.9 \
#   --out output/samples2keep_IITA_MAFpt01_prune250kb_125_pt9; 
#  24986 of 68068 variants removed. (43082 left, more than original filter)

plink --bfile data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719 \
  --keep output/samples2keep_IITA_2021May13.txt \
  --indep-pairwise 1000 'kb' 50 0.9 \
  --out output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt9; 
# Pruned 6587 variants from chromosome 1, leaving 1720.
# Pruned 2524 variants from chromosome 2, leaving 1457.
# Pruned 2904 variants from chromosome 3, leaving 1350.
# Pruned 3110 variants from chromosome 4, leaving 1006.
# Pruned 2918 variants from chromosome 5, leaving 1222.
# Pruned 2510 variants from chromosome 6, leaving 1142.
# Pruned 1288 variants from chromosome 7, leaving 800.
# Pruned 2454 variants from chromosome 8, leaving 1311.
# Pruned 2557 variants from chromosome 9, leaving 971.
# Pruned 1889 variants from chromosome 10, leaving 1123.
# Pruned 2225 variants from chromosome 11, leaving 1175.
# Pruned 1982 variants from chromosome 12, leaving 990.
# Pruned 2059 variants from chromosome 13, leaving 879.
# Pruned 4004 variants from chromosome 14, leaving 1540.
# Pruned 2877 variants from chromosome 15, leaving 1097.
# Pruned 2059 variants from chromosome 16, leaving 835.
# Pruned 1996 variants from chromosome 17, leaving 976.
# Pruned 2226 variants from chromosome 18, leaving 1051.
# Pruning complete.  48169 of 68814 variants removed. (20645 left) 

plink --bfile data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719 \
  --keep output/samples2keep_IITA_2021May13.txt \
  --indep-pairwise 1000 'kb' 50 0.8 \
  --out output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt8; 
# Pruned 7131 variants from chromosome 1, leaving 1176.
# Pruned 2888 variants from chromosome 2, leaving 1093.
# Pruned 3287 variants from chromosome 3, leaving 967.
# Pruned 3358 variants from chromosome 4, leaving 758.
# Pruned 3229 variants from chromosome 5, leaving 911.
# Pruned 2732 variants from chromosome 6, leaving 920.
# Pruned 1438 variants from chromosome 7, leaving 650.
# Pruned 2738 variants from chromosome 8, leaving 1027.
# Pruned 2826 variants from chromosome 9, leaving 702.
# Pruned 2190 variants from chromosome 10, leaving 822.
# Pruned 2506 variants from chromosome 11, leaving 894.
# Pruned 2236 variants from chromosome 12, leaving 736.
# Pruned 2268 variants from chromosome 13, leaving 670.
# Pruned 4448 variants from chromosome 14, leaving 1096.
# Pruned 3188 variants from chromosome 15, leaving 786.
# Pruned 2256 variants from chromosome 16, leaving 638.
# Pruned 2247 variants from chromosome 17, leaving 725.
# Pruned 2506 variants from chromosome 18, leaving 771.
# Pruning complete.  53472 of 68814 variants removed. (15342 left)

plink --bfile data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719 \
  --keep output/samples2keep_IITA_2021May13.txt \
  --indep-pairwise 1000 'kb' 50 0.7 \
  --out output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt7; 
# Pruned 7413 variants from chromosome 1, leaving 894.
# Pruned 3118 variants from chromosome 2, leaving 863.
# Pruned 3504 variants from chromosome 3, leaving 750.
# Pruned 3484 variants from chromosome 4, leaving 632.
# Pruned 3426 variants from chromosome 5, leaving 714.
# Pruned 2913 variants from chromosome 6, leaving 739.
# Pruned 1555 variants from chromosome 7, leaving 533.
# Pruned 2932 variants from chromosome 8, leaving 833.
# Pruned 2972 variants from chromosome 9, leaving 556.
# Pruned 2372 variants from chromosome 10, leaving 640.
# Pruned 2709 variants from chromosome 11, leaving 691.
# Pruned 2381 variants from chromosome 12, leaving 591.
# Pruned 2411 variants from chromosome 13, leaving 527.
# Pruned 4736 variants from chromosome 14, leaving 808.
# Pruned 3370 variants from chromosome 15, leaving 604.
# Pruned 2372 variants from chromosome 16, leaving 522.
# Pruned 2396 variants from chromosome 17, leaving 576.
# Pruned 2663 variants from chromosome 18, leaving 614.
# Pruning complete.  56727 of 68814 variants removed. (12087 left)

plink --bfile data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719 \
  --keep output/samples2keep_IITA_2021May13.txt \
  --indep-pairwise 1000 'kb' 50 0.6 \
  --out output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt6; 
# Pruned 7619 variants from chromosome 1, leaving 688.
# Pruned 3354 variants from chromosome 2, leaving 627.
# Pruned 3666 variants from chromosome 3, leaving 588.
# Pruned 3596 variants from chromosome 4, leaving 520.
# Pruned 3597 variants from chromosome 5, leaving 543.
# Pruned 3060 variants from chromosome 6, leaving 592.
# Pruned 1670 variants from chromosome 7, leaving 418.
# Pruned 3094 variants from chromosome 8, leaving 671.
# Pruned 3073 variants from chromosome 9, leaving 455.
# Pruned 2511 variants from chromosome 10, leaving 501.
# Pruned 2882 variants from chromosome 11, leaving 518.
# Pruned 2493 variants from chromosome 12, leaving 479.
# Pruned 2505 variants from chromosome 13, leaving 433.
# Pruned 4929 variants from chromosome 14, leaving 615.
# Pruned 3528 variants from chromosome 15, leaving 446.
# Pruned 2455 variants from chromosome 16, leaving 439.
# Pruned 2508 variants from chromosome 17, leaving 464.
# Pruned 2773 variants from chromosome 18, leaving 504.
# Pruning complete.  59313 of 68814 variants removed. (9501 left)

plink --bfile data/AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719 \
  --keep output/samples2keep_IITA_2021May13.txt \
  --indep-pairwise 1000 'kb' 50 0.5 \
  --out output/samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt5; 
# Pruned 7791 variants from chromosome 1, leaving 516.
# Pruned 3497 variants from chromosome 2, leaving 484.
# Pruned 3800 variants from chromosome 3, leaving 454.
# Pruned 3688 variants from chromosome 4, leaving 428.
# Pruned 3722 variants from chromosome 5, leaving 418.
# Pruned 3205 variants from chromosome 6, leaving 447.
# Pruned 1747 variants from chromosome 7, leaving 341.
# Pruned 3228 variants from chromosome 8, leaving 537.
# Pruned 3161 variants from chromosome 9, leaving 367.
# Pruned 2622 variants from chromosome 10, leaving 390.
# Pruned 2992 variants from chromosome 11, leaving 408.
# Pruned 2589 variants from chromosome 12, leaving 383.
# Pruned 2591 variants from chromosome 13, leaving 347.
# Pruned 5094 variants from chromosome 14, leaving 450.
# Pruned 3653 variants from chromosome 15, leaving 321.
# Pruned 2549 variants from chromosome 16, leaving 345.
# Pruned 2596 variants from chromosome 17, leaving 376.
# Pruned 2877 variants from chromosome 18, leaving 400.
# Pruning complete.  61402 of 68814 variants removed. (7412)
```

Settled on 4 subsets based on LD-pruning `plink --indep-pairwise` looking at 1 Mb windows in steps of 50 Kb pruning using the $r^2$ thresholds `0.9` (~20K), `0.8` (~15K), `0.6` (~10K), `0.5` (~7K). Skipping `0.7` because 12K just not diff enough from 10 or 15K for me.

Used plink-output list below to prune downstream input files. 

Subsetting the unfiltered version of the `haps` and `dosages` (e.g. `dosages_IITA_2021May13.rds`) can be done on read-in rather than creating new, smaller files on disk.

1. `grms` (kinship matrices) will need to be remade and pre-stored on disk.
2. Genetic map needs to interpolated that includes _all_ possible markers (in [previous version](https://wolfemd.github.io/implementGMSinCassava/04-PreprocessDataFiles.html#Genetic_Map), only did so for the ~34K "filtered" variants) and a `recombFreqMat` needs to be created and stored on disk based on that. After that the `recombFreqMat` can _also_ be subset on loading.

```{r grms, eval=F}
library(tidyverse); library(magrittr); 
dosages<-readRDS(file=here::here("data","dosages_IITA_2021May13.rds"))
source(here::here("code","gmsFunctions.R"))

snpsets<-tibble(Filter=c(5,6,8,9)) %>% 
  mutate(snps2keep=map(Filter,~read.table(here::here("output",
                                                     paste0("samples2keep_IITA_MAFpt01_prune1Mb_50kb_pt",.,".prune.in")),
                                          header = F, stringsAsFactors = F)))
snpsets %<>% 
  mutate(snps2keep=map(snps2keep,function(snps2keep,...){
    tokeep<-tibble(FULL_SNP_ID=colnames(dosages)) %>% 
      separate(FULL_SNP_ID,c("Chr","Pos","Ref","Alt"),remove = F) %>% 
      mutate(SNP_ID=paste0("S",Chr,"_",Pos)) %>% 
      filter(SNP_ID %in% snps2keep$V1)
    return(tokeep)}),
    Filter=paste0("1Mb_50kb_pt",Filter))

RhpcBLASctl::blas_set_num_threads(56)
snpsets %<>% 
  mutate(Amat=map(snps2keep,~kinship(dosages[,.$FULL_SNP_ID],type="add")),
         Dmat=map(snps2keep,~kinship(dosages[,.$FULL_SNP_ID],type="domGenotypic")))
library(qs)
qsave(snpsets,file=here::here("output","kinships_LDpruningSeries_2021Aug02.qs"),nthreads = 20)
```
Next interpolate a genetic map for _all_ 68K markers for max utility. As mentioned above, need it to make a `recombFreqMat`. 
```{r}
genmap<-read.table(here::here("data/CassavaGeneticMap",
                              "cassava_cM_pred.v6.allchr.txt"),
           header = F, stringsAsFactors = F,sep=';') %>% 
  rename(SNP_ID=V1,Pos=V2,cM=V3) %>% 
  as_tibble

snps_genmap<-tibble(DoseSNP_ID=colnames(dosages)) %>% 
  separate(DoseSNP_ID,c("Chr","Pos","Ref","Alt"),remove = F) %>% 
  mutate(SNP_ID=paste0("S",Chr,"_",Pos)) %>% 
  left_join(genmap %>% mutate(across(everything(),as.character)))
```
```{r}
interpolate_genmap<-function(data){
  # for each chromosome map
  # find and _decrements_ in the genetic map distance
  # fix them to the cumulative max to force map to be only increasing
  # fit a spline for each chromosome
  # Use it to predict values for positions not previously on the map
  # fix them AGAIN (in case) to the cumulative max, forcing map to only increase
  data_forspline<-data %>% 
    filter(!is.na(cM)) %>% 
    mutate(cumMax=cummax(cM),
           cumIncrement=cM-cumMax) %>% 
    filter(cumIncrement>=0) %>% 
    select(-cumMax,-cumIncrement)
  
  spline<-data_forspline %$% smooth.spline(x=Pos,y=cM,spar = 0.75)
  
  splinemap<-predict(spline,x = data$Pos) %>% 
    as_tibble(.) %>% 
    rename(Pos=x,cM=y) %>% 
    mutate(cumMax=cummax(cM),
           cumIncrement=cM-cumMax) %>% 
    mutate(cM=cumMax) %>% 
    select(-cumMax,-cumIncrement)
  
  return(splinemap) 
}
```
```{r genmap}
splined_snps_genmap<-snps_genmap %>% 
  select(-cM) %>% 
  mutate(Pos=as.numeric(Pos)) %>% 
  left_join(snps_genmap %>% 
              mutate(across(c(Pos,cM),as.numeric)) %>% 
              arrange(Chr,Pos) %>% 
              nest(data=-Chr) %>% 
              mutate(data=map(data,interpolate_genmap)) %>% 
              unnest(data)) %>% 
  distinct
```
```{r}
all(splined_snps_genmap$DoseSNP_ID == colnames(dosages))
# [1] TRUE

saveRDS(splined_snps_genmap,file=here::here("data","genmap_2021Aug02.rds"))
```
```{r, fig.width=9}
readRDS(here::here("data","genmap_2021Aug02.rds")) %>% 
  mutate(Map="Spline") %>% 
   ggplot(.,aes(x=Pos/1000/1000,y=cM,color=Map),alpha=0.5,size=0.75) + 
  geom_point() + 
  theme_bw() + facet_wrap(~as.integer(Chr), scales='free_x')
```


Construct a matrix of recombination frequencies at loci for all study loci. 
Pre-compute 1-2c to save time predicting cross variance.
```{r recombFreqMat, eval=F}
rm(list=ls()); gc()
library(tidyverse); library(magrittr); library(qs)
RhpcBLASctl::blas_set_num_threads(106)
source(here::here("code","predCrossVar.R"))
genmap<-readRDS(file=here::here("data","genmap_2021Aug02.rds"))

m<-genmap$cM;
names(m)<-genmap$DoseSNP_ID
recombFreqMat<-1-(2*genmap2recombfreq(m,nChr = 18))
qsave(recombFreqMat,file=here::here("data","recombFreqMat_1minus2c_2021Aug02.qs"),nthreads = 56)
```


## Parent-wise CV of marker subsets
```{r}
require(tidyverse); require(magrittr); library(qs)

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

# LD-PRUNNED SNP SETS + GRMS
snpsets<-qread(here::here("output","kinships_LDpruningSeries_2021Aug02.qs"),
               nthreads = 56)

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

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

# 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_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) 
```

**Note:** using same seed and other params as original run with 34K SNPs, which took ~34hrs.

```{r parent-wise CV - filterLevel settings}
# cbsulm15 - 112 cores, 512 GB RAM - 2021 Aug 3 - 4:20pm
filterLevel<-"1Mb_50kb_pt5"
# [1] "Accuracies predicted. Took  4.42406 hrs total.Goodbye!"
# [1] "Time elapsed: 265.469 mins"

# cbsulm17 - 112 cores, 512 GB RAM - 2021 Aug 3 - " "
filterLevel<-"1Mb_50kb_pt6"
# [1] "Accuracies predicted. Took  5.28384 hrs total.Goodbye!"
# [1] "Time elapsed: 317.053 mins"

# cbsulm29 - 104 cores, 512 GB RAM - 2021 Aug 3 - " "
filterLevel<-"1Mb_50kb_pt8"
# [1] "Accuracies predicted. Took  7.14035 hrs total.Goodbye!"
# [1] "Time elapsed: 428.437 mins"

# cbsulm31 - 112 cores, 512 GB RAM - 2021 Aug 3 - " " 
filterLevel<-"1Mb_50kb_pt9"
# [1] "Accuracies predicted. Took  10.65117 hrs total.Goodbye!"
# [1] "Time elapsed: 639.087 mins"

```

```{r code to run parent-wise CV for ea filter level}
snps2keep<-snpsets %>% 
  filter(Filter==filterLevel) %>% 
  select(snps2keep) %>% 
  unnest(snps2keep)
grms<-list(A=snpsets %>% filter(Filter==filterLevel) %$% Amat[[1]], 
           D=snpsets %>% filter(Filter==filterLevel) %$% Dmat[[1]])
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]
cvDirDom<-runParentWiseCrossVal(nrepeats=5,nfolds=5,seed=84,
                                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(cvDirDom,
        file = here::here("output",
                          paste0("cvDirDom_",filterLevel,"_predAccuracy.rds")))
endtime<-proc.time()[3]; print(paste0("Time elapsed: ",round((endtime-starttime)/60,3)," mins"))

```


## Standard CV of marker subsets
```{r}
require(tidyverse); require(magrittr); library(qs)

# 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) 
# LD-PRUNNED SNP SETS + GRMS
snpsets<-qread(here::here("output","kinships_LDpruningSeries_2021Aug02.qs"),
               nthreads = 56)

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

```

```{r standard CV - filterLevel settings}
# running...
# cbsulm18 - 88 cores, 512 GB RAM - 2021 Aug 3 - 1pm
filterLevel<-"1Mb_50kb_pt5" # [1] "Time elapsed: 131.09 mins"
# cbsulm18 - 88 cores, 512 GB RAM - 2021 Aug 3 - 4pm
filterLevel<-"1Mb_50kb_pt6" # [1] "Time elapsed: 131.54 mins"
# cbsulm18 - 88 cores, 512 GB RAM - 2021 Aug 3 - 7:15pm
filterLevel<-"1Mb_50kb_pt8" # [1] "Time elapsed: 136.943 mins"
# cbsulm18 - 88 cores, 512 GB RAM - 2021 Aug 3 - 10:20pm
filterLevel<-"1Mb_50kb_pt9"
```

```{r code to runCrossVal for ea filter level}
snps2keep<-snpsets %>% 
  filter(Filter==filterLevel) %>% 
  select(snps2keep) %>% 
  unnest(snps2keep)
grms<-list(A=snpsets %>% filter(Filter==filterLevel) %$% Amat[[1]], 
           D=snpsets %>% filter(Filter==filterLevel) %$% Dmat[[1]])
dosages<-dosages[,snps2keep$FULL_SNP_ID]
rm(snpsets); gc()

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

```


# [TO DO] PHG parent-wise cross-validation

# Next step / Results

6. [Genomic predictions](06-GenomicPredictions.html): 
  - 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](07-Results.html): Home for plots and summary tables.



