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

```{bash Renv, 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
```

Load input data

```{r inputs for GP}
# 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) 
```

Run the "DirDom" modelTypes built into `runGenomicPredictions()`. Output will contain both GBLUPs for selection of clones _and_ SNP effects to use as input for prediction of cross usefulness and subsequent mate selection. 

Get effects for: **full_set** (~31K SNPs), **medium_set** (~13K, LD-pruned) and **reduced_set** (~9K SNPs, LD-pruned).

```{r GP analysis chunks}
# cbsulm30 - 112 cores, 512 GB RAM - 2021 Aug 10 - 8:40am 
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: 24.066 mins"

# cbsulm30 - 112 cores, 512 GB RAM - 2021 Aug 10 - 8:40am
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: 22.431 mins"

# cbsulm20 - 88 cores, 512 GB RAM - 2021 Aug 10 - 2:50pm
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: 13.299 mins"

```
```{r run GP 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]
gpreds<-runGenomicPredictions(modelType="DirDom",selInd=TRUE, SIwts=SIwts,
                              getMarkEffs=TRUE,
                              returnPEV=FALSE,
                              blups=blups,grms=grms,dosages=dosages,
                              ncores=11,nBLASthreads=5)
saveRDS(gpreds,
        file = here::here("output",
                          paste0("genomicPredictions_",snpSet,"_2021Aug09.rds")))
endtime<-proc.time()[3]; print(paste0("Time elapsed: ",
                                      round((endtime-starttime)/60,3)," mins"))

```

# Genomic prediction of cross usefulness

1. Use GBLUPs on SELIND to choose top fraction of the clones. **Cross-reference with 'accessions_infield' list**.
2. For those selected parents, predict the SELIND usefulness for all pairwise matings.

```{bash Renv for crossPreds, 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 crossPred input}
# NEEDED LIBRARIES
require(tidyverse); require(magrittr); library(qs)
library(genomicMateSelectR)

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

# LIST OF ACCESSIONS LIKELY IN THE FIELD
accessions_infield<-readRDS(here::here("data",
                                       "accessions_possibly_infield_2021Aug10.rds"))

```

```{r cor of GBLUPs among snpsets, eval=F}
# gpreds_full<-readRDS(file = here::here("output","genomicPredictions_full_set_2021Aug09.rds"))
# gpreds_medium<-readRDS(file = here::here("output","genomicPredictions_medium_set_2021Aug09.rds"))
# gpreds_reduced<-readRDS(file = here::here("output","genomicPredictions_reduced_set_2021Aug09.rds"))
### Quick check that GBLUPs from full and medium_set are strongly correlated
### This will be an additional assurance that similar cross variances
### will be predicted by the reduced SNP model

### Cor between SELIND GEBV and GETGV between full_set and reduced_set of SNPs?
# left_join(gpreds_full$gblups[[1]] %>% select(GID,predOf,SELIND) %>% rename(SELIND_full=SELIND),
#           gpreds_reduced$gblups[[1]] %>% select(GID,predOf,SELIND) %>% rename(SELIND_reduced=SELIND)) %>%
#   group_by(predOf) %>% 
#   summarize(SELIND_corModels=cor(SELIND_full,SELIND_reduced))
# predOf  SELIND_corModels
# GEBV	0.0363334			
# GETGV	0.7809016	
## TERRIBLE!! :(

### Cor between SELIND GEBV and GETGV between full_set and medium_set of SNPs?
# left_join(gpreds_full$gblups[[1]] %>% select(GID,predOf,SELIND) %>% rename(SELIND_full=SELIND),
#           gpreds_medium$gblups[[1]] %>% select(GID,predOf,SELIND) %>% rename(SELIND_medium=SELIND)) %>%
#   group_by(predOf) %>% 
#   summarize(SELIND_corModels=cor(SELIND_full,SELIND_medium))
# predOf  SELIND_corModels
# GEBV	0.9901325			
# GETGV	0.9887326
## EXCELLENT!! :)
```

Choose the best parents for which to predict crosses. Use the GBLUPs from the **full_set** of SNPs.

Take the union of the top 300 clones on the SELIND in terms of GEBV and of GETGV. Probably they will be a very similar list.

