---
title: "Impute NRCRI DCas20_5440"
site: workflowr::wflow_site
date: "2020-October-8"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: inline
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
```

DArTseqLD (DCas20-5440). From Princess Onyyegbule on Sep 22, 2020: "These samples are from materials for inbreeding depression after one self-pollinated generation in two elite cassava varieties (TMS980581 and TMS070337). We want to assess the possibility of obtaining genetic gain by selecting transgressive individuals based on several productive traits, mostly high dry matter content."

Impute with RefPanelWA (W. Africa RefPanel)

# Copy data

Copy the imputation reference panel from 2019 to the `data/` folder. 

```{bash,eval = FALSE}
cp -r /home/jj332_cas/marnin/NRCRI_2020GS /workdir/mw489/
cp -r /home/jj332_cas/marnin/NRCRI_2020GS/data/Report-DCas20-5440 /workdir/mw489/NRCRI_2020GS/data/

cp -r /home/jj332_cas/CassavaGenotypeData/CassavaGeneticMap /workdir/mw489/NRCRI_2020GS/data/
cp /home/jj332_cas/CassavaGenotypeData/nextgenImputation2019/ImputationStageII_71219/chr*_ImputationReferencePanel_StageIIpartI_72219.vcf.gz /workdir/mw489/NRCRI_2020GS/data/
```

## With RefPanelWA

Impute with [Beagle V5.0](https://faculty.washington.edu/browning/beagle/b5_0.html).

Use the "imputation reference panel" dataset from 2019, e.g. `chr1_ImputationReferencePanel_StageIIpartI_72219.vcf.gz` as reference.

Used 1 large memory Cornell CBSU machine (e.g. [cbsulm16; 112 cores, 512 GB RAM](https://biohpc.cornell.edu/lab/hardware.aspx)), running 1 chromosome at a time.

R functions are stored in the `code/` sub-directory. Functions sourced from e.g. **imputationFunctions.R** are wrappers around e.g. Beagle, and other command line programs.

### Impute
```{r,eval = FALSE}
#library(tidyverse); library(magrittr);
source(here::here("code","imputationFunctions.R"))
targetVCFpath<-here::here("data/Report-DCas20-5440/") # location of the targetVCF
refVCFpath<-here::here("data/")
mapPath<-here::here("data/CassavaGeneticMap/")
outPath<-here::here("output/")
outSuffix<-"DCas20_5440"

purrr::map(1:18,~runBeagle5(targetVCF=paste0(targetVCFpath,"chr",.,"_DCas20_5440.vcf.gz"),
                            refVCF=paste0(refVCFpath,"chr",.,"_ImputationReferencePanel_StageIIpartI_72219.vcf.gz"),
                            mapFile=paste0(mapPath,"chr",.,"_cassava_cM_pred.v6_91019.map"),
                            outName=paste0(outPath,"chr",.,"_DCas20_5440_WA_REFimputed"),
                            nthreads=112))
```

Clean up Beagle log files after run. Move to sub-directory `output/BeagleLogs/`.
```{bash,eval = FALSE}
cd /workdir/mw489/NRCRI_2020GS/output/; 
mkdir BeagleLogs;
cp *_DCas20_5440_WA_REFimputed.log BeagleLogs/
cp -r BeagleLogs /home/jj332_cas/marnin/NRCRI_2020GS/output/
cp *_DCas20_5440_WA_REFimputed* /home/jj332_cas/marnin/NRCRI_2020GS/output/

```

### Post-impute filter

For now, the function will just do a fixed filter: AR2>0.75 (DR2>0.75 as of Beagle5.0), P_HWE>1e-20, MAF>0.005 [0.5%]. 

It can easily be modified in the future to include parameters to vary the filter specifications.

Input parameters
```{r,eval = FALSE}
#' @inPath path to input VCF-to-be-filtered, can be left null if path included in @inName . Must end in "/"
#' @inName name of input VCF file EXCLUDING file extension. Assumes .vcf.gz
#' @outPath path where filtered VCF and related are to be stored.Can be left null if path included in @outName . Must end in "/".
#' @outName name desired for output EXCLUDING extension. Output will be .vcf.gz 
```

Loop to filter all 18 VCF files in parallel
```{r,eval = FALSE}
inPath<-here::here("output/")
outPath<-here::here("output/")
source(here::here("code","imputationFunctions.R"))
require(furrr); options(mc.cores=ncores); plan(multiprocess)
future_map(1:18,~postImputeFilter(inPath=inPath,
                                  inName=paste0("chr",.,"_DCas20_5440_WA_REFimputed"),
                                  outPath=outPath,
                                  outName=paste0("chr",.,"_DCas20_5440_WA_REFimputedAndFiltered")))
```

Check what's left
```{r,eval = FALSE}
purrr::map(1:18,~system(paste0("zcat ",here::here("output/"),"chr",.,"_DCas20_5440_WA_REFimputedAndFiltered.vcf.gz | wc -l")))
# 6912
# 2961
# 2659
# 2317
# 2329
# 2258
# 1426
# 2229
# 2478
# 2207
# 2130
# 2290
# 1930
# 2744
# 2873
# 2015
# 1580
# 1772
```

```{bash, eval=F}
cd /workdir/mw489/NRCRI_2020GS/output/;
cp *_DCas20_5440_WA_REFimputedAndFiltered* /home/jj332_cas/marnin/NRCRI_2020GS/output/
```

### Formats for downstream analysis

```{r, eval=F}
library(tidyverse); library(magrittr);
# Make binary plink
pathIn<-"/home/jj332_cas/marnin/NRCRI_2020GS/output/"
require(furrr); options(mc.cores=ncores); plan(multiprocess)
future_map(1:18,~system(paste0("export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;",
                               "plink --vcf ",pathIn,"chr",.,
                               "_DCas20_5440_WA_REFimputedAndFiltered.vcf.gz ",
                               "--make-bed --const-fid ",
                               "--out ",pathIn,"chr",.,
                               "_DCas20_5440_WA_REFimputedAndFiltered")))
# Recode to dosage
future_map(1:18,~system(paste0("export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;",
                               "plink --bfile ",pathIn,"chr",.,
                               "_DCas20_5440_WA_REFimputedAndFiltered ",
                               "--recode A ",
                               "--out ",pathIn,"chr",.,
                               "_DCas20_5440_WA_REFimputedAndFiltered")))

# Genome-wide dosage (for use in R)
snps<-future_map(1:18,~read.table(paste0(pathIn,"chr",.,"_DCas20_5440_WA_REFimputedAndFiltered.raw"), stringsAsFactor=F, header = T) %>% 
                   dplyr::select(-FID,-PAT,-MAT,-SEX,-PHENOTYPE) %>% 
                   column_to_rownames(var = "IID") %>% 
                   as.matrix()) %>% 
  reduce(.,cbind)
dim(snps)
# [1]   251 44930
saveRDS(snps,file = paste0(pathIn,"DosageMatrix_DCas20_5440_WA_REFimputedAndFiltered.rds"))
```
