ImputeDCas20_5360.Rmd

---
title: "Impute EMBRAPA GS C2 (DCas20_5360)"
site: workflowr::wflow_site
date: "2020-August-26"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: inline
---

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

DArTseqLD (DCas20-5360) arrived on Aug. 22, 2020. Contains **GS C2** for EMBRAPA. 

Last year's reference panel for imputation had ~64K SNP. The C1 progeny imputed by it had <9K SNP _after_ post-imputation filters.

Options to proceed with imputation of C2 (DCas20-5360):

1. Last year's reference panel, without C1 (64K SNP)
2. Last year's reference panel + C1 (reduce refpanel to sites passing post-impute filter in C1 = 9K SNP)
3. Last year's ref. panel + C1 (include unfiltered C1 data = 64K refpanel SNP)

If I can use Beagle5 or latest, should be possible (fast) to create all 3 and compare the results using PCA, prediction, correlation of kinship matrices, etc.

# Impute EMBRAPA GS C2

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

```{bash,eval = FALSE}
cp /home/jj332_cas/marnin/EMBRAPA_2020GS /workdir/mw489/
cp /home/jj332_cas/CassavaGenotypeData/nextgenImputation2019/ImputationEMBRAPA_102419/chr*_ImputationReferencePanel_EMBRAPA_Phased_102619.vcf.gz /workdir/mw489/EMBRAPA_2020GS/data/
cp -r /home/jj332_cas/marnin/EMBRAPA_2020GS/code /workdir/mw489/EMBRAPA_2020GS/
cp -r /home/jj332_cas/CassavaGenotypeData/CassavaGeneticMap /workdir/mw489/EMBRAPA_2020GS/data/
cp /home/jj332_cas/CassavaGenotypeData/nextgenImputation2019/ImputationEMBRAPA_102419/chr*_EMBRAPA_C1progeny_FromDCas19_4301_*  /workdir/mw489/EMBRAPA_2020GS/data/ 
```

## RefPanel: Exclude C1

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_EMBRAPA_Phased_102619.vcf.gz` as reference.

Requires 1 large memory Cornell CBSU machine (e.g. cbsulm17; 112 cores, 512 GB RAM), 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.

```{r,eval = FALSE}
targetVCFpath<-here::here("data/Report-DCas20-5360/") # location of the targetVCF
refVCFpath<-here::here("data/")
mapPath<-here::here("data/CassavaGeneticMap/")
outPath<-here::here("output/")
samplesToExclude<-here::here("data","samplesToExclude.txt")
outSuffix<-"DCas20_5360"
```

Found a sample "A.Preto.joselito" in target VCF that is already in the refpanel, so I must exclude it. Created a file `data/samplesToExclude.txt` to give to Beagle.

```{bash, eval=F}
bcftools query \
  --list-samples /home/jj332_cas/CassavaGenotypeData/nextgenImputation2019/ImputationEMBRAPA_102419/chr1_ImputationReferencePanel_EMBRAPA_Phased_102619.vcf.gz > /home/jj332_cas/marnin/EMBRAPA_2020GS/data/chr1_ImputationReferencePanel_EMBRAPA_Phased_102619.samples
```

```{r, eval=F}
refpanel_samples<-read.table(here::here("data","chr1_ImputationReferencePanel_EMBRAPA_Phased_102619.samples"), 
                             header=F, stringsAsFactors = F)$V1
target_samples<-read.table(here::here("data/Report-DCas20-5360/","DCas20_5360.samples"), 
                             header=F, stringsAsFactors = F)$V1
samplesToExclude<-target_samples[target_samples %in% refpanel_samples]
length(samplesToExclude) # [1] 644
write.table(samplesToExclude,file=here::here("data","samplesToExclude.txt"),row.names = F, col.names = F, quote = F)
```


### Impute C2
```{r,eval = FALSE}
source(here::here("code","imputationFunctions.R"))

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

