---
title: "Genomic predictions"
author: "Marnin Wolfe"
date: "2021-July-08"
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

5.  [Parent-wise cross-validation](05-CrossValidation.html): Compute parent-wise cross-validation folds using the validated pedigree. Fit models to get marker effects and make subsequent predictions of cross means and (co)variances.

# Current steps

1.  **Genomic prediction of clone GEBV/GETGV.** Fit GBLUP model, using genotypic add-dom partition. NEW: modelType="DirDom", include genome-wide inbreeding effect in GEBV/GETGV predictions *after* backsolving SNP effects. For all models, extract GBLUPs and backsolve SNP effects for use in cross usefulness predictions (mean+variance predictions). **ALSO NEW**: selection index predictions.

2.  **Genomic prediction of cross** $UC_{parent}$ and $UC_{variety}$. Rank potential parents on SI. Predict all possible crosses of some portion of best parents.

# Genomic prediction of clone GEBV/GETGV

Operate R 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; 
# Project directory, so R will use as working dir.
cd /home/mw489/implementGMSinCassava/;
# 3) Start R
R
```

Load input data
```{r load input data for GP, eval=F}
require(tidyverse); require(magrittr); 
# 5 threads per Rsession for matrix math (openblas)
RhpcBLASctl::blas_set_num_threads(5)
# GENOMIC MATE SELECTION FUNCTIONS
source(here::here("code","gmsFunctions.R")) 
source(here::here("code","predCrossVar.R")) 

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

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

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

Run the "AD" and "DirDom" modelTypes built into `runGenomicPredictions()`. 

Output contains both GBLUPs for selection of indivuals _and_ SNP effects to use as input for prediction of cross usefulness and subsequent mate selection.

```{r run genomic predictions}
start<-proc.time()[3]
gpreds_ad<-runGenomicPredictions(modelType="AD",selInd=TRUE, SIwts=SIwts,
                                 getMarkEffs=TRUE,
                                 returnPEV=FALSE,
                                 blups=blups,grms=grms,dosages=dosages,
                                 ncores=20,nBLASthreads=5)
runtime<-proc.time()[3]-start; runtime/60
saveRDS(gpreds_ad,file = here::here("output","genomicPredictions_ModelAD.rds"))

gpreds_dirdom<-runGenomicPredictions(modelType="DirDom",selInd=TRUE, SIwts=SIwts,
                                     getMarkEffs=TRUE,
                                     returnPEV=FALSE,
                                     blups=blups,grms=grms,dosages=dosages,
                                     ncores=20,nBLASthreads=5)
runtime<-proc.time()[3]-start; runtime/60
saveRDS(gpreds_dirdom,file = here::here("output","genomicPredictions_ModelDirDom.rds"))

```

# Genomic mate predictions: cross usefulness

1. Use GBLUPs on SELIND to choose some to fraction of the clones.
2. For those selected parents, predict the SELIND usefulness for all pairwise matings.