```{r select best parents for crossPreds, eval=F}
# SELECT THE BEST PARENTS AS CROSSES-TO-BE-PREDICTED
nParentsToSelect<-300
gpreds_full<-readRDS(file = here::here("output","genomicPredictions_full_set_2021Aug09.rds"))
union_bestGEBVandGETGV<-union(gpreds_full$gblups[[1]] %>% 
                                filter(predOf=="GEBV") %>% 
                                arrange(desc(SELIND)) %>% 
                                slice(1:nParentsToSelect) %$% GID,
                              gpreds_full$gblups[[1]] %>% 
                                filter(predOf=="GETGV") %>% 
                                arrange(desc(SELIND)) %>% 
                                slice(1:nParentsToSelect) %$% GID)
rm(gpreds_full); 
length(union_bestGEBVandGETGV) 
# [1] 365 parents in top nParentsToSelect on SELIND for GEBV/GETGV 

# KEEP ONLY CANDIDATE PARENTS EXPECTED TO BE IN THE FIELD
table(union_bestGEBVandGETGV %in% accessions_infield$FullSampleName)
# FALSE  TRUE 
#   165   200
parentsToPredictCrosses<-union_bestGEBVandGETGV %>% 
  .[. %in% accessions_infield$FullSampleName]
CrossesToPredict<-crosses2predict(parentsToPredictCrosses)
nrow(CrossesToPredict)
# [1] 20100 possible crosses of 200 parents

saveRDS(parentsToPredictCrosses,
        file = here::here("output",
                          "parentsToPredictCrosses_2021Aug10.rds"))
saveRDS(CrossesToPredict,
        file = here::here("output",
                          "CrossesToPredict_2021Aug10.rds"))
```

Predict all pairwise crosses of those 200 parents.

```{r crossPred analysis chunks}
# cbsulm17 - 112 cores, 512 GB RAM - 2021 Aug 11 - 8:10am
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")))
predTheMeans<-TRUE; predTheVars<-FALSE
gpreds<-readRDS(file = here::here("output","genomicPredictions_full_set_2021Aug09.rds")); gc()
## takes ~2 minutes, predicting means only

# cbsulm - 112 cores, 512 GB RAM - 2021 Aug 10 - 8:25am
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")))
predTheMeans<-TRUE; predTheVars<-TRUE
gpreds<-readRDS(file = here::here("output","genomicPredictions_medium_set_2021Aug09.rds")); gc()

```

```{r run crossPreds 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()
snpSet; dim(dosages); predTheMeans; predTheVars

start<-proc.time()[3]
crossPreds<-predictCrosses(modelType="DirDom",stdSelInt = 2.0627128, 
                           selInd=TRUE, SIwts=SIwts,
                           CrossesToPredict=CrossesToPredict,
                           snpeffs=gpreds$genomicPredOut[[1]],
                           dosages=dosages,
                           haploMat=haploMat,recombFreqMat=recombFreqMat,
                           ncores=20,nBLASthreads=5,
                           predTheMeans = predTheMeans,
                           predTheVars = predTheVars)
runtime<-proc.time()[3]-start; runtime/60
saveRDS(crossPreds,file = here::here("output",
                                     paste0("genomicMatePredictions_",
                                            snpSet,"_2021Aug10.rds")))

```

# Write tidy breeder-friendly predictions to disk

Add genetic groups and cohort identifiers and tidy format
```{r output tidy csv of gebv/getgv}
library(tidyverse);
gpreds<-readRDS(file = here::here("output","genomicPredictions_full_set_2021Aug09.rds")) 
gpreds$gblups[[1]] %>% 
  mutate(GeneticGroup=case_when(grepl("2013_|TMS13",GID)~"C1",
                                grepl("TMS14",GID)~"C2",
                                grepl("TMS15",GID)~"C3",
                                grepl("TMS18",GID)~"C4",
                                grepl("TMS20",GID)~"C5",
                                grepl("TMS16",GID)~"TMS16",
                                grepl("TMS17",GID)~"TMS17",
                                grepl("TMS19",GID)~"TMS19",
                                !grepl("2013_|TMS13|TMS14|TMS15|TMS16|TMS17|TMS18|TMS19|TMS20",GID)~"PreGS"),
         #GeneticGroup=factor(GeneticGroup,levels = c("PreGS","C1","C2","C3","C4","C5","TMS16","TMS17","TMS19")),
         Cohort=case_when(grepl("2013_|TMS13",GID)~"TMS13",
                          grepl("TMS14",GID)~"TMS14",
                          grepl("TMS15",GID)~"TMS15",
                          grepl("TMS16",GID)~"TMS16",
                          grepl("TMS17",GID)~"TMS17",
                          grepl("TMS18",GID)~"TMS18",
                          grepl("TMS19",GID)~"TMS19",
                          grepl("TMS20",GID)~"TMS20",
                          !grepl("2013_|TMS13|TMS14|TMS15|TMS16|TMS17|TMS18|TMS19|TMS20",GID)~"PreGS")) %>% 
  relocate(GeneticGroup,.after = "predOf") %>% 
  relocate(Cohort,.after = "GeneticGroup") %>% 
  arrange(predOf,desc(SELIND)) %>% 
  write.csv(.,file = here::here("output","genomicPredictions_full_set_2021Aug09.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)
```





# Results

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