Clean up Beagle log files after run. Move to sub-directory `output/BeagleLogs/`.
```{bash,eval = FALSE}
cd /workdir/mw489/EMBRAPA_2020GS/output/; 
mkdir BeagleLogs;
cp *_DCas20_5360_REFimputed.log BeagleLogs/
cp -r BeagleLogs /home/jj332_cas/marnin/EMBRAPA_2020GS/output/
cp *_DCas20_5360_REFimputed* /home/jj332_cas/marnin/EMBRAPA_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_5360_REFimputed"),
                                  outPath=outPath,
                                  outName=paste0("chr",.,"_DCas20_5360_REFimputedAndFiltered")))
# [1] "1068 sites passing filter"
# [1] "Filtering Complete: chr1_DCas20_5360_REFimputedAndFiltered"
# [1] "437 sites passing filter"
# [1] "Filtering Complete: chr2_DCas20_5360_REFimputedAndFiltered"
# [1] "559 sites passing filter"
# [1] "Filtering Complete: chr3_DCas20_5360_REFimputedAndFiltered"
# [1] "207 sites passing filter"
# [1] "Filtering Complete: chr4_DCas20_5360_REFimputedAndFiltered"
# [1] "783 sites passing filter"
# [1] "Filtering Complete: chr5_DCas20_5360_REFimputedAndFiltered"
# [1] "447 sites passing filter"
# [1] "Filtering Complete: chr6_DCas20_5360_REFimputedAndFiltered"
# [1] "238 sites passing filter"
# [1] "Filtering Complete: chr7_DCas20_5360_REFimputedAndFiltered"
# [1] "409 sites passing filter"
# [1] "Filtering Complete: chr8_DCas20_5360_REFimputedAndFiltered"
# [1] "474 sites passing filter"
# [1] "Filtering Complete: chr9_DCas20_5360_REFimputedAndFiltered"
# [1] "243 sites passing filter"
# [1] "Filtering Complete: chr10_DCas20_5360_REFimputedAndFiltered"
# [1] "272 sites passing filter"
# [1] "Filtering Complete: chr11_DCas20_5360_REFimputedAndFiltered"
# [1] "361 sites passing filter"
# [1] "Filtering Complete: chr12_DCas20_5360_REFimputedAndFiltered"
# [1] "393 sites passing filter"
# [1] "Filtering Complete: chr13_DCas20_5360_REFimputedAndFiltered"
# [1] "755 sites passing filter"
# [1] "Filtering Complete: chr14_DCas20_5360_REFimputedAndFiltered"
# [1] "712 sites passing filter"
# [1] "Filtering Complete: chr15_DCas20_5360_REFimputedAndFiltered"
# [1] "280 sites passing filter"
# [1] "Filtering Complete: chr16_DCas20_5360_REFimputedAndFiltered"
# [1] "674 sites passing filter"
# [1] "Filtering Complete: chr17_DCas20_5360_REFimputedAndFiltered"
# [1] "610 sites passing filter"
# [1] "Filtering Complete: chr18_DCas20_5360_REFimputedAndFiltered"
```

Check what's left
```{r,eval = FALSE}
purrr::map(1:18,~system(paste0("zcat ",here::here("output/"),"chr",.,"_DCas20_5360_REFimputedAndFiltered.vcf.gz | wc -l")))
# 1078
# 447
# 569
# 217
# 793
# 457
# 248
# 419
# 484
# 253
# 282
# 371
# 403
# 765
# 722
# 290
# 684
# 620
```

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

## RefPanel: Include C1, all sites

### Merge RefPanel and C1
```{r, eval=F}
inPath<-here::here("data/")
outPath<-here::here("data/")
source(here::here("code","imputationFunctions.R"))
require(furrr); options(mc.cores=ncores); plan(multiprocess)
future_map(1:18,~mergeVCFs(inPath=inPath,
                           inVCF1=paste0("chr",.,"_ImputationReferencePanel_EMBRAPA_Phased_102619"),
                           inVCF2=paste0("chr",.,"_EMBRAPA_C1progeny_FromDCas19_4301_REFimputed_102619"),
                           outPath=outPath,
                           outName=paste0("chr",.,"_ImputationReferencePanel_C1progenyAdded_EMBRAPA")))
```

```{bash, eval=F}
cd data/ 
bcftools query --list-samples chr1_ImputationReferencePanel_C1progenyAdded_EMBRAPA.vcf.gz > chr1_ImputationReferencePanel_C1progenyAdded_EMBRAPA.samples
```

