---
title: "Very GBS-DArT sample matches for EMBRAPA (DCas19_4301)"
site: workflowr::wflow_site
date: "2019-October-24"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: inline
---

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

# Make a directory for output
```{bash}
mkdir /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419
```

# Sample Lists
```{r}
library(tidyverse); library(magrittr);
system(paste0("bcftools query --list-samples ",
              "/workdir/marnin/embrapa_dataset_102419/tassel_embrapa_newfastq_filter/",
              "tassel_embrapa_newfastq.chr13.bial.recode.vcf.gz ",
              "> /workdir/marnin/embrapa_dataset_102419/tassel_embrapa_newfastq_filter/",
              "tassel_embrapa_newfastq.chr13.bial.recode.samples"))
system(paste0("bcftools query --list-samples ",
              "/workdir/marnin/nextgenImputation2019/",
              "DCas19_4301_DArTseqLD_AllSites_AllChrom_raw_70919.vcf.gz ",
              "> /workdir/marnin/nextgenImputation2019/",
              "DCas19_4301_DArTseqLD_AllSites_AllChrom_raw_70919.samples"))

embrapa_dataset_samples<-read.table(paste0("/workdir/marnin/embrapa_dataset_102419/tassel_embrapa_newfastq_filter/",
                                           "tassel_embrapa_newfastq.chr13.bial.recode.samples"),header = F, stringsAsFactors = F)$V1
dcas19_4301_samples<-read.table(paste0("/workdir/marnin/nextgenImputation2019/",
                                       "DCas19_4301_DArTseqLD_AllSites_AllChrom_raw_70919.samples"),header = F, stringsAsFactors = F)$V1
gbs2dartNames<-read.table("/workdir/marnin/EMBRAPA/GenotypingNamesDArTandGBS_80318.txt",
                          header = T, stringsAsFactors = F,sep = '\t')

table(gbs2dartNames$NameGBS %in% embrapa_dataset_samples) # All true
table(gbs2dartNames$NameDArT %in% dcas19_4301_samples) # None
gbs2dartNames$NameDArT[1:10]
dcas19_4301_samples[1:10]
table(gbs2dartNames$NameDArT %in% dcas19_4301_samples) # None
table(gsub("-",".",gbs2dartNames$NameDArT) %in% dcas19_4301_samples) # 869 true
gsub("-",".",gbs2dartNames$NameDArT)[!gsub("-",".",gbs2dartNames$NameDArT) %in% dcas19_4301_samples]

gbs2dartNames %<>% 
    mutate(NameDArT=gsub("-",".",NameDArT)) %>% 
    filter(NameDArT %in% dcas19_4301_samples,
           NameGBS %in% embrapa_dataset_samples)
```

```{r}
path<-"/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/"
saveRDS(gbs2dartNames,file=paste0(path,"dart2gbs_NameMatchesToVerifyByIBD_102419.rds"))
write.table(unique(gbs2dartNames$NameDArT),
            file=paste0(path,
                        "dartSamples_toVerifyAgainstGBS_102419.txt"),
            row.names = F, col.names = F, quote = F)
write.table(unique(gbs2dartNames$NameGBS),
            file=paste0(path,
                        "gbsSamples_toVerifyAgainstDArT_102419.txt"),
            row.names = F, col.names = F, quote = F)
```

# Site Lists

## Get sitesWithAlleles lists

```{r}
library(tidyverse); library(magrittr); library(furrr); options(mc.cores=18); plan(multiprocess)
tibble(Chr=1:18) %>%
  mutate(SiteWithAlleles=future_map(Chr,function(Chr){ 
    system(paste0("zcat /workdir/marnin/embrapa_dataset_102419/tassel_embrapa_newfastq_filter/",
                  "tassel_embrapa_newfastq.chr",Chr,".bial.recode.vcf.gz ",
                  "| cut -f1-5 > ",
                  "/workdir/marnin/embrapa_dataset_102419/tassel_embrapa_newfastq_filter/",
                  "tassel_embrapa_newfastq.chr",Chr,".bial.recode.sitesWithAlleles"))}))

embrapa_dataset_sites<-tibble(Chr=1:18) %>% 
  mutate(SiteWithAlleles=future_map(Chr,function(Chr){ 
    sites<-read.table(file = paste0("/workdir/marnin/embrapa_dataset_102419/tassel_embrapa_newfastq_filter/",
                                    "tassel_embrapa_newfastq.chr",Chr,".bial.recode.sitesWithAlleles"),
                      stringsAsFactors = F, header = F) }))

dart_sites<-read.table(paste0("/workdir/marnin/nextgenImputation2019/",
                              "DCas19_4301_DArTseqLD_AllSites_AllChrom_raw_70919.sitesWithAlleles"),
                       stringsAsFactors = F, header = T)

```