```{bash, 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 load inputs for predictCrosses, eval=F}
require(tidyverse); require(magrittr); 
# 5 threads per Rsession for matrix math (openblas)
#RhpcBLASctl::blas_set_num_threads(5)
# GENOMIC MATE SELECTION FUNCTIONS
source(here::here("code","gmsFunctions.R")) 
source(here::here("code","predCrossVar.R")) 

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

# SNP EFFECTS - FROM PREVIOUS STEP 
### Two models AD and DirDom
gpreds_ad<-readRDS(file = here::here("output","genomicPredictions_ModelAD.rds"))
gpreds_dirdom<-readRDS(file = here::here("output","genomicPredictions_ModelDirDom.rds"))
### Cor between SELIND GEBV and GETGV between models?
# left_join(gpreds_ad$gblups[[1]] %>% select(GID,predOf,SELIND) %>% rename(SELINDad=SELIND),
#           gpreds_dirdom$gblups[[1]] %>% select(GID,predOf,SELIND) %>% rename(SELINDdirdom=SELIND)) %>% 
#   group_by(predOf) %>% summarize(SELIND_corModels=cor(SELINDad,SELINDdirdom))
# predOf SELIND_corModels
# GEBV	0.9890091			
# GETGV	0.9373170	

# SELECTION OF CROSSES-TO-BE-PREDICTED
nParentsToSelect<-100
union_bestGEBVandGETGVdirdom<-union(gpreds_dirdom$gblups[[1]] %>% 
                                         filter(predOf=="GEBV") %>% 
                                         arrange(desc(SELIND)) %>% 
                                         slice(1:nParentsToSelect) %$% GID,
                                       gpreds_dirdom$gblups[[1]] %>% 
                                         filter(predOf=="GETGV") %>% 
                                         arrange(desc(SELIND)) %>% 
                                         slice(1:nParentsToSelect) %$% GID)
union_bestGEBVandGETGVad<-union(gpreds_ad$gblups[[1]] %>% 
                                         filter(predOf=="GEBV") %>% 
                                         arrange(desc(SELIND)) %>% 
                                         slice(1:nParentsToSelect) %$% GID,
                                       gpreds_ad$gblups[[1]] %>% 
                                         filter(predOf=="GETGV") %>% 
                                         arrange(desc(SELIND)) %>% 
                                         slice(1:nParentsToSelect) %$% GID)
union_bestGEBVandGETGV<-union(union_bestGEBVandGETGVdirdom,
                              union_bestGEBVandGETGVad)
length(union_bestGEBVandGETGV) 
# [1] 121 parents in top nParentsToSelect on SELIND for GEBV/GETGV - DirDom/AD model.

CrossesToPredict<-crosses2predict(union_bestGEBVandGETGV)
nrow(CrossesToPredict)
# [1] 7381 
```

Run the "AD" and "DirDom" modelTypes.  

cbsulm15 - July 25 - 7:30pm
```{r predict UC - top 121 parents - model DirDom}
start<-proc.time()[3]
crossPreds_dirdom<-predictCrosses(modelType="DirDom",stdSelInt = 2.67, 
                                  selInd=TRUE, SIwts=SIwts,
                                  CrossesToPredict=CrossesToPredict,
                                  snpeffs=gpreds_dirdom$genomicPredOut[[1]],
                                  dosages=dosages,
                                  haploMat=haploMat,recombFreqMat=recombFreqMat,
                                  ncores=20,nBLASthreads=5)
runtime<-proc.time()[3]-start; runtime/60
saveRDS(crossPreds_dirdom,file = here::here("output","genomicMatePredictions_top121parents_ModelDirDom.rds"))
# elapsed
#  2195.6
# # A tibble: 1 x 2
#   tidyPreds              rawPreds        
#   <list>                 <list>          
# 1 <tibble [177,144 × 9]> <named list [2]>

# 36.6 hrs for 7381 crosses. Or about 5 hours per 1000 crosses
```

cbsulm17 - July 25 at 7:30pm
```{r predict UC - top 121 parents - model AD}
start<-proc.time()[3]
crossPreds_ad<-predictCrosses(modelType="AD",stdSelInt = 2.67, 
                              selInd=TRUE, SIwts=SIwts,
                              CrossesToPredict=CrossesToPredict,
                              snpeffs=gpreds_ad$genomicPredOut[[1]],
                              dosages=dosages,
                              haploMat=haploMat,recombFreqMat=recombFreqMat,
                              ncores=20,nBLASthreads=5)
runtime<-proc.time()[3]-start; runtime/60
saveRDS(crossPreds_ad,file = here::here("output","genomicMatePredictions_top121parents_ModelAD.rds"))
# elapsed
# 1511.56
# # A tibble: 1 x 2
#   tidyPreds              rawPreds        
#   <list>                 <list>          
# 1 <tibble [177,144 × 9]> <named list [2]>

# 1511.56 or 25.33 hrs for 7381 crosses 

# About 3.4 hrs per 1000 crosses?
```