```{r, eval=F}
refpanel_samples<-read.table(here::here("data","chr1_ImputationReferencePanel_C1progenyAdded_EMBRAPA.samples"), 
                             header=F, stringsAsFactors = F)$V1
target_samples<-read.table(here::here("data/Report-DCas20-5360/","DCas20_5360.samples"), 
                             header=F, stringsAsFactors = F)$V1
samplesToExclude<-target_samples[target_samples %in% refpanel_samples]
length(samplesToExclude) # [1] 713
write.table(samplesToExclude,file=here::here("data","samplesToExclude_RefPanelC1Added.txt"),
            row.names = F, col.names = F, quote = F)
```

### Impute C2
```{r,eval = FALSE}
source(here::here("code","imputationFunctions.R"))
targetVCFpath<-here::here("data/Report-DCas20-5360/") # location of the targetVCF
refVCFpath<-here::here("data/")
mapPath<-here::here("data/CassavaGeneticMap/")
outPath<-here::here("output/")
samplesToExclude<-here::here("data","samplesToExclude_RefPanelC1Added.txt")

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

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

cd /workdir/mw489/EMBRAPA_2020GS/data/;
cp *_ImputationReferencePanel_C1progenyAdded_EMBRAPA* /home/jj332_cas/marnin/EMBRAPA_2020GS/data/;
cp samplesToExclude_RefPanelC1Added.txt /home/jj332_cas/marnin/EMBRAPA_2020GS/data/;
```

### Post-impute filter
```{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_5360_REFimputedWithC1unfiltered"),
                                  outPath=outPath,
                                  outName=paste0("chr",.,"_DCas20_5360_REFimputedWithC1unfiltered_PostImputeFiltered")))
```

```{r,eval = FALSE}
purrr::map(1:18,~system(paste0("zcat ",here::here("output/"),"chr",.,"_DCas20_5360_REFimputedWithC1unfiltered_PostImputeFiltered.vcf.gz | wc -l")))
# 1568
# 768
# 917
# 10 <--- ZERO MARKERS ON CHR. 4 (10 ROWS = THE HEADER)
# 1201
# 748
# 347
# 716
# 916
# 450
# 518
# 545
# 577
# 1214
# 1179
# 406
# 1008
# 831

```

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


## RefPanel: Include C1, filtered sites only

Keep only the sites in the RefPanel that passed post-imputation filter for the C1

### Filter sites from RefPanel+C1 dataset
```{r, eval=F}
inPath<-here::here("data/")
outPath<-here::here("data/")
source(here::here("code","imputationFunctions.R"))
require(furrr); options(mc.cores=18); plan(multiprocess)
furrr::future_map(1:18,~filter_positions(inPath=inPath,
                                         inVCF=paste0("chr",.,"_ImputationReferencePanel_C1progenyAdded_EMBRAPA.vcf.gz"),
                                         positionFile=paste0("chr",.,"_EMBRAPA_C1progeny_FromDCas19_4301_REFimputed_102619.sitesPassing"),
                                         outPath=outPath,
                                         outName=paste0("chr",.,"_ImputationReferencePanel_C1progenyAddedFilteredSites_EMBRAPA")))

```

### Impute C2
```{r,eval = FALSE}
source(here::here("code","imputationFunctions.R"))
targetVCFpath<-here::here("data/Report-DCas20-5360/") # location of the targetVCF
refVCFpath<-here::here("data/")
mapPath<-here::here("data/CassavaGeneticMap/")
outPath<-here::here("output/")
samplesToExclude<-here::here("data","samplesToExclude_RefPanelC1Added.txt")

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

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

cd /workdir/mw489/EMBRAPA_2020GS/data/;
cp *_ImputationReferencePanel_C1progenyAddedFilteredSites_EMBRAPA* /home/jj332_cas/marnin/EMBRAPA_2020GS/data/;

```

### Post-impute filter
```{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_5360_REFimputedWithC1filtered"),
                                  outPath=outPath,
                                  outName=paste0("chr",.,"_DCas20_5360_REFimputedWithC1filtered_PostImputeFiltered")))
```
Check how many sites left
```{r,eval = FALSE}
purrr::map(1:18,~system(paste0("zcat ",here::here("output/"),"chr",.,"_DCas20_5360_REFimputedWithC1filtered_PostImputeFiltered.vcf.gz | wc -l")))
# 945
# 367
# 376
# 212
# 696
# 385
# 177
# 451
# 395
# 198
# 266
# 369
# 289
# 611
# 702
# 253
# 359
# 514

```

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

