---
title: "Pre-process data files for downstream cross-validation"
author: "Marnin Wolfe"
date: "2021-May-13"
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, 
                      tidy='styler', tidy.opts=list(strict=FALSE,width.cutoff=100), highlight=TRUE)
```

# Previous step

3. [Validate the pedigree obtained from cassavabase](03-validatePedigree.html): Before setting up a cross-validation scheme for  predictions that depend on a correct pedigree, add a basic verification step to the pipeline. Not trying to fill unknown or otherwise correct the pedigree. Assess evidence that relationship is correct, remove if incorrect.

# Haplotype matrix from phased VCF

Extract haps from VCF with `bcftools`

```{r, eval=F}
library(tidyverse); library(magrittr)
pathIn<-"/home/jj332_cas/marnin/implementGMSinCassava/data/"
pathOut<-pathIn
vcfName<-"AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719"
system(paste0("bcftools convert --hapsample ",
              pathOut,vcfName," ",
              pathIn,vcfName,".vcf.gz "))
```
Read haps to R
```{r, eval=F}
library(data.table)
haps<-fread(paste0(pathIn,vcfName,".hap.gz"),
            stringsAsFactors = F,header = F) %>% 
  as.data.frame
sampleids<-fread(paste0(pathIn,vcfName,".sample"),
                 stringsAsFactors = F,header = F,skip = 2) %>% 
  as.data.frame
```

**Extract needed GIDs from BLUPs and pedigree:** Subset to: (1) genotyped-plus-phenotyped and/or (2) in verified pedigree.

```{r, eval=F}
blups<-readRDS(file=here::here("output",
                               "IITA_blupsForModelTraining_twostage_asreml_2021May10.rds"))
blups %>% 
  select(Trait,blups) %>% 
  unnest(blups) %>% 
  distinct(GID) %$% GID -> gidWithBLUPs

genotypedWithBLUPs<-gidWithBLUPs[gidWithBLUPs %in% sampleids$V1]
length(genotypedWithBLUPs) # 7960

ped<-read.table(here::here("output","verified_ped.txt"),
                header = T, stringsAsFactors = F)

pednames<-union(ped$FullSampleName,
                union(ped$SireID,ped$DamID))
length(pednames) # 4384

samples2keep<-union(genotypedWithBLUPs,pednames)
length(samples2keep) # 8013

# write a sample list to disk for downstream purposes
# format suitable for subsetting with --keep in plink
write.table(tibble(FID=0,IID=samples2keep),
            file=here::here("output","samples2keep_IITA_2021May13.txt"),
            row.names = F, col.names = F, quote = F)
```

Add sample ID's
```{r, eval=F}
hapids<-sampleids %>% 
  select(V1,V2) %>% 
  mutate(SampleIndex=1:nrow(.)) %>% 
  rename(HapA=V1,HapB=V2) %>% 
  pivot_longer(cols=c(HapA,HapB),
               names_to = "Haplo",values_to = "SampleID") %>% 
  mutate(HapID=paste0(SampleID,"_",Haplo)) %>% 
  arrange(SampleIndex)
colnames(haps)<-c("Chr","HAP_ID","Pos","REF","ALT",hapids$HapID)
```

Subset haps
```{r, eval=F}
hapids2keep<-hapids %>% filter(SampleID %in% samples2keep)
hapids2keep$HapID
dim(haps) # [1] 68814 43717
haps<-haps[,c("Chr","HAP_ID","Pos","REF","ALT",hapids2keep$HapID)]
dim(haps) # [1] 68814 16031
```

Format, transpose, convert to matrix and save!
```{r, eval=F}
haps %<>% 
  mutate(HAP_ID=gsub(":","_",HAP_ID)) %>% 
  column_to_rownames(var = "HAP_ID") %>% 
  select(-Chr,-Pos,-REF,-ALT)