## DArT-GBS intersection
```{r}
dart_gbs_intersection<-embrapa_dataset_sites %>% 
    unnest(SiteWithAlleles) %>% 
    rename(CHROM=V1,
         POS=V2,
         ID=V3,
         REF=V4,
         ALT=V5) %>% 
  semi_join(dart_sites)
dart_gbs_intersection %>% dim # 1795
dart_gbs_intersection %>% 
    count(Chr)
saveRDS(dart_gbs_intersection,file=paste0("/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
                                          "dcas19_4301_embrapa_dataset_intersection_102419.rds"))
```

# Extract from raw VCFs
```{r}
library(tidyverse); library(magrittr); require(furrr); options(mc.cores=18); plan(multiprocess)
dart_gbs_intersection %>% 
    select(Chr,CHROM,POS) %>% 
    nest(isect_sites=c(CHROM,POS)) %>% 
    mutate(write_isect_sites=future_map2(Chr,isect_sites,function(Chr,isect_sites){
    write.table(isect_sites,
                file=paste0("/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                            Chr,"_dcas19_4301_embrapa_dataset_intersection_102419.sites"),
                row.names = F, col.names = F, quote = F) }))

future_map(1:18,~system(paste0("vcftools --gzvcf /workdir/marnin/embrapa_dataset_102419/tassel_embrapa_newfastq_filter/",
                               "tassel_embrapa_newfastq.chr",.,".bial.recode.vcf.gz ",
                               "--keep /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
                               "gbsSamples_toVerifyAgainstDArT_102419.txt ",
                               "--positions /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                               .,"_dcas19_4301_embrapa_dataset_intersection_102419.sites ",
                               "--recode ",
                               "--stdout | bgzip -c -@ 24 > ",
                               "/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                               .,"_gbsSamples_EMBRAPA_gbsDartIntersectingSites_102419.vcf.gz")))

future_map(1:18,~system(paste0("vcftools --gzvcf /workdir/marnin/nextgenImputation2019/",
                               "DCas19_4301_DArTseqLD_AllSites_AllChrom_raw_70919.vcf.gz ",
                               "--keep /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
                               "dartSamples_toVerifyAgainstGBS_102419.txt ",
                               "--positions /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                               .,"_dcas19_4301_embrapa_dataset_intersection_102419.sites ",
                               "--recode ",
                               "--stdout | bgzip -c -@ 24 > ",
                               "/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                               .,"_dartSamples_EMBRAPA_gbsDartIntersectingSites_102419.vcf.gz")))
```
# Index
```{r}
tibble(Chr=1:18) %>%
  mutate(Index=future_map(Chr,function(Chr){ 
    system(paste0("tabix -f -p vcf /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                  Chr,"_gbsSamples_EMBRAPA_gbsDartIntersectingSites_102419.vcf.gz")) 
      system(paste0("tabix -f -p vcf /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                    Chr,"_dartSamples_EMBRAPA_gbsDartIntersectingSites_102419.vcf.gz"))}))

```
# Merge GBS-DArT
```{r}
tibble(Chr=1:18) %>%
  mutate(Merge=future_map(Chr,function(Chr){ 
    system(paste0("bcftools merge ",
                  "--output /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                  Chr,"_AllSamplesToVerify_gbsDartIntersectingSites_102419.vcf.gz ",
                  "--merge snps --output-type z --threads 6 ",
                  "/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                  Chr,"_gbsSamples_EMBRAPA_gbsDartIntersectingSites_102419.vcf.gz ",
                  "/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                  Chr,"_dartSamples_EMBRAPA_gbsDartIntersectingSites_102419.vcf.gz")) }))
```
# Pre-analysis filter
```{r}
tibble(Chr=1:18) %>%
  mutate(PreAnalysisFilter=future_map(Chr,function(Chr){ 
    system(paste0("vcftools --gzvcf /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                  Chr,"_AllSamplesToVerify_gbsDartIntersectingSites_102419.vcf.gz ",
                  "--min-alleles 2 --max-alleles 2 --minDP 4 --maxDP 50 ", 
                  "--recode --stdout | bgzip -c -@ 24 > ",
                  "/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                  Chr,"_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419.vcf.gz"))}))

```
# Concat chroms
```{r}
tibble(Chr=1:18) %>%
  mutate(Index=future_map(Chr,function(Chr){ 
    system(paste0("tabix -f -p vcf /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",
                  Chr,"_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419.vcf.gz"))}))

system(paste0("bcftools concat --allow-overlaps ",
              "--output /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
              "AllChrom_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419.vcf.gz ",
              "--output-type z --threads 6 ",
              paste0("/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/chr",1:18,
                     "_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419.vcf.gz",collapse = " ")))
```
# Make binary blink (bed/bim/fam)
```{r}
system(paste0("export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;",
                      "plink --vcf /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
                      "AllChrom_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419.vcf.gz ", 
                      "--make-bed --const-fid ",
                      "--out /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
                      "AllChrom_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419")) 
```
# PLINK IBD
```{r}
system(paste0("export PATH=/programs/plink-1.9-x86_64-beta3.30:$PATH;",
              "plink --bfile /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
              "AllChrom_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419 ",
              "--genome --out /workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
              "AllChrom_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419;"))
```
# Choose DArT-GBS matching records
```{r}
library(data.table)
genome<-fread(paste0("/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/",
                     "AllChrom_AllSamplesToVerify_gbsDartIntersectingSites_filtered_102419.genome"),
              stringsAsFactors = F, header = T) %>%
    as_tibble
embrapa_dataset_samples<-read.table(paste0("/workdir/marnin/embrapa_dataset_102419/tassel_embrapa_newfastq_filter/",
                                           "tassel_embrapa_newfastq.chr13.bial.recode.samples"),header = F, stringsAsFactors = F)$V1
dcas19_4301_samples<-read.table(paste0("/workdir/marnin/nextgenImputation2019/",
                                       "DCas19_4301_DArTseqLD_AllSites_AllChrom_raw_70919.samples"),header = F, stringsAsFactors = F)$V1
gbs2dartNames<-read.table("/workdir/marnin/EMBRAPA/GenotypingNamesDArTandGBS_80318.txt",
                          header = T, stringsAsFactors = F,sep = '\t') %>% 
                          mutate(NameDArT=gsub("-",".",NameDArT)) %>% 
    filter(NameDArT %in% dcas19_4301_samples,
           NameGBS %in% embrapa_dataset_samples)
```

# Make "official" matches
```{r}
IBDmatches<-gbs2dartNames %>% 
  left_join(
    genome %>% 
      filter(DST>=0.9 | PI_HAT>=0.65) %>% 
      select(IID1,IID2,DST,PI_HAT) %>% 
      rename(NameGBS=IID1,
             NameDArT=IID2)
    ) %>% 
  filter(!is.na(DST)) %>% # 723 pass
  arrange(NameGBS,desc(PI_HAT)) %>% 
  group_by(NameDArT) %>% 
  slice(1) %>% 
  group_by(NameGBS) %>% 
  slice(1) %>% 
  ungroup() # 723 (of 869) unique dart-gbs matches!
IBDmatches %>% count(NameDArT) %>% arrange(desc(n))
IBDmatches %>% count(NameGBS) %>% arrange(desc(n))
```
```{r}
saveRDS(IBDmatches,file="/workdir/marnin/nextgenImputation2019/VerifyByIBD_EMBRAPA_102419/gbs2dart_SamplesVerifiedByIBD_102419.rds")
```