# Lighter post-impute filters
```{bash, eval=F}
cd /home/jj332_cas/marnin/EMBRAPA_2020GS;
```

```{r,eval = FALSE}
inPath<-here::here("output/")
outPath<-here::here("output/")
source(here::here("code","imputationFunctions.R"))
require(furrr); options(mc.cores=18); plan(multiprocess)
# Re-filter the REFimputed dataset
future_map(1:18,~postImputeFilter(inPath=inPath,
                                  inName=paste0("chr",.,"_DCas20_5360_REFimputed"),
                                  outPath=outPath,
                                  outName=paste0("chr",.,"_DCas20_5360_REFimputeLightFiltDR2pt3"),
                                  DR2thresh=0.3))
# Re-filter the REFimputedWithC1unfiltered dataset
future_map(1:18,~postImputeFilter(inPath=inPath,
                                  inName=paste0("chr",.,"_DCas20_5360_REFimputedWithC1unfiltered"),
                                  outPath=outPath,
                                  outName=paste0("chr",.,"_DCas20_5360_REFimputedWithC1unfiltered_LightPostImputeFiltDR2pt3"),
                                  DR2thresh=0.3))
# Re-filter the REFimputedWithC1filtered dataset
future_map(1:18,~postImputeFilter(inPath=inPath,
                                  inName=paste0("chr",.,"_DCas20_5360_REFimputedWithC1filtered"),
                                  outPath=outPath,
                                  outName=paste0("chr",.,"_DCas20_5360_REFimputedWithC1filtered_LightPostImputeFiltDR2pt3"),
                                  DR2thresh=0.3))
# Check what's left 
## REFimputed dataset 
purrr::map(1:18,~system(paste0("zcat ",here::here("output/"),"chr",.,"_DCas20_5360_REFimputeLightFiltDR2pt3.vcf.gz | wc -l")))
# 2399
# 1443
# 1714
# 554   <---- THAT'S BETTER THAN BELOW
# 1966
# 1319
# 618
# 1176
# 1417
# 965
# 866
# 1082
# 1019
# 1850
# 1867
# 746
# 1812
# 1347

## REFimputedWithC1unfiltered dataset 
purrr::map(1:18,~system(paste0("zcat ",here::here("output/"),"chr",.,"_DCas20_5360_REFimputedWithC1unfiltered_LightPostImputeFiltDR2pt3.vcf.gz | wc -l")))
# 2735
# 1722
# 1845
# 10   <---- THAT'S NOT GOOD. ZERO CHR. 4 MARKERS
# 2154
# 1574
# 836
# 1349
# 1563
# 1327
# 1056
# 1243
# 1149
# 2006
# 1993
# 879
# 1948
# 1530

## REFimputedWithC1filtered dataset 
purrr::map(1:18,~system(paste0("zcat ",here::here("output/"),"chr",.,"_DCas20_5360_REFimputedWithC1filtered_LightPostImputeFiltDR2pt3.vcf.gz | wc -l")))
# 970
# 416
# 406
# 228
# 755
# 396
# 179
# 487
# 405
# 222
# 274
# 394
# 316
# 697
# 733
# 259
# 382
# 552
```

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

# Formats for downstream analysis

The function below will (1) convert the input VCF to plink1.9 binary format and (2) convert the plink binary to a dosage (0,1,2) matrix with special attention to which allele gets counted in the file.

**NOTICE:** I was worried about `plink1.9` changing allele codes between files. There is some risk the counted allele could switch between e.g. the reference panel and the progeny files because of allele freq. (see plink documentation). To avoid this, went to extra trouble: write a file suffixed `*.alleleToCount` listing SNP ID (column 1) and the ALT allele from the VCF (column 2). Pass the file to `plink1.9` using the `--recode-allele` flag to ensure all output dosages count the ALT allele consistent with the VCFs. The reason to use `plink1.9` is that `Beagle5` imputed files don't have a **DS** (dosage) field that can be directly extracted. Instead, phased genotypes e.g. `0|1` need to be converted to dosages (e.g. `0|1 --> 1`, `1|1 --> 2`). An alternative might be to extract the haplotypes using `vcftools` and manually (in R) computed the dosages; that would give most control but is slow.

