---
title: "Convert IITA-EA DCas20-5261 to VCF"
site: workflowr::wflow_site
date: "2020-November-11"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: inline
---

```{bash, eval=F}
cd /home/jj332_cas/marnin/IITA_EA_2020
```

# Input Parameters
```{r, eval=F}
#' @dartvcfInput input name and path of "vcf" file from DArT
#' @dartcountsInput input name and path of counts file from DArT
#' @outName output path and name
#' @nskipvcf number of "VCF" rows to skip on read-in
#' @nskipcounts number of "counts file" rows to skip on read in
#' @ncores number of cores to use, could be VERY memory intensive
#' @dartVars chr vector, column names that _are not_ sample IDs in the read-counts file. I use this internally to assign the sampleIDs in the VCF file

library(tidyverse); library(magrittr)
dartvcfInput<-here::here("data/OrderAppendix_1_DCas20-5261","Report_5261_VCF_Ref_Version6.txt")
dartcountsInput<-here::here("data/OrderAppendix_1_DCas20-5261","Report_5261_Counts_Ref_Version6.csv")
outName<-here::here("data/OrderAppendix_1_DCas20-5261","DCas20_5261")
nskipvcf<-2 
nskipcounts<-3
ncores<-10
```

# Prelim. check format

Start manual. Check that the files read in according to previous code. Adjust code if necessary. Make a function and apply it to the input files.

```{r, eval=F}
# vcf<-read.table(dartvcfInput,
#                 stringsAsFactors = F,skip = nskipvcf, header = T, sep = "\t", comment.char = "")
# readCounts<-read.csv(dartcountsInput, stringsAsFactors = F,header = T,skip=nskipcounts)
# 
# dim(vcf)
# # [1] 13603  1844
# dim(readCounts)
# [1] 27206  1878
# 
# # the lenghts don't match?
# 
# # Initial look at names....
# colnames(readCounts)[1:100]
 #  [1] "AlleleID"                    "CloneID"                    
 #  [3] "ClusterTempIndex"            "AlleleSequence"             
 #  [5] "TrimmedSequence"             "TrimmedSequence_plus_Strand"
 #  [7] "Short"                       "Lowcomplexity"              
 #  [9] "Chrom_Cassava_v61"           "ChromPos_Cassava_v61"       
 # [11] "SNP_ChromPos_Cassava_v61"    "AlnCnt_Cassava_v61"         
 # [13] "AlnEvalue_Cassava_v61"       "Strand_Cassava_v61"         
 # [15] "SeqDiff_Cassava_v61"         "ClusterConsensusSequence"   
 # [17] "ClusterSize"                 "AlleleSeqDist"              
 # [19] "SNP"                         "SnpPosition"                
 # [21] "CallRate"                    "OneRatioRef"                
 # [23] "OneRatioSnp"                 "FreqHomRef"                 
 # [25] "FreqHomSnp"                  "FreqHets"                   
 # [27] "PICRef"                      "PICSnp"                     
 # [29] "AvgPIC"                      "AvgCountRef"                
 # [31] "AvgCountSnp"                 "RatioAvgCountRefAvgCountSnp"
 # [33] "FreqHetsMinusFreqMinHom"     "AlleleCountsCorrelation"    
 # [35] "aggregateTagsTotal"          "DerivedCorrMinusSeedCorr"   
 # [37] "RepRef"                      "RepSNP"                     
 # [39] "RepAvg"                      "PicRepRef"                  
 # [41] "PicRepSNP"                   "TotalPicRepRefTest"         
 # [43] "TotalPicRepSnpTest"          "MM161583"                   
 # [45] "MM161292"                    "MM160311"                   
 # [47] "MM161370"                    "MM160873"                   
 # [49] "MM160630"                    "MM161194"                   
 # [51] "MM160219"                    "MM160923"                   
 # [53] "MM161595"                    "MM160730"                   
 # [55] "MM160120"                    "MM160695"
#colnames(vcf)[1:30]
#  [1] "X.CHROM"  "POS"      "ID"       "REF"      "ALT"      "QUAL"    
#  [7] "FILTER"   "INFO"     "FORMAT"   "MM161104" "MM160414" "MM161519"
# [13] "MM160032" "MM161583" "MM160923" "MM160093" "MM160166" "MM160243"
# [19] "MM160638" "MM160468" "MM160092" "MM160642" "MM160867" "MM160856"
# [25] "MM160789" "MM161292" "MM161595" "MM161011" "MM161354" "MM161122"

# rm(vcf,readCounts); gc()
```

# Conversion function

Available and sourced from `code/` subdirectory: `convertDart2vcf.R`.

# Run conversion function
```{r, eval=F}
source(here::here("code/","convertDart2vcf.R"))
convertDart2vcf(dartvcfInput,dartcountsInput,outName,
                nskipvcf=2,nskipcounts=3,ncores)
```

# Genomewide to per-chrom VCFs

Split the genome-wide VCF into per-chromosome VCFs for imputation.

```{r,eval = FALSE}
require(furrr); options(mc.cores=18); plan(multiprocess)
source(here::here("code","imputationFunctions.R"))

vcfIn<-here::here("data/OrderAppendix_1_DCas20-5261","DCas20_5261.vcf.gz")
filters<-"--minDP 4 --maxDP 50" # because using GT not PL for impute (Beagle5)
outPath<-here::here("data/OrderAppendix_1_DCas20-5261/")
outSuffix<-"DCas20_5261"

future_map(1:18,
           ~splitVCFbyChr(Chr=.,
                          vcfIn=vcfIn,filters=filters,
                          outPath=outPath,outSuffix=outSuffix))
```