haps %<>% t(.) %>% as.matrix(.)
saveRDS(haps,file=here::here("data","haps_IITA_2021May13.rds")
```

# Make dosages from haps

To ensure consistency in allele counting, create dosage from haps manually.

```{r, eval=F}
dosages<-haps %>%
  as.data.frame(.) %>% 
  rownames_to_column(var = "GID") %>% 
  separate(GID,c("SampleID","Haplo"),"_Hap",remove = T) %>% 
  select(-Haplo) %>% 
  group_by(SampleID) %>% 
  summarise(across(everything(),~sum(.))) %>% 
  ungroup() %>% 
  column_to_rownames(var = "SampleID") %>% 
  as.matrix
saveRDS(dosages,file=here::here("data","dosages_IITA_2021May13.rds"))
# > dim(dosages)
# [1]  8013 68814
```

# Variant filters

**Apply a MAF filter and lightly LD prune:** The number of markers in the "raw" dataset (~68K) is ~3X the number used in the mate selection paper and I think more than is necessary. There is a burden incurred because we have to compute and store in memory (and on disk) $N_{snp} \times N_{snp}$ recombination frequency matrices. 

```{r}
# library(tidyverse); library(magrittr)
# pathIn<-"/home/jj332_cas/marnin/implementGMSinCassava/data/"
# pathOut<-pathIn
# vcfName<-"AllChrom_RefPanelAndGSprogeny_ReadyForGP_72719"
# 
# write.table(tibble(FID=0,IID=samples2keep),
#             file=here::here("output","samples2keep_IITA_2021May13.txt"),
#             row.names = F, col.names = F, quote = F)
# 
# ped2check<-read.table(file=here::here("output","ped2genos.txt"),
#                       header = F, stringsAsFactors = F)
# 
# # pednames<-union(ped2check$V1,union(ped2check$V2,ped2check$V3)) %>% 
# #   tibble(FID=0,IID=.)
# # write.table(pednames,file=here::here("output","pednames2keep.txt"), 
# #             row.names = F, col.names = F, quote = F)
```

```{plink, eval=F}
cd /home/jj332_cas/marnin/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 \
  --maf 0.01 \
  --indep-pairwise 50 25 0.98 \
  --out output/samples2keep_IITA_MAFpt01_prune50_25_pt98;
```

Used plink to output a list of pruned SNPs. 

Next, subset the columns of `haps` and `dosages` in R.

```{r, eval=F}
library(tidyverse); library(magrittr); 
haps<-readRDS(file=here::here("data","haps_IITA_2021May13.rds"))
dosages<-readRDS(file=here::here("data","dosages_IITA_2021May13.rds"))
snps2keep<-read.table(here::here("output",
                      "samples2keep_IITA_MAFpt01_prune50_25_pt98.prune.in"),
           header = F, stringsAsFactors = F)
snps2keep<-tibble(HapSNP_ID=colnames(haps)) %>% 
  separate(HapSNP_ID,c("Chr","Pos","Ref","Alt"),remove = F) %>% 
  mutate(SNP_ID=paste0("S",Chr,"_",Pos)) %>% 
  filter(SNP_ID %in% snps2keep$V1)

haps<-haps[,snps2keep$HapSNP_ID]
dosages<-dosages[,snps2keep$HapSNP_ID]

# dim(haps); dim(dosages); haps[1:5,1:10]

saveRDS(haps,file=here::here("data","haps_IITA_filtered_2021May13.rds"))
saveRDS(dosages,file=here::here("data","dosages_IITA_filtered_2021May13.rds"))
```

# Make Add and Dom GRMs from dosages
```{bash, eval=F}
# activate multithread OpenBLAS for fast matrix algebra
export OMP_NUM_THREADS=56
```
```{r, eval=F}
dosages<-readRDS(file=here::here("data","dosages_IITA_filtered_2021May13.rds"))
A<-predCrossVar::kinship(dosages,type="add")
D<-predCrossVar::kinship(dosages,type="dom")
saveRDS(A,file=here::here("output","kinship_A_IITA_2021May13.rds"))
saveRDS(D,file=here::here("output","kinship_D_IITA_2021May13.rds"))
```

```{bash, eval=F}
cd /home/mw489/implementGMSinCassava/;
screen; 
singularity shell rocker.sif; R
```

```{r, eval=F}
dosages<-readRDS(file=here::here("data","dosages_IITA_filtered_2021May13.rds"))
source(here::here("code","gsFunctions.R"))
RhpcBLASctl::blas_set_num_threads(56)
D<-kinship(dosages,type="domGenotypic")
saveRDS(D,file=here::here("output","kinship_domGenotypic_IITA_2021July5.rds"))
```


# Genetic Map

```{bash,eval = FALSE}
cp -r /home/jj332_cas/CassavaGenotypeData/CassavaGeneticMap /home/jj332_cas/marnin/implementGMSinCassava/data/
```

```{bash, eval=F}
# activate multithread OpenBLAS for fast matrix algebra
export OMP_NUM_THREADS=56
```

[Creating the map used for Beagle-imputation in 2019](https://wolfemd.github.io/NaCRRI_2020GS/Imputation_EastAfrica_StageI_82819.html#Genetic_Map_for_Beagle_-_V2): In 2019, I obtained a ICGMC-derived genetic map, I _think_ from Guillaume Bauchet and used it to create a map I've been using for imputation, which has 25K markers (Beagle interpolates the map to the markers genotyped in the panel).

*However,* the recombination frequency matrix and thus cross-variance predictions needs to have all positions for which we have marker effects. It means I have to interpolate a map from the original file `cassava_cM_pred.v6.allchr.txt`. See below:

```{r}
library(tidyverse); library(magrittr)
dosages<-readRDS(file=here::here("data","dosages_IITA_filtered_2021May13.rds"))
# genmap<-tibble(Chr=1:18) %>% 
#   mutate(geneticMap=map(Chr,~read.table(here::here("data/CassavaGeneticMap",
#                                                    paste0("chr",.,"_cassava_cM_pred.v6_91019.map")),
#                                         header = F, stringsAsFactors = F)))

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)))
# snps_genmap %>% 
#   ggplot(.,aes(x=as.integer(Pos),y=as.numeric(cM))) + 
#   geom_point() + 
#   theme_bw() + 
#   facet_wrap(~Chr)
```
```{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}
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(-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_2021May13.rds"))
```


```{r, fig.width=12}
splined_snps_genmap %>% 
  mutate(Map="Spline") %>% 
  bind_rows(snps_genmap %>% 
              mutate(across(c(Pos,cM),as.numeric)) %>% 
              arrange(Chr,Pos) %>% mutate(Map="Data")) %>% 
  ggplot(.,aes(x=Pos,y=cM,color=Map),alpha=0.5,size=0.75) + 
  geom_point() + 
  theme_bw() + facet_wrap(~as.integer(Chr), scales='free_x')
```

# Recomb. freq. matrix

Construct a matrix of recombination frequencies at loci for all study loci. 
Pre-compute 1-2c to save time predicting cross variance.
```{r, eval=F}
library(predCrossVar)
genmap<-readRDS(file=here::here("data","genmap_2021May13.rds"))
m<-genmap$cM;
names(m)<-genmap$DoseSNP_ID
recombFreqMat<-1-(2*genmap2recombfreq(m,nChr = 18))
saveRDS(recombFreqMat,file=here::here("data","recombFreqMat_1minus2c_2021May13.rds"))
```


# Pick traits to cross-validate
```{r}
# This list from Dec. 2020 GeneticGain rate estimation
# These were what Ismail/IITA/BMGF wanted to see
# Will cross-validate these traits
traits<-c("logDYLD","logFYLD","logRTNO","logTOPYLD","MCMDS","DM","BCHROMO",
          "PLTHT","BRLVLS","BRNHT1","HI")

# Full trait list = 14:
## traits<-c("MCMDS","DM","PLTHT","BRNHT1","BRLVLS","HI",
##           "logDYLD", # <-- logDYLD now included.
##           "logFYLD","logTOPYLD","logRTNO","TCHART","LCHROMO","ACHROMO","BCHROMO")

```

# Pedigree

3. [Validate the pedigree obtained from cassavabase](03-validatePedigree.html): Before setting up a cross-validation scheme for  predictions that depend on a correct pedigree, add a basic verification step to the pipeline. Not trying to fill unknown or otherwise correct the pedigree. Assess evidence that relationship is correct, remove if incorrect.

```{r, eval=F}
ped<-read.table(here::here("output","verified_ped.txt"),
                header = T, stringsAsFactors = F)
```

# BLUPs

Select traits and data to be analyzed.
```{r, eval=F}
library(tidyverse); library(magrittr);
blups<-readRDS(file=here::here("output",
                               "IITA_blupsForModelTraining_twostage_asreml_2021May10.rds"))
dosages<-readRDS(file=here::here("data","dosages_IITA_filtered_2021May13.rds"))


blups %>% 
  select(Trait,blups) %>% 
  unnest(blups) %>% 
  distinct(GID) %$% GID -> gidWithBLUPs

genotypedWithBLUPs<-gidWithBLUPs[gidWithBLUPs %in% rownames(dosages)]
length(genotypedWithBLUPs) # 7960

blups %<>% 
  filter(Trait %in% traits) %>% 
  select(Trait,blups,varcomp) %>% 
  mutate(blups=map(blups,~filter(.,GID %in% genotypedWithBLUPs)))

saveRDS(blups,file=here::here("data","blups_forCrossVal.rds"))
```

# Index weights [get from Ismail]
```{r, eval=F}
# library(tidyverse); library(magrittr); library(predCrossVar); library(BGLR);
# blups<-readRDS(here::here("data","blups_forawcdata.rds")) %>% 
#   select(Trait,blups) %>% # BLUPs long-->wide for multivar analysis
#   unnest(blups) %>% 
#   select(Trait,germplasmName,drgBLUP) %>% 
#   spread(Trait,drgBLUP)
# 
# indices<-blups %>% 
#   summarize_if(is.numeric,sd, na.rm=T) %>% 
#   pivot_longer(cols = everything(), names_to = "Trait", values_to = "blupSD") %>% 
#   left_join(tibble(Trait=c("DM","logFYLD","MCMDS","TCHART"), 
#                    stdSI_unscaled=c(5, 10, -10, -5),
#                    biofortSI_unscaled=c(10, 5, -5,10))) %>% 
#   mutate(stdSI=stdSI_unscaled/blupSD,
#          biofortSI=biofortSI_unscaled/blupSD)
# indices %>% mutate_if(is.numeric,~round(.,2))
# saveRDS(indices,file=here::here("data","selection_index_weights_4traits.rds"))
```


# Next 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.