```{bash, eval=F}
cd /home/jj332_cas/marnin/EMBRAPA_2020GS/;
```
```{r, eval=F}
library(tidyverse); library(magrittr);
source(here::here("code","imputationFunctions.R"))
require(furrr); options(mc.cores=18); plan(multiprocess)
pathOut<-here::here("output/")

# 9 different filesets to convert
## Imputation reference panels
future_map(1:18,~convertVCFtoDosage(pathIn="/home/jj332_cas/CassavaGenotypeData/nextgenImputation2019/ImputationEMBRAPA_102419/",
                                    pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_ImputationReferencePanel_EMBRAPA_Phased_102619")))
future_map(1:18,~convertVCFtoDosage(pathIn=here::here("data/"),pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_ImputationReferencePanel_C1progenyAdded_EMBRAPA")))
future_map(1:18,~convertVCFtoDosage(pathIn=here::here("data/"),pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_ImputationReferencePanel_C1progenyAddedFilteredSites_EMBRAPA")))
## Standard filtered imputed C2
future_map(1:18,~convertVCFtoDosage(pathIn=here::here("output/"),pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_DCas20_5360_REFimputedAndFiltered")))
future_map(1:18,~convertVCFtoDosage(pathIn=here::here("output/"),pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_DCas20_5360_REFimputedWithC1unfiltered_PostImputeFiltered")))
future_map(1:18,~convertVCFtoDosage(pathIn=here::here("output/"),pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_DCas20_5360_REFimputedWithC1filtered_PostImputeFiltered")))

## "Lightly" filtered imputed C2
future_map(1:18,~convertVCFtoDosage(pathIn=here::here("output/"),pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_DCas20_5360_REFimputedWithC1filtered_LightPostImputeFiltDR2pt3")))
future_map(1:18,~convertVCFtoDosage(pathIn=here::here("output/"),pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_DCas20_5360_REFimputedWithC1unfiltered_LightPostImputeFiltDR2pt3")))
future_map(1:18,~convertVCFtoDosage(pathIn=here::here("output/"),pathOut=pathOut,
                                    vcfName = paste0("chr",.,"_DCas20_5360_REFimputeLightFiltDR2pt3")))

# Genome-wide dosage (for use in R) for each dataset
## Imputation reference panels
createGenomewideDosage(pathIn = here::here("output/"), chroms=1:18, "_ImputationReferencePanel_EMBRAPA_Phased_102619")
createGenomewideDosage(pathIn = here::here("output/"), chroms=1:18, "_ImputationReferencePanel_C1progenyAdded_EMBRAPA")
createGenomewideDosage(pathIn = here::here("output/"), chroms=1:18, "_ImputationReferencePanel_C1progenyAddedFilteredSites_EMBRAPA")
## Standard filtered imputed C2
createGenomewideDosage(pathIn = here::here("output/"), chroms=1:18, "_DCas20_5360_REFimputedAndFiltered")
createGenomewideDosage(pathIn = here::here("output/"), chroms=c(1:3,5:18), # chrom 4 missing from this one 
                       "_DCas20_5360_REFimputedWithC1unfiltered_PostImputeFiltered")
createGenomewideDosage(pathIn = here::here("output/"), chroms=1:18, "_DCas20_5360_REFimputedWithC1filtered_PostImputeFiltered")
## "Lightly" filtered imputed C2
createGenomewideDosage(pathIn = here::here("output/"), chroms=1:18, "_DCas20_5360_REFimputedWithC1filtered_LightPostImputeFiltDR2pt3")
createGenomewideDosage(pathIn = here::here("output/"), chroms=c(1:3,5:18), # chrom 4 missing from this one 
                       "_DCas20_5360_REFimputedWithC1unfiltered_LightPostImputeFiltDR2pt3")
createGenomewideDosage(pathIn = here::here("output/"), chroms=1:18, "_DCas20_5360_REFimputeLightFiltDR2pt3")

```