# Write tidy breeder-friendly predictions to disk

Add genetic groups and tidy format
```{r output tidy csv of gebv/getgv}
library(tidyverse);

gpreds_ad<-readRDS(file = here::here("output","genomicPredictions_ModelAD.rds")) 
gpreds_ad$gblups[[1]] %>% 
  mutate(GeneticGroup=case_when(grepl("2013_|TMS13",GID)~"C1",
                                grepl("TMS14",GID)~"C2",
                                grepl("TMS15",GID)~"C3",
                                grepl("TMS18",GID)~"C4",
                                !grepl("2013_|TMS13|TMS14|TMS15|TMS18",GID)~"PreGS"),
         GeneticGroup=factor(GeneticGroup,levels = c("PreGS","C1","C2","C3","C4"))) %>% 
  relocate(GeneticGroup,.after = "predOf") %>% 
  arrange(predOf,desc(SELIND)) %>% 
  write.csv(.,file = here::here("output","genomicPredictions_ModelAD.csv"),
            row.names = F)

gpreds_dirdom<-readRDS(file = here::here("output","genomicPredictions_ModelDirDom.rds"))
gpreds_dirdom$gblups[[1]] %>% 
  mutate(GeneticGroup=case_when(grepl("2013_|TMS13",GID)~"C1",
                                grepl("TMS14",GID)~"C2",
                                grepl("TMS15",GID)~"C3",
                                grepl("TMS18",GID)~"C4",
                                !grepl("2013_|TMS13|TMS14|TMS15|TMS18",GID)~"PreGS"),
         GeneticGroup=factor(GeneticGroup,levels = c("PreGS","C1","C2","C3","C4"))) %>% 
  relocate(GeneticGroup,.after = "predOf") %>% 
  arrange(predOf,desc(SELIND)) %>% 
  write.csv(.,file = here::here("output","genomicPredictions_ModelDirDom.csv"),
            row.names = F)
```

```{r output tidy csv of usefulness predictions}
crossPreds_dirdom<-readRDS(file = here::here("output","genomicMatePredictions_top121parents_ModelDirDom.rds"))
crossPreds_ad<-readRDS(file = here::here("output","genomicMatePredictions_top121parents_ModelAD.rds"))

crossPreds_ad$tidyPreds[[1]] %>% 
  mutate(sireGroup=case_when(grepl("2013_|TMS13",sireID)~"C1",
                             grepl("TMS14",sireID)~"C2",
                             grepl("TMS15",sireID)~"C3",
                             grepl("TMS18",sireID)~"C4",
                             !grepl("2013_|TMS13|TMS14|TMS15|TMS18",sireID)~"PreGS"),
         damGroup=case_when(grepl("2013_|TMS13",damID)~"C1",
                            grepl("TMS14",damID)~"C2",
                            grepl("TMS15",damID)~"C3",
                            grepl("TMS18",damID)~"C4",
                            !grepl("2013_|TMS13|TMS14|TMS15|TMS18",damID)~"PreGS"),
         CrossGroup=paste0(sireGroup,"x",damGroup)) %>% 
  relocate(contains("Group"),.before = "Nsegsnps") %>% 
  write.csv(.,file = here::here("output","genomicMatePredictions_top121parents_ModelAD.csv"),
            row.names = F)

crossPreds_dirdom$tidyPreds[[1]] %>% 
  mutate(sireGroup=case_when(grepl("2013_|TMS13",sireID)~"C1",
                             grepl("TMS14",sireID)~"C2",
                             grepl("TMS15",sireID)~"C3",
                             grepl("TMS18",sireID)~"C4",
                             !grepl("2013_|TMS13|TMS14|TMS15|TMS18",sireID)~"PreGS"),
         damGroup=case_when(grepl("2013_|TMS13",damID)~"C1",
                            grepl("TMS14",damID)~"C2",
                            grepl("TMS15",damID)~"C3",
                            grepl("TMS18",damID)~"C4",
                            !grepl("2013_|TMS13|TMS14|TMS15|TMS18",damID)~"PreGS"),
         CrossGroup=paste0(sireGroup,"x",damGroup)) %>% 
  relocate(contains("Group"),.before = "Nsegsnps") %>% 
  write.csv(.,file = here::here("output","genomicMatePredictions_top121parents_ModelDirDom.csv"),
            row.names = F)
```

Write CSV's

```{r}
# ## Format and write GEBV
# predModelA %>% 
#   select(Trait,genomicPredOut) %>% 
#   unnest(genomicPredOut) %>% 
#   select(-varcomps) %>% 
#   unnest(gblups) %>% 
#   select(-GETGV,-contains("PEV")) %>%
#   spread(Trait,GEBV) %>% 
#   mutate(Group=case_when(GID %in% nrTP ~ "nrTP",
#                          GID %in% c1a ~ "C1a",
#                          GID %in% c1b ~ "C1b",
#                          GID %in% c2a ~ "C2a",
#                          GID %in% c2b ~ "C2b",
#                          GID %in% c3a ~ "C3a",
#                          GID %in% c3b ~ "C3b")) %>% 
#   select(Group,GID,any_of(traits)) %>% 
#   arrange(desc(Group)) %>% 
#   write.csv(., file = here::here("output","GEBV_NRCRI_ModelA_2021May03.csv"), row.names = F)
# 
# ## Format and write GETGV
# predModelADE %>% 
#   select(Trait,genomicPredOut) %>% 
#   unnest(genomicPredOut) %>% 
#   select(-varcomps) %>% 
#   unnest(gblups) %>% 
#   select(GID,Trait,GETGV) %>% 
#   spread(Trait,GETGV) %>% 
#   mutate(Group=case_when(GID %in% nrTP ~ "nrTP",
#                          GID %in% c1a ~ "C1a",
#                          GID %in% c1b ~ "C1b",
#                          GID %in% c2a ~ "C2a",
#                          GID %in% c2b ~ "C2b",
#                          GID %in% c3a ~ "C3a",
#                          GID %in% c3b ~ "C3b")) %>% 
#   select(Group,GID,any_of(traits)) %>% 
#   arrange(desc(Group)) %>% 
#   write.csv(., file = here::here("output","GETGV_NRCRI_ModelADE_2021May03.csv"), row.names = F)
# 
# ### Make a unified "tidy" long-form: 
# predModelA %>% 
#   select(Trait,genomicPredOut) %>% 
#   unnest(genomicPredOut) %>% 
#   select(-varcomps) %>% 
#   unnest(gblups) %>% 
#   select(-GETGV) %>% 
#   full_join(predModelADE %>% 
#               select(Trait,genomicPredOut) %>% 
#               unnest(genomicPredOut) %>% 
#               select(-varcomps) %>% 
#               unnest(gblups) %>% 
#               rename(GEBV_modelADE=GEBV,
#                      PEV_modelADE=PEVa) %>% 
#               select(-genomicPredOut)) %>% 
#   mutate(Group=case_when(GID %in% nrTP ~ "nrTP",
#                          GID %in% c1a ~ "C1a",
#                          GID %in% c1b ~ "C1b",
#                          GID %in% c2a ~ "C2a",
#                          GID %in% c2b ~ "C2b",
#                          GID %in% c3a ~ "C3a",
#                          GID %in% c3b ~ "C3b")) %>% 
#   relocate(Group,.before = GID) %>% 
#   write.csv(., file = here::here("output","genomicPredictions_NRCRI_2021May03.csv"), row.names = F)
```




# Results

See [Results](07-Results.html): Home for plots and summary tables.
