KIAT_SNP.Rmd

---
title: "KIAT_SNP"
author: "Ruijuan Li"
date: "11/2/2016"
output: 
  html_document: 
    keep_md: yes
---
# vcfR is another tool to work with VCF file. 

# SNP calling freebayes 
```{r}
# 1) uniquely mapped reads, remove PCR duplicate, add read group & freebayes 
# /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/analysis/freebayes/freebayes_SNP_calling.sh 
# 2) call SNPs (found that to filter based on RPL & RPR, I need to call SNPs from Da-Ae & Da-Ol-1 seperately...) 
# /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/analysis/freebayes/run_freebayes.sh   
# 3) process raw vcf file
# /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/analysis/freebayes/process_vcf.sh 
# 4) extract only SNPs 

# 5) get only biallelic SNPs 
# extract_biallelic.sh input.vcf output.vcf 

# data location: /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/result/freebayes/Ae_Ol_modified_SNPonly_biallelic.recode.ann.vcf 
``` 

# work on Da-Ae & Da-Ol-1 seperately... 
```{r}
source("/Users/ruijuanli/Desktop/Brassica_project/KIAT_RNA_seq/analysis/function_BnRNAseq.R")

# import 
vcf.freebayes.Ae <- read.table("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/as_diploid/no_bolting_unique_mapped/Ae_SNP_biallelic.recode.vcf",as.is=T,na.strings = ".:.:.:.:.:.:.") 
dim(vcf.freebayes.Ae) # 1019803      10
# save(vcf.freebayes.Ae, file ="~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.Ae.Rdata")

vcf.freebayes.Ol <- read.table("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/as_diploid/no_bolting_unique_mapped/Ol_SNP_biallelic.recode.vcf",as.is=T,na.strings = ".:.:.:.:.:.:.") 
dim(vcf.freebayes.Ol) # 997209     10
save(vcf.freebayes.Ol, file  ="~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.Ol.Rdata")

# reformat 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.Ae.Rdata")
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.Ol.Rdata")

# Da-Ae 
# get header 
vcf.header.freebayes.Ae <- system("grep '#C' ~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/as_diploid/no_bolting_unique_mapped/Ae_SNP_biallelic.recode.vcf",intern = TRUE) 
vcf.header.freebayes.Ae <- sub("#","",vcf.header.freebayes.Ae) #get rid of the pound sign
vcf.header.freebayes.Ae <- unlist(strsplit(vcf.header.freebayes.Ae,split="\t"))

# use function to reformat vcf file 
vcf.freebayes.reform.Ae <- SNP.freebayes.reformat.Ae(vcf.freebayes.Ae, vcf.header.freebayes.Ae)
dim(vcf.freebayes.reform.Ae) # 1019803      17   
table(vcf.freebayes.reform.Ae$Ae_gt)
#    0/0    0/1    1/1     NA 
# 423018 531949 522721 201104 

save(vcf.freebayes.reform.Ae, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.reform.Ae.Rdata")

load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.reform.Ae.Rdata")

# this way no 0/0.... but still this can be used to understand RPL & RLR... 
# filter based on QUAL 
vcf.freebayes.HQ.Ae <- vcf.freebayes.reform.Ae[vcf.freebayes.reform.Ae$QUAL>40,] # QUAL score of 20, which means less than 0.01 probability that it isn't polymorphic  
dim(vcf.freebayes.HQ.Ae) # 798423     17  
nrow(vcf.freebayes.HQ.Ae) / nrow(vcf.freebayes.reform.Ae) # 0.7829189 of SNPs were retained with QUAL > 40

min(vcf.freebayes.HQ.Ae[vcf.freebayes.HQ.Ae$Ae_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.Ae[vcf.freebayes.HQ.Ae$Ae_gt=="1/1",]$QUAL)

# count the number 
table(vcf.freebayes.HQ.Ae$Ae_gt)
#    0/1    1/1 
# 285815 512608 

mean(vcf.freebayes.HQ.Ae$Ae_tot.depth, na.rm=T) # 146

# distribution of depth
png("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/read.depth.png", width=4, height=3, units="in", res=300)
hist(log10(vcf.freebayes.HQ.Ae$Ae_tot.depth), xlab = "log10(read depth)", ylab="Number of SNPs", main="read depth")  
abline(v=c(1,3), col="red")
dev.off()

save(vcf.freebayes.HQ.Ae, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.HQ.Ae.Rdata")
```

```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.HQ.Ae.Rdata")
# use 1000 as the max, 10 as min 
vcf.freebayes.HQ.1.Ae <- vcf.freebayes.HQ.Ae[(!is.na(vcf.freebayes.HQ.Ae$Ae_tot.depth) & 
                              vcf.freebayes.HQ.Ae$Ae_tot.depth > 10 & 
                              vcf.freebayes.HQ.Ae$Ae_tot.depth < 1000),]
dim(vcf.freebayes.HQ.1.Ae) # 441950     17  

table(vcf.freebayes.HQ.1.Ae$Ae_gt)
#    0/1    1/1 
# 224945 217005 

# # get genotype quality score, Freebayes has its own way to calculate GQ, it uses population info as the prior, different from GATK. 
# head(vcf.freebayes.HQ.1.Ae)
# gen.lik <- matrix(
#   unlist(strsplit(vcf.freebayes.HQ.1.Ae$Ae_gen.lik,split = ",")),
#   nrow=nrow(vcf.freebayes.HQ.1.Ae),
#   byrow=TRUE
#   )
# 
# head(gen.lik)
# 
# colnames(gen.lik) <- c("ref.lik","het.lik","alt.lik")
# 
# head(gen.lik)
# gen.lik <- as.data.frame(gen.lik)
# colnames(gen.lik)
# 
# vcf.freebayes.HQ.1.Ae <- cbind(vcf.freebayes.HQ.1.Ae,gen.lik,stringsAsFactors=FALSE)
# head(vcf.freebayes.HQ.1.Ae)
# dim(vcf.freebayes.HQ.1.Ae) # 441950     20 
# 
# for (i in 1:nrow(vcf.freebayes.HQ.1.Ae)){
#   if (median(c(vcf.freebayes.HQ.1.Ae[i,]$ref.lik, vcf.freebayes.HQ.1.Ae[i,]$het.lik, vcf.freebayes.HQ.1.Ae[i,]$alt.lik)) * (-10) >= 99) {
#     vcf.freebayes.HQ.1.Ae[i,]$Ae_GQ = 99
#   } else { 
#     vcf.freebayes.HQ.1.Ae[i,]$Ae_GQ = median(c(vcf.freebayes.HQ.1.Ae[i,]$ref.lik, vcf.freebayes.HQ.1.Ae[i,]$het.lik, vcf.freebayes.HQ.1.Ae[i,]$alt.lik) * (-10))
#   }
# }
# 
# vcf.freebayes.HQ.1.Ae$alt.lik[1]

#####

tmp.list <- strsplit(vcf.freebayes.HQ.1.Ae$INFO,split = ";")
length(tmp.list) # 493019 
length(tmp.list[[1]]) 

for (i in 1:length(tmp.list)){
    tmp.list[[i]] <- tmp.list[[i]][c(1:41)] 
}

# 1) seperate INFO column  
head(vcf.freebayes.HQ.1.Ae)
info <- matrix(
  unlist(tmp.list),
  nrow=nrow(vcf.freebayes.HQ.1.Ae),  
  byrow=TRUE
  )

head(info)

colnames(info) <- c("AB","ABP","AC","AF","AN","AO","CIGAR","DP","DPB","DPRA","EPP","EPPR","GTI","LEN","MEANALT",
                    "MQM","MQMR","NS","NUMALT","ODDS","PAIRED","PAIREDR","PAO","PQA","PQR","PRO","QA","QR","RO",
                    "RPL","RPP","RPPR","RPR","RUN","SAF","SAP","SAR","SRF","SRP","SRR","TYPE")

head(info)
info <- as.data.frame(info)
colnames(info)

for (i in 1:ncol(info)){
  info[,i] <- gsub("([[:print:]]+)(=)([[:print:]])","\\3",info[,i])
}

head(info)

vcf.freebayes.HQ.1.Ae <- cbind(vcf.freebayes.HQ.1.Ae,info,stringsAsFactors=FALSE)
head(vcf.freebayes.HQ.1.Ae)
dim(vcf.freebayes.HQ.1.Ae) # 493019     67 
vcf.freebayes.HQ.1.reform.Ae <- vcf.freebayes.HQ.1.Ae

vcf.freebayes.HQ.1.reform.Ae$MAF_Ae <- round(pmin(vcf.freebayes.HQ.1.reform.Ae$Ae_alt.depth, vcf.freebayes.HQ.1.reform.Ae$Ae_ref.depth)/vcf.freebayes.HQ.1.reform.Ae$Ae_tot.depth, digits = 2) 

hist(vcf.freebayes.HQ.1.reform.Ae$MAF_Ae)

vcf.freebayes.HQ.1.reform.Ae$RPR <- as.numeric(vcf.freebayes.HQ.1.reform.Ae$RPR)
vcf.freebayes.HQ.1.reform.Ae$RPL <- as.numeric(vcf.freebayes.HQ.1.reform.Ae$RPL)

## tail bias 
hist(log10(vcf.freebayes.HQ.1.reform.Ae$RPR))
hist(log10(vcf.freebayes.HQ.1.reform.Ae$RPL))
min(vcf.freebayes.HQ.1.reform.Ae$RPR)
min(vcf.freebayes.HQ.1.reform.Ae$RPL)
max(vcf.freebayes.HQ.1.reform.Ae$RPR + vcf.freebayes.HQ.1.reform.Ae$RPL)

vcf.freebayes.HQ.1.reform.Ae$MBR <- round(pmin(vcf.freebayes.HQ.1.reform.Ae$RPR, vcf.freebayes.HQ.1.reform.Ae$RPL)/(vcf.freebayes.HQ.1.reform.Ae$RPR + vcf.freebayes.HQ.1.reform.Ae$RPL), digits = 2) 


hist(vcf.freebayes.HQ.1.reform.Ae$MBR, breaks = 50)

png("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/MBR_Ae.png", width=4, height=3, units="in", res=300)
hist(vcf.freebayes.HQ.1.reform.Ae$MBR, breaks = 50, xlab = "minor reads balance ratio (MBR)", ylab="number of SNPs", main="")
abline(v=c(0.01), col="red")
dev.off() 

# cor(vcf.freebayes.HQ.1.reform.Ae$MBR, vcf.freebayes.HQ.1.reform.Ae$MAF_Ae)
# vcf.freebayes.HQ.1.reform.Ae.het <- vcf.freebayes.HQ.1.reform.Ae[vcf.freebayes.HQ.1.reform.Ae$Ae_gt=="1/1",]
# 
# cor(vcf.freebayes.HQ.1.reform.Ae.het$MBR, vcf.freebayes.HQ.1.reform.Ae.het$MAF_Ae)

vcf.freebayes.HQ.2.tmp.Ae <- vcf.freebayes.HQ.1.reform.Ae

vcf.freebayes.HQ.2.test1.Ae <- 
vcf.freebayes.HQ.2.tmp.Ae[(vcf.freebayes.HQ.2.tmp.Ae$MBR>0.01),]

# vcf.freebayes.HQ.2.test2 <- 
# vcf.freebayes.HQ.2.tmp[(vcf.freebayes.HQ.2.tmp$RPR>=0 & 
#                            vcf.freebayes.HQ.2.tmp$RPL>0),]

min(as.numeric(vcf.freebayes.HQ.2.test1.Ae$RPL))
min(as.numeric(vcf.freebayes.HQ.2.test1.Ae$RPR))

table(vcf.freebayes.HQ.2.test1.Ae$Ae_gt) 
#     0/1    1/1 
# 127827 187049 

head(vcf.freebayes.HQ.1.reform.Ae) 
# save(vcf.freebayes.HQ.1.reform.Ae, file = "~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.HQ.1.reform.Ae.Rdata")

### het problem 
### compare QUAL, read depth, between 0/1 and 1/1

head(vcf.freebayes.HQ.1.Ae)

min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.reform.Ae$Ae_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ae$Ae_gt=="1/1",]$QUAL)

min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ae$Ae_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ae$Ae_gt=="1/1",]$QUAL)

min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ol$Ol_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ol$Ol_gt=="1/1",]$QUAL)

min(vcf.freebayes.HQ.2.test1.Ol[vcf.freebayes.HQ.2.test1.Ol$Ol_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.2.test1.Ol[vcf.freebayes.HQ.2.test1.Ol$Ol_gt=="1/1",]$QUAL) 

library(ggplot2) 
p.QUAL.Ol <- ggplot(data = vcf.freebayes.HQ.2.test1.Ol) 
p.QUAL.Ol <- p.QUAL.Ol + geom_histogram(aes(log10(QUAL), fill=Ol_gt)) 
p.QUAL.Ol <- p.QUAL.Ol + facet_wrap(~Ol_gt, ncol = 1) 
p.QUAL.Ol <- p.QUAL.Ol +labs(list(x="log10(QUAL)", y="Number of SNPs"))  
p.QUAL.Ol <- p.QUAL.Ol + theme(legend.position = "none") 
p.QUAL.Ol 

# plot QUAL by Ae_gt  
library(ggplot2) 
p.QUAL.Ae <- ggplot(data = vcf.freebayes.HQ.2.test1.Ae) 
p.QUAL.Ae <- p.QUAL.Ae + geom_histogram(aes(log10(QUAL), fill=Ae_gt)) 
p.QUAL.Ae <- p.QUAL.Ae + facet_wrap(~Ae_gt, ncol = 1) 
p.QUAL.Ae <- p.QUAL.Ae +labs(list(x="log10(QUAL)", y="Number of SNPs"))  
p.QUAL.Ae <- p.QUAL.Ae + theme(legend.position = "none")
p.QUAL.Ae 

ggsave(p.QUAL.Ae, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/p.QUAL.Ae.png", height = 4, width = 6)

# plot total depth by Ae_gt 
p.tot.depth.Ae <- ggplot(data = vcf.freebayes.HQ.2.test1.Ae) 
p.tot.depth.Ae <- p.tot.depth.Ae + geom_histogram(aes(Ae_tot.depth, fill=Ae_gt)) + scale_x_log10() 
p.tot.depth.Ae <- p.tot.depth.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.tot.depth.Ae

# plot MBR
p.MBR.Ae <- ggplot(data = vcf.freebayes.HQ.2.test1.Ae) 
p.MBR.Ae <- p.MBR.Ae + geom_histogram(aes(MBR, fill=Ae_gt))  
p.MBR.Ae <- p.MBR.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.MBR.Ae  
ggsave(p.MBR.Ae, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/p.MBR.Ae.png", height = 4, width = 6)

# plot GQ
p.GQ.Ae <- ggplot(data = vcf.freebayes.HQ.1.reform.Ae) 
p.GQ.Ae <- p.GQ.Ae + geom_histogram(aes(Ae_gt.qual, fill=Ae_gt))  
p.GQ.Ae <- p.GQ.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.GQ.Ae 

### check GQ 
head(vcf.freebayes.HQ.1.reform.Ae[vcf.freebayes.HQ.1.reform.Ae$Ae_gt.qual<10,])[,c("Ae_gt.qual","Ae_gen.lik")]
max(vcf.freebayes.HQ.1.reform.Ae$Ae_gt.qual)


# alt/ref allele ratio
vcf.freebayes.HQ.2.test1.Ae.het <- 
vcf.freebayes.HQ.2.test1.Ae[vcf.freebayes.HQ.2.test1.Ae$Ae_gt=="0/1",]

vcf.freebayes.HQ.2.test1.Ae.het$allele_ratio <- round(vcf.freebayes.HQ.2.test1.Ae.het$Ae_alt.depth/vcf.freebayes.HQ.2.test1.Ae.het$Ae_ref.depth, digits = 2)

hist(log10(vcf.freebayes.HQ.2.test1.Ae.het$allele_ratio), breaks= 150) 


nrow(vcf.freebayes.HQ.2.test1.Ae.het[vcf.freebayes.HQ.2.test1.Ae.het$allele_ratio>0,])
nrow(vcf.freebayes.HQ.2.test1.Ae.het) 

pl <- ggplot(data = vcf.freebayes.HQ.2.test1.Ae.het) 
pl <- pl + geom_histogram(aes(allele_ratio)) 
pl <- pl +labs(list(x="alt.depth/ref.depth", y="nunber of SNPs"))
pl

# didn't see correlation between MAF & MBR... 
cor.test(vcf.freebayes.HQ.2.test1.Ae.het)
cor.test(vcf.freebayes.HQ.2.test1.Ae.het$MAF_Ae, vcf.freebayes.HQ.2.test1.Ae.het$MBR)
cor.test(vcf.freebayes.HQ.2.test1.Ae$MAF_Ae, vcf.freebayes.HQ.2.test1.Ae$MBR)

cor.test(vcf.freebayes.HQ.2.test1.Ae.het$MAF_Ae, vcf.freebayes.HQ.2.test1.Ae.het$QUAL)
cor.test(vcf.freebayes.HQ.2.test1.Ae$MAF_Ae, vcf.freebayes.HQ.2.test1.Ae$QUAL) 

hist(vcf.freebayes.HQ.2.test1.Ae.het$Ae_alt.depth/vcf.freebayes.HQ.2.test1.Ae.het$Ae_ref.depth, breaks = 50) 

#### 

# Da-Ol-1 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.Ol.Rdata")
# get header 
vcf.header.freebayes.Ol <- system("grep '#C' ~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/as_diploid/no_bolting_unique_mapped/Ol_SNP_biallelic.recode.vcf",intern = TRUE) 
vcf.header.freebayes.Ol <- sub("#","",vcf.header.freebayes.Ol) #get rid of the pound sign
vcf.header.freebayes.Ol <- unlist(strsplit(vcf.header.freebayes.Ol,split="\t"))

# use function to reformat vcf file 
vcf.freebayes.reform.Ol <- SNP.freebayes.reformat.Ol(vcf.freebayes.Ol, vcf.header.freebayes.Ol)
dim(vcf.freebayes.reform.Ol) # 1019803      17   
table(vcf.freebayes.reform.Ol$Ol_gt)
#    0/0    0/1    1/1     NA 
# 423018 531949 522721 201104 

save(vcf.freebayes.reform.Ol, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.reform.Ol.Rdata")

# filter based on QUAL 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.reform.Ol.Rdata") 
vcf.freebayes.HQ.Ol <- vcf.freebayes.reform.Ol[vcf.freebayes.reform.Ol$QUAL>40,] # QUAL score of 20, which means less than 0.01 probability that it isn't polymorphic  
dim(vcf.freebayes.HQ.Ol) # 777458     17 
nrow(vcf.freebayes.HQ.Ol) / nrow(vcf.freebayes.reform.Ol) # 70% of SNPs were retained with QUAL > 40

# count the number 
table(vcf.freebayes.HQ.Ol$Ol_gt)

mean(vcf.freebayes.HQ.Ol$Ol_tot.depth, na.rm=T) # 146 

# distribution of depth 
hist(log10(vcf.freebayes.HQ.Ol$Ol_tot.depth))

# use 1000 as the max, 10 as min 
vcf.freebayes.HQ.1.Ol <- vcf.freebayes.HQ.Ol[(!is.na(vcf.freebayes.HQ.Ol$Ol_tot.depth) & 
                              vcf.freebayes.HQ.Ol$Ol_tot.depth > 10 & 
                              vcf.freebayes.HQ.Ol$Ol_tot.depth < 1000),]
dim(vcf.freebayes.HQ.1.Ol) # 441950     17  

table(vcf.freebayes.HQ.1.Ol$Ol_gt)

#### compare 0/0 & 0/1
# plot QUAL by Ol_gt 
p.QUAL.Ol <- ggplot(data = vcf.freebayes.HQ.1.Ol) 
p.QUAL.Ol <- p.QUAL.Ol + geom_histogram(aes(QUAL, fill=Ol_gt)) + scale_x_log10() 
p.QUAL.Ol <- p.QUAL.Ol + facet_wrap(~Ol_gt, ncol = 1)
p.QUAL.Ol

# plot total depth by Ol_gt 
p.tot.depth.Ol <- ggplot(data = vcf.freebayes.HQ.1.Ol) 
p.tot.depth.Ol <- p.tot.depth.Ol + geom_histogram(aes(Ol_tot.depth, fill=Ol_gt)) + scale_x_log10() 
p.tot.depth.Ol <- p.tot.depth.Ol + facet_wrap(~Ol_gt, ncol = 1)
p.tot.depth.Ol 

######

tmp.list <- strsplit(vcf.freebayes.HQ.1.Ol$INFO,split = ";")
length(tmp.list) # 493019 
length(tmp.list[[1]]) 

for (i in 1:length(tmp.list)){
    tmp.list[[i]] <- tmp.list[[i]][c(1:41)] 
}

# 1) seperate INFO column  
head(vcf.freebayes.HQ.1.Ol)
info <- matrix(
  unlist(tmp.list),
  nrow=nrow(vcf.freebayes.HQ.1.Ol),  
  byrow=TRUE
  )

head(info)

colnames(info) <- c("AB","ABP","AC","AF","AN","AO","CIGAR","DP","DPB","DPRA","EPP","EPPR","GTI","LEN","MEANALT",
                    "MQM","MQMR","NS","NUMALT","ODDS","PAIRED","PAIREDR","PAO","PQA","PQR","PRO","QA","QR","RO",
                    "RPL","RPP","RPPR","RPR","RUN","SAF","SAP","SAR","SRF","SRP","SRR","TYPE")

head(info)
info <- as.data.frame(info)
colnames(info) 

for (i in 1:ncol(info)){
  info[,i] <- gsub("([[:print:]]+)(=)([[:print:]])","\\3",info[,i])
}

head(info)

vcf.freebayes.HQ.1.Ol <- cbind(vcf.freebayes.HQ.1.Ol,info,stringsAsFactors=FALSE)
head(vcf.freebayes.HQ.1.Ol)
dim(vcf.freebayes.HQ.1.Ol) # 493019     67 
vcf.freebayes.HQ.1.reform.Ol <- vcf.freebayes.HQ.1.Ol

vcf.freebayes.HQ.1.reform.Ol$MAF_Ol <- round(pmin(vcf.freebayes.HQ.1.reform.Ol$Ol_alt.depth, vcf.freebayes.HQ.1.reform.Ol$Ol_ref.depth)/vcf.freebayes.HQ.1.reform.Ol$Ol_tot.depth, digits = 2) 

hist(vcf.freebayes.HQ.1.reform.Ol$MAF_Ol)

vcf.freebayes.HQ.1.reform.Ol$RPR <- as.numeric(vcf.freebayes.HQ.1.reform.Ol$RPR)
vcf.freebayes.HQ.1.reform.Ol$RPL <- as.numeric(vcf.freebayes.HQ.1.reform.Ol$RPL)

## tail bias 
hist(log10(vcf.freebayes.HQ.1.reform.Ol$RPR)) 
hist(log10(vcf.freebayes.HQ.1.reform.Ol$RPL))
min(vcf.freebayes.HQ.1.reform.Ol$RPR)
min(vcf.freebayes.HQ.1.reform.Ol$RPL)

vcf.freebayes.HQ.1.reform.Ol$MBR <- round(pmin(vcf.freebayes.HQ.1.reform.Ol$RPR, vcf.freebayes.HQ.1.reform.Ol$RPL)/(vcf.freebayes.HQ.1.reform.Ol$RPR + vcf.freebayes.HQ.1.reform.Ol$RPL), digits = 2) 

hist(vcf.freebayes.HQ.1.reform.Ol$MBR, xlim = c(0, 0.1), breaks = 50)
save(vcf.freebayes.HQ.1.reform.Ol, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.freebayes.HQ.1.reform.Ol.Rdata")

# cor(vcf.freebayes.HQ.1.reform.Ae$MBR, vcf.freebayes.HQ.1.reform.Ae$MAF_Ae)
# vcf.freebayes.HQ.1.reform.Ae.het <- vcf.freebayes.HQ.1.reform.Ae[vcf.freebayes.HQ.1.reform.Ae$Ae_gt=="1/1",]
# 
# cor(vcf.freebayes.HQ.1.reform.Ae.het$MBR, vcf.freebayes.HQ.1.reform.Ae.het$MAF_Ae)
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.freebayes.HQ.1.reform.Ol.Rdata")

vcf.freebayes.HQ.2.tmp.Ol <- vcf.freebayes.HQ.1.reform.Ol

vcf.freebayes.HQ.2.test1.Ol <- 
vcf.freebayes.HQ.2.tmp.Ol[(vcf.freebayes.HQ.2.tmp.Ol$MBR>0.01),] 

min(as.numeric(vcf.freebayes.HQ.2.test1.Ol$RPL))
min(as.numeric(vcf.freebayes.HQ.2.test1.Ol$RPR))

table(vcf.freebayes.reform.Ol$Ol_gt)
table(vcf.freebayes.HQ.2.test1.Ol$Ol_gt) 
 #   0/0    0/1    1/1 
 # 96721 112112 145987 
 
save(vcf.freebayes.HQ.1.reform.Ol, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.HQ.1.reform.Ol.Rdata") 
save(vcf.freebayes.HQ.2.test1.Ol, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.HQ.2.test1.Ol")
save(vcf.freebayes.HQ.2.test1.Ae, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.HQ.2.test1.Ae")

head(vcf.freebayes.HQ.1.reform.Ae[,c("CHROM", "POS", "RPL", "RPR")])

######## visualize in IGV ###################################### 
vcf.freebayes.HQ.2.test3.Ae <-
vcf.freebayes.HQ.2.tmp.Ae[(vcf.freebayes.HQ.2.tmp.Ae$RPL>0),] 

# find filtered SNPs because of tail bias, and visulize them in IGV 
test1 <- vcf.freebayes.HQ.2.tmp.Ae[,c("CHROM", "POS")]
test2 <- vcf.freebayes.HQ.2.test3.Ae[,c("CHROM", "POS")]
library(dplyr)
filtered_out_snp_tail_bias <- anti_join(test1, test2)
head(filtered_out_snp_tail_bias[filtered_out_snp_tail_bias$CHROM=="chrA01",], 50)
```  

# combine Da-Ae & Da-Ol data  
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.HQ.2.test1.Ae")
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.HQ.2.test1.Ol")

dim(vcf.freebayes.HQ.2.test1.Ae) 
dim(vcf.freebayes.HQ.2.test1.Ol)

table(vcf.freebayes.HQ.2.test1.Ae$Ae_gt)
table(vcf.freebayes.HQ.2.test1.Ol$Ol_gt)

# double check QUAL, depth, and MBR
min(vcf.freebayes.HQ.2.test1.Ae$QUAL)
min(vcf.freebayes.HQ.2.test1.Ae$Ae_tot.depth)
max(vcf.freebayes.HQ.2.test1.Ae$Ae_tot.depth)
min(vcf.freebayes.HQ.2.test1.Ae$MBR)

min(vcf.freebayes.HQ.2.test1.Ol$QUAL)
min(vcf.freebayes.HQ.2.test1.Ol$Ol_tot.depth)
max(vcf.freebayes.HQ.2.test1.Ol$Ol_tot.depth)
min(vcf.freebayes.HQ.2.test1.Ol$MBR)

# combine 
colnames(vcf.freebayes.HQ.2.test1.Ol)

Ae <- vcf.freebayes.HQ.2.test1.Ae[,c("CHROM", "POS", "REF", "ALT", "Ae_gt")]
Ol <- vcf.freebayes.HQ.2.test1.Ol[,c("CHROM", "POS", "REF", "ALT", "Ol_gt")]

head(Ae)
Ae$feature <- paste(Ae$CHROM, Ae$POS, Ae$REF, Ae$ALT, sep = "_")
Ol$feature <- paste(Ol$CHROM, Ol$POS, Ol$REF, Ol$ALT, sep = "_")

freebayes_result <- merge(Ae, Ol, by="feature", all=T)
?merge
dim(freebayes_result)
head(freebayes_result) 

freebayes_result_reform <- freebayes_result[,c("feature", "Ae_gt", "Ol_gt")]
head(freebayes_result_reform)  
freebayes_result_reform$Ae_gt <- gsub("NA", "1/1", freebayes_result_reform$Ae_gt)
freebayes_result_reform[is.na(freebayes_result_reform$Ae_gt)] <- "0/0"

freebayes_result_reform[which(is.na(freebayes_result_reform$Ae_gt)),]$Ae_gt <- "0/0"
freebayes_result_reform[which(is.na(freebayes_result_reform$Ol_gt)),]$Ol_gt <- "0/0"
head(freebayes_result_reform)  

freebayes_result_1 <- freebayes_result_reform[freebayes_result_reform$Ae_gt!=freebayes_result_reform$Ol_gt & freebayes_result_reform$Ae_gt!="0/1" & freebayes_result_reform$Ol_gt!="0/1",]  

dim(freebayes_result_1) # 215702
head(freebayes_result_1, 100)
table(freebayes_result_1$Ae_gt)
table(freebayes_result_1$Ol_gt)

save(freebayes_result_1, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/freebayes_result_1.Rdata")
###### save each step afterwards... 
```

# test correlation between MAF & MBR (whether skewed distribution of hetergyzougs is due to realignment...)
```{r}
cor.test(vcf.freebayes.HQ.1.reform.Ol$MBR, vcf.freebayes.HQ.1.reform.Ol$MAF_Ol)
cor.test(vcf.freebayes.HQ.1.reform.Ae$MBR, vcf.freebayes.HQ.1.reform.Ae$MAF_Ae)

vcf.freebayes.HQ.1.reform.Ae.het <-  vcf.freebayes.HQ.1.reform.Ae[vcf.freebayes.HQ.1.reform.Ae$Ae_gt=="0/1",]

vcf.freebayes.HQ.1.reform.Ae.homo <-  vcf.freebayes.HQ.1.reform.Ae[vcf.freebayes.HQ.1.reform.Ae$Ae_gt!="0/1",]

head(vcf.freebayes.HQ.1.reform.Ae.het[,c("CHROM", "POS")],100)

hist(vcf.freebayes.HQ.1.reform.Ae.het$MBR)
hist(vcf.freebayes.HQ.1.reform.Ae.homo$MBR, ylim = c(1, 1e+5))

### great, seems like most of the het calls are due to realignment problem
# but there are still some het calls, why? allele specific expression? still need to test that ... 
```

### analyze combined Da-Ae & Da-Ol-1 data  
```{r}
# source("/Users/ruijuanli/Desktop/Brassica_project/KIAT_RNA_seq/analysis/function_BnRNAseq.R")
# vcf.freebayes.no.ANN <- read.table("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/as_diploid/no_bolting_unique_mapped/Ae_Ol_modified_SNPonly_biallelic.recode.vcf",as.is=T,na.strings = ".:.:.:.:.:.:.") 
# dim(vcf.freebayes.no.ANN) # 1678792      11
# # save(vcf.freebayes, file ="~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.no.ANN.Rdata")
# load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.Rdata")
# 
# vcf.header.freebayes <- system("grep '#C' ~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/as_diploid/no_bolting_unique_mapped/Ae_Ol_modified_SNPonly_biallelic.recode.vcf",intern = TRUE) 
# vcf.header.freebayes <- sub("#","",vcf.header.freebayes) #get rid of the pound sign
# vcf.header.freebayes <- unlist(strsplit(vcf.header.freebayes,split="\t"))
# 
# # use function to reformat vcf file 
# vcf.freebayes.no.ANN.reform <- SNP.freebayes.reformat(vcf.freebayes.no.ANN, vcf.header.freebayes)
# dim(vcf.freebayes.no.ANN.reform) # 1678792      25  
# table(vcf.freebayes.no.ANN.reform$Ae_gt)
# #    0/0    0/1    1/1     NA 
# # 423018 531949 522721 201104 
# 
# table(vcf.freebayes.no.ANN.reform$Ol_gt)
# #    0/0    0/1    1/1     NA 
# # 459651 517580 518072 183489 
# 
# # 1) seperate INFO column  
# info <- matrix(
#   unlist(strsplit(vcf.freebayes.no.ANN.reform$INFO,split = ";")),
#   nrow=nrow(vcf.freebayes.no.ANN.reform),  
#   byrow=TRUE
#   )
# 
# head(info)
# 
# colnames(info) <- c("AB","ABP","AC","AF","AN","AO","CIGAR","DP","DPB","DPRA","EPP","EPPR","GTI","LEN","MEANALT",
#                     "MQM","MQMR","NS","NUMALT","ODDS","PAIRED","PAIREDR","PAO","PQA","PQR","PRO","QA","QR","RO",
#                     "RPL","RPP","RPPR","RPR","RUN","SAF","SAP","SAR","SRF","SRP","SRR","TYPE")
# 
# info <- as.data.frame(info)
# colnames(info)
# 
# for (i in 1:ncol(info)){
#   info[,i] <- gsub("([[:print:]]+)(=)([[:print:]])","\\3",info[,i])
# }
# 
# vcf.freebayes.no.ANN.reform <- cbind(vcf.freebayes.no.ANN,info,stringsAsFactors=FALSE)
# unique(vcf.freebayes.no.ANN.reform$CIGAR) 
# 
# save(vcf.freebayes.no.ANN.reform, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.no.ANN.reform.Rdata")

# It turns out SNPs all get CIGAR score of "1X". 
```

### Freebayes w/ reported 
```{r}
# data location: /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/result/freebayes/from_cabernet_call_separately 
``` 

# work on Da-Ae & Da-Ol-1 seperately... 
```{r}
source("~/KIAT/function_BnRNAseq.R")

# import 
vcf.freebayes.Ae <- read.table("~/SNP_parent/result/freebayes/from_cabernet_call_separately/Ae_SNP_bi.recode.vcf",as.is=T,na.strings = ".") 
dim(vcf.freebayes.Ae) # 1019617      10
# save(vcf.freebayes.Ae, file ="~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.Ae.Rdata")

vcf.freebayes.Ol <- read.table("~/SNP_parent/result/freebayes/from_cabernet_call_separately/Ol_SNP_bi.recode.vcf",as.is=T,na.strings = ".") 
dim(vcf.freebayes.Ol) # 997357     10
# save(vcf.freebayes.Ol, file  ="~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.Ol.Rdata")

# reformat 
# load("~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.Ae.Rdata")
# load("~/SNP_parent/result/freebayes/from_cabernet_call_separately/output//vcf.freebayes.Ol.Rdata")

# Da-Ae 
# get header 
vcf.header.freebayes.Ae <- system("grep '#C' ~/SNP_parent/result/freebayes/from_cabernet_call_separately/Ae_SNP_bi.recode.vcf",intern = TRUE) 
vcf.header.freebayes.Ae <- sub("#","",vcf.header.freebayes.Ae) #get rid of the pound sign
vcf.header.freebayes.Ae <- unlist(strsplit(vcf.header.freebayes.Ae,split="\t"))

# use function to reformat vcf file 
vcf.freebayes.reform.Ae <- SNP.freebayes.reformat.Ae.GQ(vcf.freebayes.Ae, vcf.header.freebayes.Ae)
dim(vcf.freebayes.reform.Ae) # 1019617      19 
table(vcf.freebayes.reform.Ae$Ae_gt)
#    0/0    0/1    1/1     NA 
# 423018 531949 522721 201104 

save(vcf.freebayes.reform.Ae, file = "~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.reform.Ae.Rdata")

load("~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.reform.Ae.Rdata")

# this way no 0/0.... but still this can be used to understand RPL & RLR... 
# filter based on QUAL 
vcf.freebayes.HQ.Ae <- vcf.freebayes.reform.Ae[vcf.freebayes.reform.Ae$QUAL>40,] # QUAL score of 20, which means less than 0.01 probability that it isn't polymorphic  
dim(vcf.freebayes.HQ.Ae) # 798730     19  
nrow(vcf.freebayes.HQ.Ae) / nrow(vcf.freebayes.reform.Ae) # 0.7829189 of SNPs were retained with QUAL > 40

min(vcf.freebayes.HQ.Ae[vcf.freebayes.HQ.Ae$Ae_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.Ae[vcf.freebayes.HQ.Ae$Ae_gt=="1/1",]$QUAL)

# count the number 
table(vcf.freebayes.HQ.Ae$Ae_gt)
#    0/1    1/1 
# 285815 512608 

mean(vcf.freebayes.HQ.Ae$Ae_tot.depth, na.rm=T) # 146

# distribution of depth
png("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/read.depth.png", width=4, height=3, units="in", res=300)
hist(log10(vcf.freebayes.HQ.Ae$Ae_tot.depth), xlab = "log10(read depth)", ylab="Number of SNPs", main="read depth")  
abline(v=c(1,3), col="red")
dev.off()

save(vcf.freebayes.HQ.Ae, file = "~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.HQ.Ae.Rdata")
```

```{r}
load("~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.HQ.Ae.Rdata")
# use 1000 as the max, 10 as min 
vcf.freebayes.HQ.1.Ae <- vcf.freebayes.HQ.Ae[(!is.na(vcf.freebayes.HQ.Ae$Ae_tot.depth) & 
                              vcf.freebayes.HQ.Ae$Ae_tot.depth > 10 & 
                              vcf.freebayes.HQ.Ae$Ae_tot.depth < 1000),]
dim(vcf.freebayes.HQ.1.Ae) # 442050     19

table(vcf.freebayes.HQ.1.Ae$Ae_gt)
#    0/1    1/1 
# 224945 217005 

# # get genotype quality score, Freebayes has its own way to calculate GQ, it uses population info as the prior, different from GATK. 
# head(vcf.freebayes.HQ.1.Ae)
# gen.lik <- matrix(
#   unlist(strsplit(vcf.freebayes.HQ.1.Ae$Ae_gen.lik,split = ",")),
#   nrow=nrow(vcf.freebayes.HQ.1.Ae),
#   byrow=TRUE
#   )
# 
# head(gen.lik)
# 
# colnames(gen.lik) <- c("ref.lik","het.lik","alt.lik")
# 
# head(gen.lik)
# gen.lik <- as.data.frame(gen.lik)
# colnames(gen.lik)
# 
# vcf.freebayes.HQ.1.Ae <- cbind(vcf.freebayes.HQ.1.Ae,gen.lik,stringsAsFactors=FALSE)
# head(vcf.freebayes.HQ.1.Ae)
# dim(vcf.freebayes.HQ.1.Ae) # 441950     20 
# 
# for (i in 1:nrow(vcf.freebayes.HQ.1.Ae)){
#   if (median(c(vcf.freebayes.HQ.1.Ae[i,]$ref.lik, vcf.freebayes.HQ.1.Ae[i,]$het.lik, vcf.freebayes.HQ.1.Ae[i,]$alt.lik)) * (-10) >= 99) {
#     vcf.freebayes.HQ.1.Ae[i,]$Ae_GQ = 99
#   } else { 
#     vcf.freebayes.HQ.1.Ae[i,]$Ae_GQ = median(c(vcf.freebayes.HQ.1.Ae[i,]$ref.lik, vcf.freebayes.HQ.1.Ae[i,]$het.lik, vcf.freebayes.HQ.1.Ae[i,]$alt.lik) * (-10))
#   }
# }
# 
# vcf.freebayes.HQ.1.Ae$alt.lik[1]

#####

tmp.list <- strsplit(vcf.freebayes.HQ.1.Ae$INFO,split = ";")
length(tmp.list) # 442050
length(tmp.list[[1]]) 

for (i in 1:length(tmp.list)){
    tmp.list[[i]] <- tmp.list[[i]][c(1:41)] 
}

# 1) seperate INFO column  
head(vcf.freebayes.HQ.1.Ae)
info <- matrix(
  unlist(tmp.list),
  nrow=nrow(vcf.freebayes.HQ.1.Ae),  
  byrow=TRUE
  )

head(info)

colnames(info) <- c("AB","ABP","AC","AF","AN","AO","CIGAR","DP","DPB","DPRA","EPP","EPPR","GTI","LEN","MEANALT",
                    "MQM","MQMR","NS","NUMALT","ODDS","PAIRED","PAIREDR","PAO","PQA","PQR","PRO","QA","QR","RO",
                    "RPL","RPP","RPPR","RPR","RUN","SAF","SAP","SAR","SRF","SRP","SRR","TYPE")

head(info)
info <- as.data.frame(info)
colnames(info)

for (i in 1:ncol(info)){
  info[,i] <- gsub("([[:print:]]+)(=)([[:print:]])","\\3",info[,i])
}

head(info)

vcf.freebayes.HQ.1.Ae <- cbind(vcf.freebayes.HQ.1.Ae,info,stringsAsFactors=FALSE)
head(vcf.freebayes.HQ.1.Ae)
dim(vcf.freebayes.HQ.1.Ae) # 442050     60 
vcf.freebayes.HQ.1.reform.Ae <- vcf.freebayes.HQ.1.Ae

vcf.freebayes.HQ.1.reform.Ae$MAF_Ae <- round(pmin(vcf.freebayes.HQ.1.reform.Ae$Ae_alt.depth, vcf.freebayes.HQ.1.reform.Ae$Ae_ref.depth)/vcf.freebayes.HQ.1.reform.Ae$Ae_tot.depth, digits = 2) 

hist(vcf.freebayes.HQ.1.reform.Ae$MAF_Ae)

vcf.freebayes.HQ.1.reform.Ae$RPR <- as.numeric(vcf.freebayes.HQ.1.reform.Ae$RPR)
vcf.freebayes.HQ.1.reform.Ae$RPL <- as.numeric(vcf.freebayes.HQ.1.reform.Ae$RPL)

## tail bias 
hist(log10(vcf.freebayes.HQ.1.reform.Ae$RPR))
hist(log10(vcf.freebayes.HQ.1.reform.Ae$RPL))
min(vcf.freebayes.HQ.1.reform.Ae$RPR)
min(vcf.freebayes.HQ.1.reform.Ae$RPL)
max(vcf.freebayes.HQ.1.reform.Ae$RPR + vcf.freebayes.HQ.1.reform.Ae$RPL)

vcf.freebayes.HQ.1.reform.Ae$MBR <- round(pmin(vcf.freebayes.HQ.1.reform.Ae$RPR, vcf.freebayes.HQ.1.reform.Ae$RPL)/(vcf.freebayes.HQ.1.reform.Ae$RPR + vcf.freebayes.HQ.1.reform.Ae$RPL), digits = 2) 


hist(vcf.freebayes.HQ.1.reform.Ae$MBR, breaks = 50)

png("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/MBR_Ae.png", width=4, height=3, units="in", res=300)
hist(vcf.freebayes.HQ.1.reform.Ae$MBR, breaks = 50, xlab = "minor reads balance ratio (MBR)", ylab="number of SNPs", main="")
abline(v=c(0.01), col="red")
dev.off() 

# cor(vcf.freebayes.HQ.1.reform.Ae$MBR, vcf.freebayes.HQ.1.reform.Ae$MAF_Ae)
# vcf.freebayes.HQ.1.reform.Ae.het <- vcf.freebayes.HQ.1.reform.Ae[vcf.freebayes.HQ.1.reform.Ae$Ae_gt=="1/1",]
# 
# cor(vcf.freebayes.HQ.1.reform.Ae.het$MBR, vcf.freebayes.HQ.1.reform.Ae.het$MAF_Ae)

vcf.freebayes.HQ.2.tmp.Ae <- vcf.freebayes.HQ.1.reform.Ae

vcf.freebayes.HQ.2.test1.Ae <- 
vcf.freebayes.HQ.2.tmp.Ae[(vcf.freebayes.HQ.2.tmp.Ae$MBR>0.01),]

# vcf.freebayes.HQ.2.test2 <- 
# vcf.freebayes.HQ.2.tmp[(vcf.freebayes.HQ.2.tmp$RPR>=0 & 
#                            vcf.freebayes.HQ.2.tmp$RPL>0),]

min(as.numeric(vcf.freebayes.HQ.2.test1.Ae$RPL))
min(as.numeric(vcf.freebayes.HQ.2.test1.Ae$RPR))

table(vcf.freebayes.HQ.2.test1.Ae$Ae_gt) 
#     0/1    1/1 
# 127827 187049 

head(vcf.freebayes.HQ.1.reform.Ae) 
save(vcf.freebayes.HQ.1.reform.Ae, file = "~/SNP_parent/result/freebayes/from_cabernet_call_separately/output/vcf.freebayes.HQ.1.reform.Ae.Rdata")

### het problem 
### compare QUAL, read depth, between 0/1 and 1/1

head(vcf.freebayes.HQ.1.Ae)

min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.reform.Ae$Ae_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ae$Ae_gt=="1/1",]$QUAL)

min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ae$Ae_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ae$Ae_gt=="1/1",]$QUAL)

min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ol$Ol_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.1.Ae[vcf.freebayes.HQ.1.Ol$Ol_gt=="1/1",]$QUAL)

min(vcf.freebayes.HQ.2.test1.Ol[vcf.freebayes.HQ.2.test1.Ol$Ol_gt=="0/1",]$QUAL)
min(vcf.freebayes.HQ.2.test1.Ol[vcf.freebayes.HQ.2.test1.Ol$Ol_gt=="1/1",]$QUAL) 

library(ggplot2) 
p.QUAL.Ol <- ggplot(data = vcf.freebayes.HQ.2.test1.Ol) 
p.QUAL.Ol <- p.QUAL.Ol + geom_histogram(aes(log10(QUAL), fill=Ol_gt)) 
p.QUAL.Ol <- p.QUAL.Ol + facet_wrap(~Ol_gt, ncol = 1) 
p.QUAL.Ol <- p.QUAL.Ol +labs(list(x="log10(QUAL)", y="Number of SNPs"))  
p.QUAL.Ol <- p.QUAL.Ol + theme(legend.position = "none") 
p.QUAL.Ol 

# plot QUAL by Ae_gt  
library(ggplot2) 
p.QUAL.Ae <- ggplot(data = vcf.freebayes.HQ.Ae) 
p.QUAL.Ae <- p.QUAL.Ae + geom_histogram(aes(log10(QUAL), fill=Ae_gt)) 
p.QUAL.Ae <- p.QUAL.Ae + facet_wrap(~Ae_gt, ncol = 1) 
p.QUAL.Ae <- p.QUAL.Ae +labs(list(x="log10(QUAL)", y="Number of SNPs"))  
p.QUAL.Ae <- p.QUAL.Ae + theme(legend.position = "none")
p.QUAL.Ae 

ggsave(p.QUAL.Ae, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/p.QUAL.Ae.png", height = 4, width = 6)

# plot total depth by Ae_gt 
p.tot.depth.Ae <- ggplot(data = vcf.freebayes.HQ.2.test1.Ae) 
p.tot.depth.Ae <- p.tot.depth.Ae + geom_histogram(aes(Ae_tot.depth, fill=Ae_gt)) + scale_x_log10() 
p.tot.depth.Ae <- p.tot.depth.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.tot.depth.Ae

# plot MBR
p.MBR.Ae <- ggplot(data = vcf.freebayes.HQ.1.reform.Ae) 
p.MBR.Ae <- p.MBR.Ae + geom_histogram(aes(MBR, fill=Ae_gt))  
p.MBR.Ae <- p.MBR.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.MBR.Ae  
ggsave(p.MBR.Ae, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/p.MBR.Ae.png", height = 4, width = 6)

# alt/ref allele ratio
vcf.freebayes.HQ.2.test1.Ae.het <- 
vcf.freebayes.HQ.2.test1.Ae[vcf.freebayes.HQ.2.test1.Ae$Ae_gt=="0/1",]

vcf.freebayes.HQ.2.test1.Ae.het$allele_ratio <- round(vcf.freebayes.HQ.2.test1.Ae.het$Ae_alt.depth/vcf.freebayes.HQ.2.test1.Ae.het$Ae_ref.depth, digits = 2)

hist(log10(vcf.freebayes.HQ.2.test1.Ae.het$allele_ratio), breaks= 150) 


nrow(vcf.freebayes.HQ.2.test1.Ae.het[vcf.freebayes.HQ.2.test1.Ae.het$allele_ratio>0,])
nrow(vcf.freebayes.HQ.2.test1.Ae.het) 

pl <- ggplot(data = vcf.freebayes.HQ.2.test1.Ae.het) 
pl <- pl + geom_histogram(aes(allele_ratio)) 
pl <- pl +labs(list(x="alt.depth/ref.depth", y="nunber of SNPs"))
pl

# didn't see correlation between MAF & MBR... 
cor.test(vcf.freebayes.HQ.2.test1.Ae.het)
cor.test(vcf.freebayes.HQ.2.test1.Ae.het$MAF_Ae, vcf.freebayes.HQ.2.test1.Ae.het$MBR)
cor.test(vcf.freebayes.HQ.2.test1.Ae$MAF_Ae, vcf.freebayes.HQ.2.test1.Ae$MBR)

cor.test(vcf.freebayes.HQ.2.test1.Ae.het$MAF_Ae, vcf.freebayes.HQ.2.test1.Ae.het$QUAL)
cor.test(vcf.freebayes.HQ.2.test1.Ae$MAF_Ae, vcf.freebayes.HQ.2.test1.Ae$QUAL) 

hist(vcf.freebayes.HQ.2.test1.Ae.het$Ae_alt.depth/vcf.freebayes.HQ.2.test1.Ae.het$Ae_ref.depth, breaks = 50) 
```


# using freebayes to call SNP (now pretty sure is from no-bolting uniquely mapped result)
# import data 
```{r}
source("/Users/ruijuanli/Desktop/Brassica_project/KIAT_RNA_seq/analysis/function_BnRNAseq.R")
vcf.freebayes <- read.table("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/as_diploid/no_bolting_unique_mapped/Ae_Ol_modified_SNPonly_biallelic.recode.ann.vcf",as.is=T,na.strings = ".:.:.:.:.:.:.") 
dim(vcf.freebayes) # 1678792      11
vcf.freebayes$V11 
# save(vcf.freebayes, file ="~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.Rdata")
```

# analysis 
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/vcf.freebayes.Rdata")
# get header 
vcf.header.freebayes <- system("grep '#C' ~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/as_diploid/no_bolting_unique_mapped/Ae_Ol_modified_SNPonly_biallelic.recode.ann.vcf",intern = TRUE) 
vcf.header.freebayes <- sub("#","",vcf.header.freebayes) #get rid of the pound sign
vcf.header.freebayes <- unlist(strsplit(vcf.header.freebayes,split="\t"))

# use function to reformat vcf file 
vcf.freebayes.reform <- SNP.freebayes.reformat(vcf.freebayes, vcf.header.freebayes)
dim(vcf.freebayes.reform) # 1678792      25   
table(vcf.freebayes.reform$Ae_gt)
#    0/0    0/1    1/1     NA 
# 423018 531949 522721 201104 

table(vcf.freebayes.reform$Ol_gt)
#    0/0    0/1    1/1     NA 
# 459651 517580 518072 183489 

# save(vcf.freebayes.reform, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.reform.Rdata") 
```

# 1st round of filtering 
```{r}
#########
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.reform.Rdata")

### 1) filter based on QUAl score 

# subset the data to keep positions where the quality score is 40 or higher 
vcf.freebayes.HQ <- vcf.freebayes.reform[vcf.freebayes.reform$QUAL>40,] # QUAL score of 20, which means less than 0.01 probability that it isn't polymorphic  
dim(vcf.freebayes.HQ) # 1174345      25  
nrow(vcf.freebayes.HQ) / nrow(vcf.freebayes.reform) # 70% of SNPs were retained with QUAL > 40

# count the number 
table(vcf.freebayes.HQ$Ae_gt)
#    0/0    0/1    1/1     NA 
# 229905 314228 473435 156777

table(vcf.freebayes.HQ$Ol_gt) 
#    0/0    0/1    1/1     NA 
# 253687 300596 474418 145644  

# head(vcf.Ae.Ol.no.bolting.unique.HQ[which(vcf.Ae.Ol.no.bolting.unique.HQ$Ae_gt=="0/1"),], 20)

# tmp.unique <- ftable(vcf.Ae.Ol.no.bolting.unique.HQ[,c("Ae_gt","Ol_gt")])
# head(tmp.unique) 
# write.table(vcf.Ae.Ol.no.bolting.unique.HQ, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/output/gt_Ae_Ol.vcf") 
# save(vcf.freebayes.HQ, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.Rdata")
```

# check why there are so many heterozygous, plot by read depth & QUAL score 
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.Rdata")
# 1) Julin suggested plot by mapping score, total depth, and QUAL score... 
# 2) also visualize SNPs in igv to spot the realignment problem, SNP calls from end of reads & intron-exon insertion is bad, don’t believe… 
# 3) 3, 4, and 5 visualize them in IGV…  

# vcf.freebayes.HQ.Ae.het <- vcf.freebayes.HQ[vcf.freebayes.HQ$Ae_gt=="0/1",]
# vcf.freebayes.HQ

head(vcf.freebayes.HQ)

# plot QUAL by Ae_gt 
library(ggplot2)
p.QUAL.Ae <- ggplot(data = vcf.freebayes.HQ) 
p.QUAL.Ae <- p.QUAL.Ae + geom_histogram(aes(QUAL, fill=Ae_gt)) + scale_x_log10() 
p.QUAL.Ae <- p.QUAL.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.QUAL.Ae

# plot QUAL by Ol_gt 
p.QUAL.Ol <- ggplot(data = vcf.freebayes.HQ) 
p.QUAL.Ol <- p.QUAL.Ol + geom_histogram(aes(QUAL, fill=Ol_gt)) + scale_x_log10() 
p.QUAL.Ol <- p.QUAL.Ol + facet_wrap(~Ol_gt, ncol = 1)
p.QUAL.Ol

# plot total depth by Ae_gt 
p.tot.depth.Ae <- ggplot(data = vcf.freebayes.HQ) 
p.tot.depth.Ae <- p.tot.depth.Ae + geom_histogram(aes(Ae_tot.depth, fill=Ae_gt)) + scale_x_log10() 
p.tot.depth.Ae <- p.tot.depth.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.tot.depth.Ae

# plot total depth by Ol_gt 
p.tot.depth.Ol <- ggplot(data = vcf.freebayes.HQ) 
p.tot.depth.Ol <- p.tot.depth.Ol + geom_histogram(aes(Ol_tot.depth, fill=Ol_gt)) + scale_x_log10() 
p.tot.depth.Ol <- p.tot.depth.Ol + facet_wrap(~Ol_gt, ncol = 1)
p.tot.depth.Ol

# no difference in mapping score, because we used only mapping score of 255 for SNP calling. mapping score of 255 means unique mapping, check stats didn't tell me anything... 

#### export vcf file to igv for visulization 
# dim(vcf.Ae.Ol.no.bolting.unique.HQ) # 1502350      25 
# vcf.Ae.Ol.no.bolting.unique.HQ.export <- vcf.Ae.Ol.no.bolting.unique.HQ[,c(1:11)]
# 
# write.table(vcf.Ae.Ol.no.bolting.unique.HQ.export, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.no.bolting.unique.HQ.vcf", row.names = F)

# conclusion: 2,3,4,5 are caused by re-alignment problem, they are not SNPs
```

# continue explore why many 0/1 in my inbred lines (arround 30%), already filtered by QUAL 
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.Rdata")

# 1) filter by max depth 
mean(vcf.freebayes.HQ$Ae_tot.depth, na.rm=T) # 159 
mean(vcf.freebayes.HQ$Ol_tot.depth, na.rm=T) # 156 

# distribution of depth 
?hist
hist(log10(vcf.freebayes.HQ$Ae_tot.depth))
hist(log10(vcf.freebayes.HQ$Ol_tot.depth)) 

# use 1000 as the max, 10 as min 
vcf.freebayes.HQ.1 <- vcf.freebayes.HQ[(!is.na(vcf.freebayes.HQ$Ae_tot.depth) & 
                              !is.na(vcf.freebayes.HQ$Ol_tot.depth) & 
                              vcf.freebayes.HQ$Ae_tot.depth > 10 & 
                              vcf.freebayes.HQ$Ol_tot.depth > 10 & 
                              vcf.freebayes.HQ$Ae_tot.depth < 1000 & 
                              vcf.freebayes.HQ$Ol_tot.depth < 1000),]
dim(vcf.freebayes.HQ.1) # 493019     25 
head(vcf.freebayes.HQ.1)  

####### 
# 2) filter by genotype quality (figure genotype quality in freebayes vcf file)
# ... need to figure this out... 

table(vcf.freebayes.HQ.1$Ae_gt)
table(vcf.freebayes.HQ.1$Ol_gt)

# function to draw plot 
pl.allele.ratio.Ae <- function(vcf.subset){
  vcf.subset$tmp <- with(vcf.subset,Ae_alt.depth/Ae_ref.depth)
  pl <- ggplot(data = vcf.subset) 
  pl <- pl + geom_histogram(aes(x=tmp)) 
  pl <- pl + geom_histogram(aes(Ae_alt.depth/Ae_ref.depth)) + scale_x_log10() 
  # pl <- pl +labs(list(x="log10(alt.depth/ref.depth)"))  + xlim(-1,5)
  return(pl)  
}   

pl.allele.ratio.Ol <- function(vcf.subset){
  pl <- ggplot(data = vcf.subset)
  pl <- pl + geom_histogram(aes(Ol_alt.depth/Ol_ref.depth)) + scale_x_log10()
  pl <- pl +labs(list(x="log10(alt.depth/ref.depth)"))
  return(pl)
} 

pl <- ggplot(data = vcf.GATK.HQ.2.Ae.het)
pl <- pl + geom_histogram(aes(Ae_alt.depth/Ae_ref.depth)) 
pl <- pl +labs(list(x="log10(alt.depth/ref.depth)"))
pl  

pl.allele.ratio.Ae(vcf.freebayes.HQ.1.Ae.het)  
pl.allele.ratio.Ol(vcf.freebayes.HQ.1.Ol.het)   

######
pl.allele.ratio.Ae.orginal <- function(vcf.subset){
  vcf.subset$tmp <- with(vcf.subset,Ae_alt.depth/Ae_ref.depth)
  pl <- ggplot(data = vcf.subset) 
  pl <- pl + geom_histogram(aes(x=tmp)) 
  pl <- pl + geom_histogram(aes(Ae_alt.depth/Ae_ref.depth)) 
  pl <- pl +labs(list(x="(alt.depth/ref.depth)"))  + xlim(0,3) + ylim(0,40000)
  return(pl)  
}   
  
allele.ratio.Ae.het <- pl.allele.ratio.Ae.orginal(vcf.freebayes.HQ.1.Ae.het)
allele.ratio <- pl.allele.ratio.Ae.orginal(vcf.freebayes.HQ.1)

allele.ratio.Ae.het
allele.ratio

# ggsave(allele.ratio.Ae.het, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/allele_ratio_Ae_het.png", height = 8, width = 11)
# ggsave(allele.ratio, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/allele_ratio.png", height = 8, width = 11)

# very similar to what we see from the review paper... 

save(vcf.freebayes.HQ.1, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.1.Rdata")

####### 

# a lot of het calls for Ae & Ol, also there is skewed distribution for the allele depth ratio, is this a result of allele specific expression? If yes, SNPs from the same gene will have skewness. So to test this hypothesis, check the correlation of SNP position & allele depth ratio. 
library(IRanges)
library(GenomicRanges)

# gff3 file 
library(GenomicFeatures)
library("rtracklayer")
# get only mRNA from gff3: cat Brassica_napus.annotation_v5.gff3 | awk '$3=="mRNA"' > Brassica_napus.annotation_v5_mRNAonly.gff3 

gffRangeData <- import.gff("~/Desktop/Brassica_project/IGV_file/Brassica_napus.annotation_v5.gff3") 
gff<-as(gffRangeData, "GRanges")

length(gff)
head(gff, 3)

# get only mRNA 
gff.mRNA <- gff[elementMetadata(gff)[,2] == "mRNA"] 
?disjoin # ordered by range position
gff.mRNA.df <- as.data.frame(disjoin(gff.mRNA))
dim(gff.mRNA.df) # 101040 5
head(gff.mRNA.df)
gff.mRNA.df.chrA01 <- gff.mRNA.df[gff.mRNA.df$seqnames=="chrA01",] 
head(gff.mRNA.df.chrA01)
tail(gff.mRNA.df.chrA01)

gff.mRNA.df$Name <- elementMetadata(gff[elementMetadata(gff)[,2] == "mRNA"])$Name
head(gff.mRNA.df)

#######
# something wrong here .... disjoin gives me sorted data based on ranges, whereas metadata was not sorted...
# combing the two dataset gave me wrong result... 
#######

# vcf.freebayes.HQ.1.Ae.het.tmp <- vcf.freebayes.HQ.1.Ae.het[c("CHROM", "POS", "Ae_gt", "Ae_ref.depth", "Ae_alt.depth")]
# head(vcf.freebayes.HQ.1.Ae.het.tmp)
# 
# vcf.freebayes.HQ.1.Ae.het.tmp.chrA01 <- 
# vcf.freebayes.HQ.1.Ae.het.tmp[vcf.freebayes.HQ.1.Ae.het.tmp$CHROM=="chrA01",]
# dim(vcf.freebayes.HQ.1.Ae.het.tmp[vcf.freebayes.HQ.1.Ae.het.tmp$CHROM=="chrA01",])
# gff.mRNA.df.chrA01 <- 
# gff.mRNA.df[gff.mRNA.df$seqnames=="chrA01",]
# 
# # assign gene name to each SNP based on range... 
# for (i in 1:nrow(vcf.freebayes.HQ.1.Ae.het.tmp.chrA01)){
#   for (j in 1:nrow(gff.mRNA.df.chrA01)){
#   if (
#       (vcf.freebayes.HQ.1.Ae.het.tmp.chrA01$CHROM[i]==gff.mRNA.df.chrA01$seqnames[j]) &  
#       (vcf.freebayes.HQ.1.Ae.het.tmp.chrA01$POS[i] <= gff.mRNA.df.chrA01$end[j]) & 
#       (vcf.freebayes.HQ.1.Ae.het.tmp.chrA01$POS[i] >= gff.mRNA.df.chrA01$start[j])
#       ){
#     vcf.freebayes.HQ.1.Ae.het.tmp.chrA01$Name[i] <- gff.mRNA.df.chrA01$Name[j]
#     }
#   }
# }
# 
# vcf.freebayes.HQ.1.Ae.het.tmp.chrA01
# dim(vcf.freebayes.HQ.1.Ae.het.tmp.chrA01)
# 
# save(vcf.freebayes.HQ.1.Ae.het.tmp.chrA01, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.freebayes.HQ.1.Ae.het.tmp.chrA01.Rdata")
# load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.freebayes.HQ.1.Ae.het.tmp.chrA01.Rdata")
# View(vcf.freebayes.HQ.1.Ae.het.tmp.chrA01)
# vcf.freebayes.HQ.1.Ae.het.tmp.chrA01$allele_ratio <- round(vcf.freebayes.HQ.1.Ae.het.tmp.chrA01$Ae_alt.depth/vcf.freebayes.HQ.1.Ae.het.tmp.chrA01$Ae_ref.depth, digits = 2) 

# only keep 42 columns of INFO, 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.1.Rdata")

tmp.list <- strsplit(vcf.freebayes.HQ.1$INFO,split = ";")
length(tmp.list) # 493019 
length(tmp.list[[1]]) 

# check before modify 
for (i in 1:length(tmp.list)){
    if(length(tmp.list[[i]]) != 42){
      print(i) 
    }
} 

for (i in 1:length(tmp.list)){
    tmp.list[[i]] <- tmp.list[[i]][c(1:42)] 
}

for (i in 1:length(tmp.list)){
    if(length(tmp.list[[i]]) != 42){
      print(i) 
    }
} 

# 1) seperate INFO column  
head(vcf.freebayes.HQ.1)
info <- matrix(
  unlist(tmp.list),
  nrow=nrow(vcf.freebayes.HQ.1),  
  byrow=TRUE
  )

head(info)

colnames(info) <- c("AB","ABP","AC","AF","AN","AO","CIGAR","DP","DPB","DPRA","EPP","EPPR","GTI","LEN","MEANALT",
                    "MQM","MQMR","NS","NUMALT","ODDS","PAIRED","PAIREDR","PAO","PQA","PQR","PRO","QA","QR","RO",
                    "RPL","RPP","RPPR","RPR","RUN","SAF","SAP","SAR","SRF","SRP","SRR","TYPE", "ANN")

head(info)
info <- as.data.frame(info)
colnames(info)

for (i in 1:ncol(info)){
  info[,i] <- gsub("([[:print:]]+)(=)([[:print:]])","\\3",info[,i])
}

head(info)

vcf.freebayes.HQ.1 <- cbind(vcf.freebayes.HQ.1,info,stringsAsFactors=FALSE)
head(vcf.freebayes.HQ.1)
dim(vcf.freebayes.HQ.1) # 493019     67 
tail(vcf.freebayes.HQ.1)
vcf.freebayes.HQ.1.reform <- vcf.freebayes.HQ.1
save(vcf.freebayes.HQ.1.reform, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.1.reform.Rdata")

write.csv(vcf.freebayes.HQ.1$ANN, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/ANN.csv")
# cat ANN.csv | tail -493019 | awk 'BEGIN{FS="|"}{print $4}' | sed 's/^$/intergenic/g' > ANN_gene_ID

ANN_gene_ID <- read.table("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/ANN_gene_ID")

head(ANN_gene_ID)
tail(ANN_gene_ID) 
dim(ANN_gene_ID) # 493019 
class(ANN_gene_ID)
ANN_gene_ID$V1

vcf.freebayes.HQ.1.reform$gene_ID <- ANN_gene_ID$V1
length(vcf.freebayes.HQ.1.reform$gene_ID)
head(vcf.freebayes.HQ.1.reform)
tail(vcf.freebayes.HQ.1.reform)
dim(vcf.freebayes.HQ.1.reform) # 493019     68 

# save(vcf.freebayes.HQ.1.reform, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.1.reform.Rdata") 
```

### check correlation between gene ID and allele ratio, this is the way to check allele specific expression 
```{r} 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.1.reform.Rdata")

# get MAF 
vcf.freebayes.HQ.1.reform$MAF_Ae <- round(pmin(vcf.freebayes.HQ.1.reform$Ae_alt.depth, vcf.freebayes.HQ.1.reform$Ae_ref.depth)/vcf.freebayes.HQ.1.reform$Ae_tot.depth, digits = 2) 

vcf.freebayes.HQ.1.reform$MAF_Ol <- round(pmin(vcf.freebayes.HQ.1.reform$Ol_alt.depth, vcf.freebayes.HQ.1.reform$Ol_ref.depth)/vcf.freebayes.HQ.1.reform$Ol_tot.depth, digits = 2) 

hist(vcf.freebayes.HQ.1.reform$MAF_Ae)
hist(vcf.freebayes.HQ.1.reform$MAF_Ol)

# extract specific columns for better visualization   
colnames(vcf.freebayes.HQ.1.reform)
for_allele_specific_expression <- vcf.freebayes.HQ.1.reform[,c("CHROM", "POS", "gene_ID", "MAF_Ae", "Ae_tot.depth", "MAF_Ol", "Ol_tot.depth")] 
View(for_allele_specific_expression) 

for_allele_specific_expression.Ae.het <- 
  vcf.freebayes.HQ.1.reform[vcf.freebayes.HQ.1.reform$Ae_gt=="0/1",][,c("CHROM", "POS", "gene_ID", "MAF_Ae", "Ae_tot.depth")]

for_allele_specific_expression.Ol.het <- 
  vcf.freebayes.HQ.1.reform[vcf.freebayes.HQ.1.reform$Ol_gt=="0/1",][,c("CHROM", "POS", "gene_ID", "MAF_Ol", "Ol_tot.depth")]

for_allele_specific_expression.Ae.homo <- 
  vcf.freebayes.HQ.1.reform[vcf.freebayes.HQ.1.reform$Ae_gt!="0/1",][,c("CHROM", "POS", "gene_ID", "MAF_Ae", "Ae_tot.depth")]

for_allele_specific_expression.Ol.homo <- 
  vcf.freebayes.HQ.1.reform[vcf.freebayes.HQ.1.reform$Ol_gt!="0/1",][,c("CHROM", "POS", "gene_ID", "MAF_Ol", "Ol_tot.depth")]

hist(for_allele_specific_expression.Ae.het$MAF_Ae)
hist(for_allele_specific_expression.Ae.homo$MAF_Ae)
hist(for_allele_specific_expression.Ol.homo$MAF_Ol)

min(for_allele_specific_expression.Ae.het$MAF_Ae)
min(for_allele_specific_expression.Ol.het$MAF_Ol)
max(for_allele_specific_expression.Ae.homo$MAF_Ae)
max(for_allele_specific_expression.Ol.homo$MAF_Ol)

save(for_allele_specific_expression, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/for_allele_specific_expression.Rdata") 

save(for_allele_specific_expression, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/for_allele_specific_expression_Ae.het.Rdata")

save(for_allele_specific_expression, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/for_allele_specific_expression.Ol.het.Rdata") 

# check by visulization doesn't work.... need some kind of statistical test... 
### one way ANOVA OR auto-correlation to check allele specific expression 
### one way ANOVA, null hypothesis, MAF are the same across all genes? one way ANOVA doesn't work because the dataset is too big... 
class(for_allele_specific_expression.Ae.het$CHROM)
class(for_allele_specific_expression.Ae.het$gene_ID)

# fit <- aov(for_allele_specific_expression.Ae.het$MAF_Ae ~ for_allele_specific_expression.Ae.het$gene_ID) 

# for_ANOVA <- 
# aggregate(for_allele_specific_expression.Ae.het$MAF_Ae, 
#           list(for_allele_specific_expression.Ae.het$gene_ID),
#           mean) 

### check auto-correlation.... leave for later... 

```

# more sophisticated filtering (SNP tail bias, position bias, strand bias, genotype quality...)
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.1.reform.Rdata")
# based on tail bias and strand bias 
vcf.freebayes.HQ.1.reform$SAF <- as.numeric(vcf.freebayes.HQ.1.reform$SAF)
vcf.freebayes.HQ.1.reform$SAR <- as.numeric(vcf.freebayes.HQ.1.reform$SAR)
vcf.freebayes.HQ.1.reform$RPR <- as.numeric(vcf.freebayes.HQ.1.reform$RPR)
vcf.freebayes.HQ.1.reform$RPL <- as.numeric(vcf.freebayes.HQ.1.reform$RPL)

# strand bias, present on both strands
head(vcf.freebayes.HQ.1.reform)

hist(log10(vcf.freebayes.HQ.1.reform$SAF))    
hist(log10(vcf.freebayes.HQ.1.reform$SAR))

# don't need to filter based on strand bias because 
vcf.freebayes.HQ.2.tmp <- 
vcf.freebayes.HQ.1.reform[(vcf.freebayes.HQ.1.reform$SAF>0 & 
                           vcf.freebayes.HQ.1.reform$SAR>0) ,] 

min(vcf.freebayes.HQ.1.reform$SAR); min(vcf.freebayes.HQ.1.reform$SAF)
min(vcf.freebayes.HQ.2$SAR); min(vcf.freebayes.HQ.2$SAF)

dim(vcf.freebayes.HQ.2.tmp) # 460336     68
dim(vcf.freebayes.HQ.1.reform) # 493019     68 

## tail bias 
hist(log10(vcf.freebayes.HQ.2.tmp$RPR))
hist(log10(vcf.freebayes.HQ.2.tmp$RPL))
min(vcf.freebayes.HQ.2.tmp$RPR)
min(vcf.freebayes.HQ.2.tmp$RPL)

vcf.freebayes.HQ.2 <- 
vcf.freebayes.HQ.2.tmp[(vcf.freebayes.HQ.2.tmp$RPR>1 & 
                           vcf.freebayes.HQ.2.tmp$RPL>1),]  


vcf.freebayes.HQ.2.test1 <- 
vcf.freebayes.HQ.2.tmp[(vcf.freebayes.HQ.2.tmp$RPR>0 & 
                           vcf.freebayes.HQ.2.tmp$RPL>=0),]

vcf.freebayes.HQ.2.test2 <- 
vcf.freebayes.HQ.2.tmp[(vcf.freebayes.HQ.2.tmp$RPR>=0 & 
                           vcf.freebayes.HQ.2.tmp$RPL>0),]


min(as.numeric(vcf.freebayes.HQ.2$RPL))
min(as.numeric(vcf.freebayes.HQ.2$RPR))

min(as.numeric(vcf.freebayes.HQ.2.test1$RPL))
min(as.numeric(vcf.freebayes.HQ.2.test1$RPR))

min(as.numeric(vcf.freebayes.HQ.2.test2$RPL))
min(as.numeric(vcf.freebayes.HQ.2.test2$RPR))

########
# need to do one way ANOVA OR auto-correlation to test for correlation between gene & MAF... 
########

dim(vcf.freebayes.HQ.2) # 354820     68 

table(vcf.freebayes.HQ.2$Ae_gt) 
 #   0/0    0/1    1/1 
 # 96721 112112 145987 

table(vcf.freebayes.HQ.2$Ol_gt)
#    0/0    0/1    1/1 
# 100508 109057 145255

# find filtered SNPs because of tail bias, and visulize them in IGV 
test1 <- vcf.freebayes.HQ.2.tmp[,c("CHROM", "POS")]
test2 <- vcf.freebayes.HQ.2[,c("CHROM", "POS")]
library(dplyr)
filtered_out_snp_tail_bias <- anti_join(test1, test2)
head(filtered_out_snp_tail_bias[filtered_out_snp_tail_bias$CHROM=="chrA01",], 30)

save(vcf.freebayes.HQ.2, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.2.Rdata")
```

# correlation between minor read balance ratio & MAF, because the assumption is that read balance ratio measures the tail bias, the small the ratio is, the more probable that SNPs from that site is not a true SNP,  
```{r} 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.2.Rdata")
head(vcf.freebayes.HQ.2) 


```

# export the filtered dataset as VCF for IGV visualization
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.2.Rdata")
 
head(vcf.freebayes.HQ.2)
colnames(vcf.freebayes.reform) 
tmp <- c("CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT", "Ae", "Ol")

vcf.freebayes.HQ.2.reform <- vcf.freebayes.HQ.2[,tmp]
colnames(vcf.freebayes.HQ.2.reform)
head(vcf.freebayes.HQ.2.reform) 

write.table(vcf.freebayes.HQ.2.reform, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.2.reform.txt")

# get header from original vcf file 
# run ./reconstruct_vcf.sh (https://github.com/leejimmy93/KIAT_whitney/tree/master/parent_SNP) 
# import into IGV for better visualization... 

```

# SNPs between Da-Ae & Da-Ol-1 
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.1.reform.Rdata")

vcf.tmp <- vcf.freebayes.HQ.1.reform # assign a variable so I don't have to change the name every time... 

dim(vcf.tmp) # 493019     68  

# SNPs between Ae & Ol 
sum(vcf.tmp$Ae_gt=="1/1" & vcf.tmp$Ol_gt=="0/0") +  sum(vcf.tmp$Ae_gt=="0/0" & vcf.tmp$Ol_gt=="1/1") # 137069 SNPs between Ae & Ol 

(sum(vcf.tmp$Ae_gt=="1/1" & vcf.tmp$Ol_gt=="0/0") +  sum(vcf.tmp$Ae_gt=="0/0" & vcf.tmp$Ol_gt=="1/1"))/length(vcf.freebayes.HQ.1.reform$CHROM) # 27% of the total SNPs 
# subset SNPs between Ae & Ol 
vcf.Ae.Ol.freebayes <- vcf.tmp[((vcf.tmp$Ae_gt=="1/1" & vcf.tmp$Ol_gt=="0/0") | (vcf.tmp$Ae_gt=="0/0" & vcf.tmp$Ol_gt=="1/1")),] 

dim(vcf.Ae.Ol.freebayes) # 137069     68 
table(vcf.Ae.Ol.freebayes$Ae_gt)
table(vcf.Ae.Ol.freebayes$Ol_gt)

save(vcf.Ae.Ol.freebayes, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.Ae.Ol.freebayes.Rdata")

######### using IGV to check SNP 
# tmp <- c("CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT", "Ae", "Ol")
# 
# vcf.Ae.Ol.freebayes <- vcf.Ae.Ol.freebayes[,tmp]
# colnames(vcf.freebayes.HQ.2.reform)
# head(vcf.freebayes.HQ.2.reform) 
# 
# write.table(vcf.Ae.Ol.freebayes, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.Ae.Ol.freebayes.txt")

# scp to whitney & run it on whitney
# scp vcf.Ae.Ol.freebayes.txt ruijuanli@whitney.plb.ucdavis.edu:/Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/result/freebayes 
# get header from original vcf file 

# run ./reconstruct_vcf.sh (https://github.com/leejimmy93/KIAT_whitney/tree/master/parent_SNP) 
# import into IGV for visualization... 

head(vcf.Ae.Ol.freebayes[,c("CHROM", "POS")], 50)
```


# stats on SNPs between Da-Ae & Da-Ol-1 
```{r}
# number of SNPs on per chromosome # subgenome group/chr ID/main or random chromosome 
chr_ID <- unique(Ae_Ol_after_filter.unique$CHROM)
number.unique <- list() 

for (chr in chr_ID) {
  number.unique[chr] <- sum(sum(Ae_Ol_after_filter.unique$CHROM==chr))
  print(number.unique[chr])
}

SNP_chr.unique <- as.data.frame(t(as.data.frame(number.unique)))
SNP_chr.unique$ID <- rownames(SNP_chr.unique)
head(SNP_chr.unique)
# split to main & random chromosomes 
# main 
SNP_chr_main.unique <- SNP_chr.unique[!grepl("random", SNP_chr.unique$ID),]  
SNP_chr_main.unique$subgenome <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\2", SNP_chr_main.unique$ID)   
SNP_chr_main.unique$chr_ID <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\3", SNP_chr_main.unique$ID) 

# random 
SNP_chr_random.unique <- SNP_chr.unique[grep("random", SNP_chr.unique$ID),]  
SNP_chr_random.unique$subgenome <- gsub("(chr)(A|C|U)(01|02|03|04|05|06|07|08|09|10|nn)(_)(random)", "\\2", SNP_chr_random.unique$ID)  
SNP_chr_random.unique$chr_ID <- gsub("(chr)(A|C|U)(01|02|03|04|05|06|07|08|09|10|nn)(_)(random)", "\\3", SNP_chr_random.unique$ID) 

# do a little calculation: how many SNPs on main & random chromosomes? 
sum(SNP_chr_main.unique$V1) # 200317    
sum(SNP_chr_random.unique$V1) # 32160    
sum(SNP_chr_main.unique$V1)/sum(sum(SNP_chr_main.unique$V1),sum(SNP_chr_random.unique$V1)) # 86.2% 

# make a plot to see how many SNPs on each chromosome 
library(ggplot2)
pl.SNP.main.1.unique <- ggplot(data = SNP_chr_main.unique) 
pl.SNP.main.1.unique <- pl.SNP.main.1.unique + geom_bar(aes(x=factor(chr_ID), y = V1, fill=subgenome), stat = "identity")
pl.SNP.main.1.unique <- pl.SNP.main.1.unique + facet_wrap(~subgenome, ncol = 1)
# pl.SNP.main.1 <- pl.SNP.main.1 + geom_text(aes(x=factor(chr_ID), y = number, label=factor(number), size=1)) 
pl.SNP.main.1.unique <- pl.SNP.main.1.unique + labs(list(title = "", x = "chromosome ID", y = "number of SNPs"))
pl.SNP.main.1.unique <- pl.SNP.main.1.unique + theme(legend.position = "none")
pl.SNP.main.1.unique 

ggsave(pl.SNP.main.1.unique, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/output/SNP.main.1.unique.01_26.png", height = 6, width = 8)

########### plot the position along the chromosome of each SNP, read depth 
head(Ae_Ol_after_filter.unique) 
# write.table(Ae_Ol_after_filter.unique, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/output/Ae_Ol_after_filter.unique.vcf")
# get main chromosome and random chromsome seperately 
Ae_Ol_after_filter.main.unique <- Ae_Ol_after_filter.unique[grep("random", Ae_Ol_after_filter.unique$CHROM, invert = T),]
dim(Ae_Ol_after_filter.main.unique) # 200317     25   

Ae_Ol_after_filter.main.unique$subgenome <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\2", Ae_Ol_after_filter.main.unique$CHROM)
Ae_Ol_after_filter.main.unique$chr_ID <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\3", Ae_Ol_after_filter.main.unique$CHROM)

# number of SNPs on A & C subgenome
sum(Ae_Ol_after_filter.main.unique$subgenome=="A")/nrow(Ae_Ol_after_filter.main.unique) # 69.6%
sum(Ae_Ol_after_filter.main.unique$subgenome=="C") # 61464

## make plot for main & random chromosome seperately (distribution of SNPs per Mb)
library(ggplot2)
pl.SNP.main.2.unique <- ggplot(data = Ae_Ol_after_filter.main.unique)
pl.SNP.main.2.unique <- pl.SNP.main.2.unique + geom_histogram(aes(x=POS, fill=subgenome), binwidth = 1000000) 
pl.SNP.main.2.unique <- pl.SNP.main.2.unique + facet_grid(chr_ID ~subgenome)
pl.SNP.main.2.unique <- pl.SNP.main.2.unique + labs(list(title = "", x = "chromosome ID", y = "number of SNPs"))
pl.SNP.main.2.unique <- pl.SNP.main.2.unique + theme(legend.position = "none")
pl.SNP.main.2.unique

ggsave("~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/output/SNP.main.2.unique.01.26.2017.png", width = 23, height = 13)

Ae_Ol_after_filter.random.unique <- Ae_Ol_after_filter.unique[grep("random", Ae_Ol_after_filter.unique$CHROM),] 
dim(Ae_Ol_after_filter.random.unique) # 32160    25   

### check Da-Ae VS ref, do I get the same pattern? 

chr_ID_Ae_ref.unique <- unique(Ae_ref_after_filter.unique$CHROM)
number_Ae_ref.unique <- c()

for (chr in chr_ID_Ae_ref.unique) {
  number_Ae_ref.unique[chr] <- sum(sum(Ae_ref_after_filter.unique$CHROM==chr))
  print(number_Ae_ref.unique[chr])
}

SNP_chr_Ae_ref.unique <- as.data.frame(number_Ae_ref.unique)
SNP_chr_Ae_ref.unique
SNP_chr_Ae_ref.unique$ID <- rownames(SNP_chr_Ae_ref.unique)
# split to main & random chromosomes
# main
SNP_chr_main_Ae_ref.unique <- SNP_chr_Ae_ref.unique[!grepl("random", SNP_chr_Ae_ref.unique$ID),]
SNP_chr_main_Ae_ref.unique$subgenome <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\2", SNP_chr_main_Ae_ref.unique$ID)
SNP_chr_main_Ae_ref.unique$subgenome
SNP_chr_main_Ae_ref.unique$chr_ID <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\3", SNP_chr_main_Ae_ref.unique$ID)
SNP_chr_main_Ae_ref.unique$chr_ID
SNP_chr_main_Ae_ref.unique$number_Ae_ref
# random
SNP_chr_random_Ae_ref.unique <- SNP_chr_Ae_ref.unique[grep("random", SNP_chr_Ae_ref.unique$ID),]
SNP_chr_random_Ae_ref.unique$subgenome <- gsub("(chr)(A|C|U)(01|02|03|04|05|06|07|08|09|10|nn)(_)(random)", "\\2", SNP_chr_random_Ae_ref.unique$ID)
SNP_chr_random_Ae_ref.unique$chr_ID <- gsub("(chr)(A|C|U)(01|02|03|04|05|06|07|08|09|10|nn)(_)(random)", "\\3", SNP_chr_random_Ae_ref.unique$ID)

# do a little calculation: how many SNPs on main & random chromosomes?
sum(SNP_chr_main_Ae_ref.unique$number) # 325903  
sum(SNP_chr_random_Ae_ref.unique$number) # 68834  
sum(SNP_chr_main_Ae_ref.unique$number)/sum(sum(SNP_chr_main_Ae_ref.unique$number),sum(SNP_chr_random_Ae_ref.unique$number)) # 82.6% 

sum(SNP_chr_main_Ae_ref.unique[SNP_chr_main_Ae_ref.unique$subgenome=="A",]$number_Ae_ref) # 203780  
sum(SNP_chr_main_Ae_ref.unique[SNP_chr_main_Ae_ref.unique$subgenome=="C",]$number_Ae_ref) # 122123    

SNP_chr_main_Ae_ref.unique$chr_ID_2 <- paste(SNP_chr_main_Ae_ref.unique$subgenome, SNP_chr_main_Ae_ref.unique$chr_ID, sep = "")
SNP_chr_main_Ae_ref.unique$chr_ID_2
 
# make a plot to see how many SNPs on each chromosome 
pl.SNP.main.Ae_ref.unique <- ggplot(data = SNP_chr_main_Ae_ref.unique) 
pl.SNP.main.Ae_ref.unique <- pl.SNP.main.Ae_ref.unique + geom_bar(aes(x=factor(chr_ID_2), y = number_Ae_ref.unique, fill=subgenome), stat = "identity")
# pl.SNP.main.Ae_ref <- pl.SNP.main.Ae_ref + facet_wrap(~subgenome, ncol = 1)
# pl.SNP.main.1 <- pl.SNP.main.1 + geom_text(aes(x=factor(chr_ID), y = number, label=factor(number), size=1))
pl.SNP.main.Ae_ref.unique <- pl.SNP.main.Ae_ref.unique + labs(list(title = "", x = "chromosome ID", y = "number of SNPs"))
# pl.SNP.main.Ae_ref <- pl.SNP.main.Ae_ref + theme(legend.position = "none")
pl.SNP.main.Ae_ref.unique

# ggsave(pl.SNP.main.Ae_ref, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/output/SNP.main.Ae_ref.png", height = 6, width = 10) 

chr_ID_Ol_ref.unique <- unique(Ol_ref_after_filter.unique$CHROM)
number_Ol_ref.unique <- c()  

for (chr in chr_ID_Ol_ref.unique) {
  number_Ol_ref.unique[chr] <- sum(sum(Ol_ref_after_filter.unique$CHROM==chr))
  print(number_Ol_ref.unique[chr])s
}

SNP_chr_Ol_ref.unique <- as.data.frame(number_Ol_ref.unique)
SNP_chr_Ol_ref.unique
SNP_chr_Ol_ref.unique$ID <- rownames(SNP_chr_Ol_ref.unique)
# split to main & random chromosomes
# main
SNP_chr_main_Ol_ref.unique <- SNP_chr_Ol_ref.unique[!grepl("random", SNP_chr_Ol_ref.unique$ID),]
SNP_chr_main_Ol_ref.unique$subgenome <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\2", SNP_chr_main_Ol_ref.unique$ID)
SNP_chr_main_Ol_ref.unique$subgenome
SNP_chr_main_Ol_ref.unique$chr_ID <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\3", SNP_chr_main_Ol_ref.unique$ID)
SNP_chr_main_Ol_ref.unique$chr_ID
SNP_chr_main_Ol_ref.unique$number_Ol_ref
# random
SNP_chr_random_Ol_ref.unique <- SNP_chr_Ol_ref.unique[grep("random", SNP_chr_Ol_ref.unique$ID),]
SNP_chr_random_Ol_ref.unique$subgenome <- gsub("(chr)(A|C|U)(01|02|03|04|05|06|07|08|09|10|nn)(_)(random)", "\\2", SNP_chr_random_Ol_ref.unique$ID)
SNP_chr_random_Ol_ref.unique$chr_ID <- gsub("(chr)(A|C|U)(01|02|03|04|05|06|07|08|09|10|nn)(_)(random)", "\\3", SNP_chr_random_Ol_ref.unique$ID)

# do a little calculation: how many SNPs on main & random chromosomes?
sum(SNP_chr_main_Ol_ref.unique$number) #322119  
sum(SNP_chr_random_Ol_ref.unique$number) # 69825   

sum(SNP_chr_main_Ol_ref.unique[SNP_chr_main_Ol_ref.unique$subgenome=="A",]$number_Ol_ref) # 189154   
sum(SNP_chr_main_Ol_ref.unique[SNP_chr_main_Ol_ref.unique$subgenome=="C",]$number_Ol_ref) # 132965  

SNP_chr_main_Ol_ref.unique$chr_ID_2 <- paste(SNP_chr_main_Ol_ref.unique$subgenome, SNP_chr_main_Ol_ref.unique$chr_ID, sep = "")

# make a plot to see how many SNPs on each chromosome
pl.SNP.main.Ol_ref.unique <- ggplot(data = SNP_chr_main_Ol_ref.unique)
pl.SNP.main.Ol_ref.unique <- pl.SNP.main.Ol_ref.unique + geom_bar(aes(x=factor(chr_ID_2), y = number_Ol_ref.unique, fill=subgenome), stat = "identity")
# pl.SNP.main.Ol_ref <- pl.SNP.main.Ol_ref + facet_wrap(~subgenome, ncol = 1)
# pl.SNP.main.1 <- pl.SNP.main.1 + geom_text(aes(x=factor(chr_ID), y = number, label=factor(number), size=1))
pl.SNP.main.Ol_ref.unique <- pl.SNP.main.Ol_ref.unique + labs(list(title = "", x = "chromosome ID", y = "number of SNPs"))
# pl.SNP.main.Ol_ref <- pl.SNP.main.Ol_ref + theme(legend.position = "none")
pl.SNP.main.Ol_ref.unique

ggsave(pl.SNP.main.Ol_ref.unique, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/output/SNP.main.Ol_ref.01_26_2017.png", height = 6, width = 10)

# combine 2 graphs
library(cowplot)
plot.w.ref.unique<-plot_grid(
  pl.SNP.main.Ae_ref.unique+labs(title="Da-Ae VS. reference genome"),
  pl.SNP.main.Ol_ref.unique+labs(title="Da-Ol-1 VS. reference genome"),
  ncol=1, nrow = 2,labels=c("","","",""))

plot.w.ref.unique 
save_plot("~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/output/Ae_Ol_ref", plot.w.ref, ncol = 1, nrow = 2,base_aspect_ratio = 0.8)  

```

# using GATK to call SNP (with bolting data) 
# SNP from no bolting data GATK 
```{r}
# 1) prepare for SNP calling using SNP_calling_GATK_pipeline.sh (on github https://github.com/leejimmy93/KIAT_whitney/blob/master/parent_SNP/SNP_calling_GATK_pipeline.sh)  

# 2) call SNPs using SNP_calling_array_no_bolting.slurm (on github https://github.com/leejimmy93/KIAT_cabernet/blob/master/parent_SNP/SNP_calling_array_no_bolting.slurm) 

# 3) SNP filtering based on no more than 3 SNPs within a cluster 35nt 
# SNP_filter_array_02_09_2017.slurm (https://github.com/leejimmy93/KIAT_cabernet/blob/master/parent_SNP/SNP_filter_array_02_09_2017.slurm) 

# 4) combine SNPs from different chromosomes 
# GATK_combine.sh (https://github.com/leejimmy93/KIAT_whitney/blob/master/parent_SNP/GATK_combine.sh)

# 5) VCF format check/validate 
# GATK -T ValidateVariants -R ~/Reference/B.napus/Brassica_napus_v4.1.chromosomes.fa -V Ae.vcf --validationTypeToExclude ALL
 
# 6) extract only biallelic SNPs from vcf files 
# extract_SNP.sh & extract_biallelic.sh (https://github.com/leejimmy93/KIAT_whitney/tree/master/parent_SNP)

# data location: /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/result/GATK/raw_vcf/vcf_split_by_chr/GATK_SNP_biallelic.recode.ann.vcf  

# 7) import into R for analysis 
# SNP.GATK.reformat(vcf, header) can be used to reformat vcf file from GATK 
# SNP.GATK.basic.filter(vcf) can be used to filter based on snpcluster & QUAL, also give some stats on genotyping matrxi result, and filtering stats 
```

# SNP filtering from GATK (filtered based on sliding window & annotated)
```{r}
### import data 
vcf.GATK <- read.table("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/GATK/as_diploid/no_bolting/GATK_SNP_biallelic.recode.ann.vcf",as.is=T,na.strings = ".") 

head(vcf.GATK, 100)
# save(vcf.GATK, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/GATK/as_diploid/no_bolting/vcf.GATK.Rdata")
```

# analysis 
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/GATK/as_diploid/no_bolting/vcf.GATK.Rdata") 
dim(vcf.GATK) # 1375360      11 
vcf.header.GATK <- system("grep '#C' ~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/GATK/as_diploid/no_bolting/GATK_SNP_biallelic.recode.ann.vcf",intern = TRUE) 
vcf.header.GATK <- sub("#","",vcf.header.GATK) #get rid of the pound sign
vcf.header.GATK <- unlist(strsplit(vcf.header.GATK,split="\t"))

vcf.GATK.reform <- SNP.GATK.reformat(vcf.GATK, vcf.header.GATK)
dim(vcf.GATK.reform) # 1375116      21 
head(vcf.GATK.reform)

table(vcf.GATK.reform$Ae_gt)
table(vcf.GATK.reform$Ol_gt)

#### filter (based on QUAL score and snpcluster)
# basic filtering based on snpcluster & QUAL 
vcf.GATK.HQ <- SNP.GATK.basic.filter(vcf.GATK.reform)
dim(vcf.GATK.HQ) # 617141     21   

# save(vcf.GATK.reform, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/GATK/vcf.GATK.reform.Rdata")
# save(vcf.GATK.HQ, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/GATK/vcf.GATK.HQ.Rdata")
``` 

# 2nd round of filtering 
```{r}
# load data 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/GATK/vcf.GATK.HQ.Rdata")

# 1) max read depth, should not be much higher than average read depth 
mean(vcf.GATK.HQ$Ae_approx.depth, na.rm=T) # 80.38611 
mean(vcf.GATK.HQ$Ol_approx.depth, na.rm=T) # 82.96372 

# plot read depth 
hist(log10(vcf.GATK.HQ$Ae_approx.depth))
hist(log10(vcf.GATK.HQ$Ol_approx.depth))

# based on the distribution, choose 10-1000 as the threshold 

# use 10-1000 as the depth range  

vcf.GATK.HQ.1 <- vcf.GATK.HQ[(!is.na(vcf.GATK.HQ$Ae_approx.depth) & 
                              !is.na(vcf.GATK.HQ$Ol_approx.depth) & 
                              vcf.GATK.HQ$Ae_approx.depth > 10 & 
                              vcf.GATK.HQ$Ol_approx.depth > 10 & 
                              vcf.GATK.HQ$Ae_approx.depth < 1000 & 
                              vcf.GATK.HQ$Ol_approx.depth < 1000),]
dim(vcf.GATK.HQ.1) # 252439     21 
unique(vcf.GATK.HQ.1$CHROM)

# 2) genotype quality, greater than 30   
png("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/genotype.quality.png", width=4, height=3, units="in", res=300)
hist(vcf.GATK.HQ.1$Ae_genotype.qual, xlab = "genotype quality", ylab="number of SNPs", main="genotype quality")
abline(v=c(30), col="red")
dev.off() 

hist(vcf.GATK.HQ.1$Ol_genotype.qual) 
head(vcf.GATK.HQ.1)

## figure 
p.gtqual.Ae <- ggplot(data = vcf.GATK.HQ.1) 
p.gtqual.Ae <- p.gtqual.Ae + geom_histogram(aes(Ae_genotype.qual))  
# p.gtqual.Ae <- p.gtqual.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.gtqual.Ae

vcf.GATK.HQ.2 <- vcf.GATK.HQ.1[(!is.na(vcf.GATK.HQ.1$Ae_genotype.qual) &
                                !is.na(vcf.GATK.HQ.1$Ae_genotype.qual) & 
                                vcf.GATK.HQ.1$Ae_genotype.qual > 30 & 
                                vcf.GATK.HQ.1$Ol_genotype.qual > 30),]

dim(vcf.GATK.HQ.2) # 232294     21 
table(vcf.GATK.HQ.2$Ae_gt)
table(vcf.GATK.HQ.2$Ol_gt)
save(vcf.GATK.HQ.2, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/GATK/vcf.GATK.HQ.2.Rdata")

# 3) check MAF 
#### reformat get seperate depth for alt & ref alleles 
Ae.gt <- matrix(
  unlist(strsplit(vcf.GATK.HQ.2$Ae_ref.alt.depth,split = ",")),
  nrow=nrow(vcf.GATK.HQ.2),  
  byrow=TRUE 
  ) 

colnames(Ae.gt) <- paste("Ae",c("ref.depth","alt.depth"),sep="_")

Ol.gt <- matrix(
  unlist(strsplit(vcf.GATK.HQ.2$Ol_ref.alt.depth,split = ",")),
  nrow=nrow(vcf.GATK.HQ.2),  
  byrow=TRUE
  )

colnames(Ol.gt) <- paste("Ol",c("ref.depth","alt.depth"),sep="_")

vcf.GATK.HQ.2 <- cbind(vcf.GATK.HQ.2, Ae.gt, Ol.gt, stringsAsFactors=FALSE)

vcf.GATK.HQ.2$Ae_alt.depth <- as.numeric(vcf.GATK.HQ.2$Ae_alt.depth)
vcf.GATK.HQ.2$Ae_ref.depth <- as.numeric(vcf.GATK.HQ.2$Ae_ref.depth)
vcf.GATK.HQ.2$Ol_alt.depth <- as.numeric(vcf.GATK.HQ.2$Ol_alt.depth)
vcf.GATK.HQ.2$Ol_ref.depth <- as.numeric(vcf.GATK.HQ.2$Ol_ref.depth)

head(vcf.GATK.HQ.2)

# calculate MAF for each allele # total read depth or specific to Da-Ae and Da-Ol-1??? 
vcf.GATK.HQ.2$Ae_MAF <- round(pmin(vcf.GATK.HQ.2$Ae_alt.depth, vcf.GATK.HQ.2$Ae_ref.depth)/vcf.GATK.HQ.2$Ae_approx.depth, digits = 2)

vcf.GATK.HQ.2$Ol_MAF <- round(pmin(vcf.GATK.HQ.2$Ol_alt.depth, vcf.GATK.HQ.2$Ol_ref.depth)/vcf.GATK.HQ.2$Ol_approx.depth, digits = 2)

hist(vcf.GATK.HQ.2$Ae_MAF)
hist(vcf.GATK.HQ.2$Ol_MAF)

vcf.GATK.HQ.2.reform <- vcf.GATK.HQ.2
save(vcf.GATK.HQ.2.reform, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.GATK.HQ.2.reform.Rdata") 
``` 

# figure out why many 0/1 
```{r} 
# plot the distribution of genotype quality for Ae & Ol 
library(ggplot2)
p.gtqual.Ae <- ggplot(data = vcf.GATK.HQ.2) 
p.gtqual.Ae <- p.gtqual.Ae + geom_histogram(aes(Ae_genotype.qual, fill=Ae_gt))  
p.gtqual.Ae <- p.gtqual.Ae + facet_wrap(~Ae_gt, ncol = 1)
p.gtqual.Ae

p.gtqual.Ol <- ggplot(data = vcf.GATK.HQ.2) 
p.gtqual.Ol <- p.gtqual.Ol + geom_histogram(aes(Ol_genotype.qual, fill=Ol_gt))  
p.gtqual.Ol <- p.gtqual.Ol + facet_wrap(~Ol_gt, ncol = 1)
p.gtqual.Ol

head(vcf.GATK.HQ.2) 
dim(vcf.GATK.HQ.2) # 232294     21 
vcf.GATK.HQ.2$Ae_ref.alt.depth 

# reformat to get seperate read depth for Ae & Ol 
Ae.gt <- matrix(
  unlist(strsplit(vcf.GATK.HQ.2$Ae_ref.alt.depth,split = ",")),
  nrow=nrow(vcf.GATK.HQ.2),  
  byrow=TRUE 
  ) 

# check 
# for (i in 1:length(strsplit(vcf.GATK.HQ$Ae_ref.alt.depth,split = ","))){
#   if (length(strsplit(vcf.GATK.HQ$Ae_ref.alt.depth,split = ",")[[i]])!=2){
#   print(i)    
#   }
# }
# 
# for (i in 1:length(strsplit(vcf.GATK.reform$Ae_ref.alt.depth,split = ","))){
#   if (length(strsplit(vcf.GATK.reform$Ae_ref.alt.depth,split = ",")[[i]])!=2){
#   print(i)    
#   }
# }

colnames(Ae.gt) <- paste("Ae",c("ref.depth","alt.depth"),sep="_")

Ol.gt <- matrix(
  unlist(strsplit(vcf.GATK.HQ.2$Ol_ref.alt.depth,split = ",")),
  nrow=nrow(vcf.GATK.HQ.2),  
  byrow=TRUE
  )

colnames(Ol.gt) <- paste("Ol",c("ref.depth","alt.depth"),sep="_")

vcf.GATK.HQ.2.reform <- cbind(vcf.GATK.HQ.2, Ae.gt, Ol.gt, stringsAsFactors=FALSE)

vcf.GATK.HQ.2.reform$Ae_alt.depth <- as.numeric(vcf.GATK.HQ.2.reform$Ae_alt.depth)
vcf.GATK.HQ.2.reform$Ae_ref.depth <- as.numeric(vcf.GATK.HQ.2.reform$Ae_ref.depth)
vcf.GATK.HQ.2.reform$Ol_alt.depth <- as.numeric(vcf.GATK.HQ.2.reform$Ol_alt.depth)
vcf.GATK.HQ.2.reform$Ol_ref.depth <- as.numeric(vcf.GATK.HQ.2.reform$Ol_ref.depth)

head(vcf.GATK.HQ.2.reform)
save(vcf.GATK.HQ.2.reform, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.GATK.HQ.2.reform.Rdata")

# after filtered by max read depth & genotype quality, high het problem still exist 

######## below, work to figure the cause... 

# 1) check the ref/alt allele ratio for the het calls, distribution plot for Ae & Ol 
# get het & homo call for Ae & Ol 
vcf.GATK.HQ.2.Ae.het <- vcf.GATK.HQ.2[(vcf.GATK.HQ.2$Ae_gt=="0/1"),]
sum(vcf.GATK.HQ.2.Ae.het$Ae_alt.depth=="0") 
sum(vcf.GATK.HQ.2.Ae.het$Ae_ref.depth=="0") 
vcf.GATK.HQ.2.Ae.het[vcf.GATK.HQ.2.Ae.het$Ae_ref.depth=="0",]


vcf.GATK.HQ.2.Ae.het[vcf.GATK.HQ.2.Ae.het$Ae_ref.depth=="0",]
# vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ae.ref <- vcf.Ae.Ol.no.bolting.unique.HQ.filtered[(vcf.Ae.Ol.no.bolting.unique.HQ.filtered$Ae_gt=="0/0"),]
# vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ae.alt <- vcf.Ae.Ol.no.bolting.unique.HQ.filtered[(vcf.Ae.Ol.no.bolting.unique.HQ.filtered$Ae_gt=="1/1"),]

# dim(vcf.GATK.reform.HQ.filtered.Ae.het); dim(vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ae.ref); dim(vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ae.alt)

vcf.GATK.HQ.2.Ol.het <- vcf.GATK.HQ.2[(vcf.GATK.HQ.2$Ol_gt=="0/1"),]
sum(vcf.GATK.HQ.2.Ol.het$Ol_alt.depth=="0") 
sum(vcf.GATK.HQ.2.Ol.het$Ol_ref.depth=="0") 
vcf.GATK.HQ.2.Ol.het[vcf.GATK.HQ.2.Ol.het$Ol_ref.depth=="0",]

vcf.GATK.HQ.2.Ae.het[vcf.GATK.HQ.2.Ae.het$Ae_ref.depth=="0",]
# vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ol.ref <- vcf.Ae.Ol.no.bolting.unique.HQ.filtered[(vcf.Ae.Ol.no.bolting.unique.HQ.filtered$Ol_gt=="0/0"),]
# vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ol.alt <- vcf.Ae.Ol.no.bolting.unique.HQ.filtered[(vcf.Ae.Ol.no.bolting.unique.HQ.filtered$Ol_gt=="1/1"),]

# table(vcf.Ae.Ol.no.bolting.unique.HQ.filtered$Ol_gt) 
# dim(vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ol.het); dim(vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ol.ref); dim(vcf.Ae.Ol.no.bolting.unique.HQ.filtered.Ol.alt) 

# function to draw plot 
pl.allele.ratio.Ae <- function(vcf.subset){
  vcf.subset$tmp <- with(vcf.subset,Ae_alt.depth/Ae_ref.depth)
  pl <- ggplot(data = vcf.subset) 
  pl <- pl + geom_histogram(aes(x=tmp)) 
  pl <- pl + geom_histogram(aes(Ae_alt.depth/Ae_ref.depth)) + scale_x_log10() 
  # pl <- pl +labs(list(x="log10(alt.depth/ref.depth)"))  + xlim(-1,5)
  return(pl)  
}   

pl.allele.ratio.Ol <- function(vcf.subset){
  pl <- ggplot(data = vcf.subset)
  pl <- pl + geom_histogram(aes(Ol_alt.depth/Ol_ref.depth)) + scale_x_log10()
  pl <- pl +labs(list(x="log10(alt.depth/ref.depth)"))
  return(pl)
} 

colnames(vcf.GATK.HQ.2.Ae.het)

pl <- ggplot(data = vcf.GATK.HQ.2.Ae.het)
pl <- pl + geom_histogram(aes(Ae_alt.depth/Ae_ref.depth)) 
pl <- pl +labs(list(x="log10(alt.depth/ref.depth)"))
pl  

pl.allele.ratio.Ae(vcf.GATK.HQ.2.Ae.het)  

### zero existed in my het? 
# vcf.GATK.reform.HQ.filtered.Ae.het$Ae_alt.depth==0, 
vcf.GATK.reform.HQ.filtered[vcf.GATK.reform.HQ.filtered$Ae_gt=="0/1",]$Ae_ref.alt.depth

length(which(vcf.GATK.reform.HQ.filtered.Ae.het$Ae_alt.depth!="NA" & vcf.GATK.reform.HQ.filtered.Ae.het$Ae_alt.depth=="0")) / nrow(vcf.GATK.reform.HQ.filtered.Ae.het)

head(vcf.GATK.reform.HQ.filtered.Ae.het)

vcf.GATK.reform.HQ.filtered.Ae.het$Ae_alt.depth

vcf.GATK.reform.HQ.filtered.Ae.het$Ae_alt.depth==0 
vcf.Ae.Ol.no.bolting.unique.diploid.Ae.het <- vcf.Ae.Ol.no.bolting.unique.diploid[vcf.Ae.Ol.no.bolting.unique.diploid$Ae_gt=="0/1",]

vcf.Ae.Ol.no.bolting.unique.diploid$Ae_alt.depth==0  
``` 

# SNPs between Ae & ref
```{r}
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.GATK.HQ.2.reform.Rdata")

dim(vcf.GATK.HQ.2.reform) # 232294     25 
vcf.tmp <- vcf.GATK.HQ.2.reform 

# subset SNPs between Ae & Ol 
vcf.Ae.Ol.GATK <- vcf.tmp[((vcf.tmp$Ae_gt=="1/1" & vcf.tmp$Ol_gt=="0/0") | (vcf.tmp$Ae_gt=="0/0" & vcf.tmp$Ol_gt=="1/1")),] 

dim(vcf.Ae.Ol.GATK) # 78555    25
table(vcf.Ae.Ol.GATK$Ae_gt)
table(vcf.Ae.Ol.GATK$Ol_gt)
save(vcf.Ae.Ol.GATK, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.GATK.Rdata")
```  

# compare SNPs from freebayes & GATK (all SNPs)
```{r}
# compare all SNPs w/o consideration on genotypes of the two parents 
# load data 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.freebayes.HQ.1.reform.Rdata")
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.GATK.HQ.2.reform.Rdata") 
dim(vcf.freebayes.HQ.1.reform) # 493019     70
dim(vcf.GATK.HQ.2.reform) # 232294     27 
dim(intersection.test) # 160740      3 
table(intersection.test$Ae_gt)
#   0/0   0/1   1/1 
# 47500 43180 70060
table(intersection.test$Ol_gt)
#   0/0   0/1   1/1 
# 50324 40827 69589 

colnames(vcf.freebayes.HQ.1.reform) 
colnames(vcf.GATK.HQ.2.reform)
head(vcf.GATK.HQ.2.reform)

vcf.freebayes.for.compare <- vcf.freebayes.HQ.1.reform[,c("CHROM", "POS", "REF", "ALT", "Ae_gt", "Ol_gt", "Ae_tot.depth", "Ol_tot.depth", "Ae_ref.depth", "Ae_alt.depth", "Ol_ref.depth", "Ol_alt.depth")]
vcf.GATK.for.compare <- vcf.GATK.HQ.2.reform[, c("CHROM", "POS", "REF", "ALT", "Ae_gt", "Ol_gt", "Ae_approx.depth", "Ol_approx.depth", "Ae_genotype.qual", "Ol_genotype.qual", "Ae_ref.depth", "Ae_alt.depth", "Ol_ref.depth", "Ol_alt.depth")]  

head(vcf.freebayes.for.compare)
head(vcf.GATK.for.compare)

vcf.freebayes.for.compare$feature <- paste(vcf.freebayes.for.compare$CHROM, vcf.freebayes.for.compare$POS, vcf.freebayes.for.compare$Ae_gt, vcf.freebayes.for.compare$Ol_gt, sep = ":")

vcf.GATK.for.compare$feature <- paste(vcf.GATK.for.compare$CHROM, vcf.GATK.for.compare$POS, vcf.GATK.for.compare$Ae_gt, vcf.GATK.for.compare$Ol_gt, sep = ":")

intersection.test <- as.data.frame(intersect(vcf.freebayes.for.compare$feature, vcf.GATK.for.compare$feature))
head(intersection.test)
colnames(intersection.test) <- "feature"
intersection.test$Ae_gt <- gsub("([[:print:]]+)(:)([[:digit:]]+)(:)(1/1|0/0|0/1)(:)(1/1|0/0|0/1)", "\\5", intersection.test$feature)
intersection.test$Ol_gt <- gsub("([[:print:]]+)(:)([[:digit:]]+)(:)(1/1|0/0|0/1)(:)(1/1|0/0|0/1)", "\\7", intersection.test$feature)

table(intersection.test$Ae_gt)
table(intersection.test$Ol_gt)
# the intersect also get very high proportion of het, so this can prove that the biological material is not highly homozyougs... 

colnames(vcf.freebayes.for.compare) <- paste("freebayes", c("CHROM", "POS", "REF", "ALT", "Ae_gt", "Ol_gt", "Ae_tot.depth", "Ol_tot.depth", "Ae_ref.depth", "Ae_alt.depth", "Ol_ref.depth", "Ol_alt.depth"),  sep = "|") 
colnames(vcf.GATK.for.compare) <- paste("GATK", c("CHROM", "POS", "REF", "ALT", "Ae_gt", "Ol_gt", "Ae_approx.depth", "Ol_approx.depth", "Ae_genotype.qual", "Ol_genotype.qual", "Ae_ref.depth", "Ae_alt.depth", "Ol_ref.depth", "Ol_alt.depth"), sep = "|")

### find intersection between freebayes & GATK based on CHROM, POS, Ae_gt, and Ol_gt
compare(vcf.freebayes.for.compare, vcf.free)

View(vcf.freebayes.for.compare)
View(vcf.GATK.for.compare) # even not same read count as ref and alt... 
```

# intersection between freebayes & GATK (SNPs between Da-Ae & Da-Ol-1)
```{r}
########### freebayes with Ae & Ol combined 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/freebayes/vcf.Ae.Ol.freebayes.Rdata") 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.GATK.Rdata") 
dim(vcf.Ae.Ol.freebayes) # 137069     68
dim(vcf.Ae.Ol.GATK) # 78555    25

head(vcf.Ae.Ol.freebayes)
head(vcf.Ae.Ol.GATK)

vcf.Ae.Ol.freebayes.tmp <- vcf.Ae.Ol.freebayes
vcf.Ae.Ol.GATK.tmp <- vcf.Ae.Ol.GATK

# goal: get intersections where there SNP calls for both Da-Ae & Da-Ol-1 are the same. 

vcf.freebayes <- paste(vcf.Ae.Ol.freebayes.tmp$CHROM, vcf.Ae.Ol.freebayes.tmp$POS, 
                       vcf.Ae.Ol.freebayes.tmp$REF, vcf.Ae.Ol.freebayes.tmp$ALT,
                       vcf.Ae.Ol.freebayes.tmp$Ae_gt, vcf.Ae.Ol.freebayes.tmp$Ol_gt, 
                       sep = "_")

vcf.GATK <- paste(vcf.Ae.Ol.GATK.tmp$CHROM, vcf.Ae.Ol.GATK.tmp$POS, 
                  vcf.Ae.Ol.GATK.tmp$REF, vcf.Ae.Ol.GATK.tmp$ALT,
                  vcf.Ae.Ol.GATK.tmp$Ae_gt, vcf.Ae.Ol.GATK.tmp$Ol_gt,
                  sep = "_")  

length(intersect(vcf.freebayes, vcf.GATK)) # 62365 same SNP positions between freebayes & GATK 

vcf.Ae.Ol.intersect <- intersect(vcf.freebayes, vcf.GATK)
length(vcf.Ae.Ol.intersect) #  62365 

vcf.Ae.Ol.intersect.df <- data.frame(row.names = 1:length(vcf.Ae.Ol.intersect),
                                     CHROM = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\1", vcf.Ae.Ol.intersect), 
                                     POS = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\3", vcf.Ae.Ol.intersect), 
                                     REF = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\5", vcf.Ae.Ol.intersect), 
                                     ALT = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\7", vcf.Ae.Ol.intersect),
                                     Ae.gt = gsub("([[:print:]]+)(_)(0/0|1/1)(_)((0/0|1/1))", "\\3", vcf.Ae.Ol.intersect),
                                     Ol.gt = gsub("([[:print:]]+)(_)(0/0|1/1)(_)((0/0|1/1))", "\\5", vcf.Ae.Ol.intersect))

head(vcf.Ae.Ol.intersect.df)
ftable(vcf.Ae.Ol.intersect.df$Ae.gt, vcf.Ae.Ol.intersect.df$Ol.gt)  
save(vcf.Ae.Ol.intersect.df, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.intersect.df.Rdata")
write.csv(vcf.Ae.Ol.intersect.df, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.intersect.csv")

########## freebayes with Ae & Ol seperate ######### 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/freebayes/freebayes_result_1.Rdata")
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.GATK.Rdata") 
dim(freebayes_result_1) #  215702      
dim(vcf.Ae.Ol.GATK) # 78555    25

head(freebayes_result_1)
head(vcf.Ae.Ol.GATK)

freebayes_result_1$CHROM <- gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]])", "\\1", freebayes_result_1$feature)
freebayes_result_1$POS <- gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]])", "\\3", freebayes_result_1$feature)
freebayes_result_1$REF <- gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]])", "\\5", freebayes_result_1$feature)
freebayes_result_1$ALT <- gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]])", "\\7", freebayes_result_1$feature)

vcf.Ae.Ol.freebayes.tmp <- freebayes_result_1
vcf.Ae.Ol.GATK.tmp <- vcf.Ae.Ol.GATK

# goal: get intersections where there SNP calls for both Da-Ae & Da-Ol-1 are the same. 

vcf.freebayes <- paste(vcf.Ae.Ol.freebayes.tmp$CHROM, vcf.Ae.Ol.freebayes.tmp$POS, 
                       vcf.Ae.Ol.freebayes.tmp$REF, vcf.Ae.Ol.freebayes.tmp$ALT,
                       vcf.Ae.Ol.freebayes.tmp$Ae_gt, vcf.Ae.Ol.freebayes.tmp$Ol_gt, 
                       sep = "_")

vcf.GATK <- paste(vcf.Ae.Ol.GATK.tmp$CHROM, vcf.Ae.Ol.GATK.tmp$POS, 
                  vcf.Ae.Ol.GATK.tmp$REF, vcf.Ae.Ol.GATK.tmp$ALT,
                  vcf.Ae.Ol.GATK.tmp$Ae_gt, vcf.Ae.Ol.GATK.tmp$Ol_gt,
                  sep = "_")  
head(vcf.freebayes)
head(vcf.GATK)

length(intersect(vcf.freebayes, vcf.GATK)) # 63981 same SNP positions between freebayes & GATK 

vcf.Ae.Ol.intersect.2 <- intersect(vcf.freebayes, vcf.GATK)
length(vcf.Ae.Ol.intersect.2) # 63981
head(vcf.Ae.Ol.intersect.2)

vcf.Ae.Ol.intersect.df.2 <- data.frame(row.names = 1:length(vcf.Ae.Ol.intersect.2),
                                     CHROM = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\1", vcf.Ae.Ol.intersect.2), 
                                     POS = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\3", vcf.Ae.Ol.intersect.2), 
                                     REF = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\5", vcf.Ae.Ol.intersect.2), 
                                     ALT = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\7", vcf.Ae.Ol.intersect.2),
                                     Ae.gt = gsub("([[:print:]]+)(_)(0/0|1/1)(_)((0/0|1/1))", "\\3", vcf.Ae.Ol.intersect.2),
                                     Ol.gt = gsub("([[:print:]]+)(_)(0/0|1/1)(_)((0/0|1/1))", "\\5", vcf.Ae.Ol.intersect.2))

head(vcf.Ae.Ol.intersect.df.2)
head(vcf.Ae.Ol.intersect.df)

ftable(vcf.Ae.Ol.intersect.df.2$Ae.gt, vcf.Ae.Ol.intersect.df.2$Ol.gt)  
save(vcf.Ae.Ol.intersect.df.2, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.intersect.df.2.Rdata")
write.csv(vcf.Ae.Ol.intersect.df.2, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.intersect.2.csv") 

####### make venn diagram between freebayes & GATK
length(vcf.freebayes)  # 215702 
length(vcf.GATK) # 78555    

length(intersect(vcf.GATK, vcf.freebayes)) # 63981  

# #### plot of venn diagram 
library(gplots)
library(VennDiagram)

grid.newpage()
venn.freebayes_GATK <-
draw.pairwise.venn(length(vcf.freebayes), length(vcf.GATK), length(intersect(vcf.GATK, vcf.freebayes)),
                 # category = c("Freebayes:215702", "GATK:78555"),
                 lty = "blank",
                 fill = c("skyblue", "pink1"))
save_plot("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/venn_freebayes_GATK.png", venn.freebayes_GATK,ncol = 1, nrow = 1,base_aspect_ratio = 1, base_height = 4)

#### plot of distribution of final dataset per subgenome per Mb 
########### check SNP distribution by subgenome & chromosome ################
data <- vcf.Ae.Ol.intersect.df.2

# below the block of code should go into function... 
chr_ID <- unique(data$CHROM)
number.unique <- list() 

for (chr in chr_ID) {
  number.unique[chr] <- sum(sum(data$CHROM==chr))
  print(number.unique[chr])
}

SNP_chr.unique <- as.data.frame(t(as.data.frame(number.unique)))
SNP_chr.unique$ID <- rownames(SNP_chr.unique)
head(SNP_chr.unique)
# split to main & random chromosomes 
# main 
SNP_chr_main.unique <- SNP_chr.unique[!grepl("random", SNP_chr.unique$ID),]  
head(SNP_chr_main.unique)
SNP_chr_main.unique$subgenome <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\2", SNP_chr_main.unique$ID)   
SNP_chr_main.unique$chr_ID <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\3", SNP_chr_main.unique$ID) 

# do a little calculation: how many SNPs on main & random chromosomes? 
sum(SNP_chr_main.unique$V1) # 62958     
head(SNP_chr_main.unique) # 200317    

sum(SNP_chr_main.unique[SNP_chr_main.unique$subgenome=="A",]$V1)/sum(SNP_chr_main.unique$V1) 

# make a plot to see how many SNPs on each chromosome 
library(ggplot2)
pl.SNP.main.1.unique <- ggplot(data = SNP_chr_main.unique) 
pl.SNP.main.1.unique <- pl.SNP.main.1.unique + geom_bar(aes(x=factor(chr_ID), y = V1, fill=subgenome), stat = "identity")
pl.SNP.main.1.unique <- pl.SNP.main.1.unique + facet_wrap(~subgenome, ncol = 1)
# pl.SNP.main.1 <- pl.SNP.main.1 + geom_text(aes(x=factor(chr_ID), y = number, label=factor(number), size=1)) 
pl.SNP.main.1.unique <- pl.SNP.main.1.unique + labs(list(title = "", x = "chromosome ID", y = "number of SNPs"))
pl.SNP.main.1.unique <- pl.SNP.main.1.unique + theme(legend.position = "none")
pl.SNP.main.1.unique 
ggsave("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/SNP.main.png", width = 9, height = 7) 

########## Per Mb ... 
head(vcf.Ae.Ol.intersect.df.2)
vcf.Ae.Ol.intersect.df.2 <- vcf.Ae.Ol.intersect.df.2[!grepl("random", vcf.Ae.Ol.intersect.df.2$CHROM),]  
vcf.Ae.Ol.intersect.df.2$subgenome <-  gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\2", vcf.Ae.Ol.intersect.df.2$CHROM) 
vcf.Ae.Ol.intersect.df.2$Chr_ID <- gsub("(chr)(A|C)(01|02|03|04|05|06|07|08|09|10)", "\\3", vcf.Ae.Ol.intersect.df.2$CHROM) 
head(vcf.Ae.Ol.intersect.df.2) 

write.csv(vcf.Ae.Ol.intersect.df.2, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.intersect.df.2.csv") # sort in excel manually 

vcf.Ae.Ol.intersect.df.2.sorted <-  read.csv("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.intersect.df.2.csv")
head(vcf.Ae.Ol.intersect.df.2.sorted) 

library(ggplot2)
pl.SNP.main.2.unique <- ggplot(data = vcf.Ae.Ol.intersect.df.2.sorted)
pl.SNP.main.2.unique <- pl.SNP.main.2.unique + geom_histogram(aes(x=POS, fill=subgenome), binwidth = 1000000) 
pl.SNP.main.2.unique <- pl.SNP.main.2.unique + facet_grid(Chr_ID ~subgenome)
pl.SNP.main.2.unique <- pl.SNP.main.2.unique + labs(list(title = "", x = "chromosome ID", y = "number of SNPs"))
pl.SNP.main.2.unique <- pl.SNP.main.2.unique + theme(legend.position = "none")
pl.SNP.main.2.unique 

ggsave("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/SNP.main.PerMb.png", width = 9, height = 7) 

######
# library(ggplot2)
# pl.SNP.main.2.unique <- ggplot(data = Ae_Ol_after_filter.main.unique)
# pl.SNP.main.2.unique <- pl.SNP.main.2.unique + geom_histogram(aes(x=POS, fill=subgenome), binwidth = 1000000) 
# pl.SNP.main.2.unique <- pl.SNP.main.2.unique + facet_grid(chr_ID ~subgenome)
# pl.SNP.main.2.unique <- pl.SNP.main.2.unique + labs(list(title = "", x = "chromosome ID", y = "number of SNPs"))
# pl.SNP.main.2.unique <- pl.SNP.main.2.unique + theme(legend.position = "none")
# pl.SNP.main.2.unique 
######

######### Ae_ref 
# vcf.Ae.ref.freebayes <- Ae_ref_after_filter.unique[Ae_ref_after_filter.unique$CHROM=="chrA01",]
# nrow(vcf.Ae.ref.freebayes)  # 18059  
# dim(Ae_ref.GATK) # 12502 
# 
# vcf.Ae.ref.freebayes.SNP <- paste(vcf.Ae.ref.freebayes$CHROM, vcf.Ae.ref.freebayes$POS, sep = "_")
# vcf.Ae.ref.GATK.SNP <- paste(Ae_ref.GATK$CHROM, Ae_ref.GATK$POS, sep = "_") 
# 
# intersect(vcf.Ae.ref.GATK.SNP, vcf.Ae.ref.freebayes.SNP) 
# 
# #### plot of venn diagram 
# 
# grid.newpage() 
# venn.freebayes_GATK.Ae_ref <- 
# draw.pairwise.venn(nrow(vcf.Ae.ref.freebayes), nrow(Ae_ref.GATK), length(intersect(vcf.Ae.ref.GATK.SNP, vcf.Ae.ref.freebayes.SNP)), 
#                  category = c("freebayes (11379)", "GATK (8491)"), 
#                  lty = "blank", 
#                  fill = c("skyblue", "pink1"))
# 
# ######### Ol_ref 
# vcf.Ol.ref.freebayes <- Ol_ref_after_filter.unique[Ol_ref_after_filter.unique$CHROM=="chrA01",]
# nrow(vcf.Ol.ref.freebayes)  # 16225   
# dim(Ol_ref.GATK) # 11352  
# 
# vcf.Ol.ref.freebayes.SNP <- paste(vcf.Ol.ref.freebayes$CHROM, vcf.Ol.ref.freebayes$POS, sep = "_")
# vcf.Ol.ref.GATK.SNP <- paste(Ol_ref.GATK$CHROM, Ol_ref.GATK$POS, sep = "_") 
# 
# intersect(vcf.Ol.ref.GATK.SNP, vcf.Ol.ref.freebayes.SNP) 
# 
# #### plot of venn diagram 
# 
# grid.newpage() 
# venn.freebayes_GATK.Ol_ref <- 
# draw.pairwise.venn(nrow(vcf.Ol.ref.freebayes), nrow(Ol_ref.GATK), length(intersect(vcf.Ol.ref.GATK.SNP, vcf.Ol.ref.freebayes.SNP)), 
#                  category = c("freebayes (11379)", "GATK (8491)"), 
#                  lty = "blank", 
#                  fill = c("skyblue", "pink1")) 


######## compare the two intersect dataset... use the intersect...  
vcf.Ae.Ol.intersect.final <- intersect(vcf.Ae.Ol.intersect.2, vcf.Ae.Ol.intersect)
length(vcf.Ae.Ol.intersect.final) # 59343 # common between combined & sepereate freebayes call 
head(vcf.Ae.Ol.intersect.final) 

test.combined <- vcf.Ae.Ol.intersect.df.2[,c("CHROM", "POS")]
test.sep <- vcf.Ae.Ol.intersect.df[,c("CHROM", "POS")]

tmp2 <- anti_join(test.combined, test.sep)
dim(tmp2)

tmp3 <- anti_join(test.sep, test.combined)
dim(tmp3) 

vcf.Ae.Ol.intersect.df.final <- data.frame(row.names = 1:length(vcf.Ae.Ol.intersect.final),
                                     CHROM = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\1", vcf.Ae.Ol.intersect.final), 
                                     POS = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\3", vcf.Ae.Ol.intersect.final), 
                                     REF = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\5", vcf.Ae.Ol.intersect.final), 
                                     ALT = gsub("([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)(_)([[:print:]]+)", "\\7", vcf.Ae.Ol.intersect.final),
                                     Ae.gt = gsub("([[:print:]]+)(_)(0/0|1/1)(_)((0/0|1/1))", "\\3", vcf.Ae.Ol.intersect.final),
                                     Ol.gt = gsub("([[:print:]]+)(_)(0/0|1/1)(_)((0/0|1/1))", "\\5", vcf.Ae.Ol.intersect.final))

dim(vcf.Ae.Ol.intersect.df.final) # 59343     6 

head(vcf.Ae.Ol.intersect.df.final)
ftable(vcf.Ae.Ol.intersect.df.final$Ae.gt, vcf.Ae.Ol.intersect.df.final$Ol.gt)  
save(vcf.Ae.Ol.intersect.df.final, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.intersect.df.final.Rdata")
write.csv(vcf.Ae.Ol.intersect.df.final, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.intersect.final.csv") 

``` 

# compare mapping result to napus, rapa, and juncea for Da-Ae & Da-Ol-1 
```{r}
# 1) get material, flowering stage sample for Da-Ae & Da-Ol-1 
# done

# 2) index reference rapa & juncea genome 
# rapa genome 
# juncea genome 

# 3) map 

# 4) check & compare result  
#### total mapping summary  
total_mapping <- read.csv("~/Desktop/Brassica_project/KIAT_RNA_seq/mapping_napus_rapa_juncea/mapping_summary.csv", header = F)
head(total_mapping)
total_mapping$value <- as.numeric(gsub("%", "", total_mapping$V3, fixed=T))/100
total_mapping$ref <- gsub("([[:print:]]+)(rapa|napus|juncea)", "\\2", total_mapping$V1)
total_mapping$lib <- gsub("([[:print:]]+)(_)(paired.star.dir.rapa|paired.star.dir.napus|paired.star.dir.juncea)", "\\1", total_mapping$V1)
total_mapping$gt <- c(rep("Da-Ae",24), rep("Da-Ol-1", 24))
total_mapping$lib

p.total_mapping <- ggplot(data = total_mapping) + theme_grey()
p.total_mapping <- p.total_mapping + geom_bar(aes(x=factor(V2), y = value, fill=V2), stat="identity")
p.total_mapping <- p.total_mapping + facet_grid(ref~lib)
p.total_mapping <- p.total_mapping + geom_text(data = total_mapping, aes(x=factor(V2), y=value, label=factor(value)))
# geom_text(data=DEGs_number_between_gt.melt.parent,aes(x=DE,y=number*1.05,label=factor(number)))
p.total_mapping 
ggsave(p.total_mapping, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/mapping_napus_rapa_juncea/p.total_mapping.png", width = 12, height=8)

# conlussions: no difference between Da-Ae and Da-Ol-1's mapping ratio to different ref genomes. napus has the highest mapping ratio among the 3 ref genomes, followed by juncea and rapa.  

#### subgenome mapping summary (only unique mapped reads)
subgenome_mapping <- read.csv("~/Desktop/Brassica_project/KIAT_RNA_seq/mapping_napus_rapa_juncea/subgenome_mapping_summary.csv", header = F)
head(subgenome_mapping, 20)
subgenome_mapping <- subgenome_mapping[grep("^A$|^B$|^C$", subgenome_mapping$V2),]
unique_mapping_sum <- aggregate(subgenome_mapping$V1, list(subgenome_mapping$V3), sum)
names(unique_mapping_sum) <- c("V3", "total")
head(unique_mapping_sum)
subgenome_mapping.2 <- merge(subgenome_mapping, unique_mapping_sum, by="V3")
subgenome_mapping.2$ratio <- subgenome_mapping.2$V1/subgenome_mapping.2$total

subgenome_mapping.2$ref <- gsub("([[:print:]]+)(_paired.star.dir.)(napus|rapa|juncea)([[:punct:]]+)(subgenome_mapping_summary_main)", "\\3", subgenome_mapping.2$V3)
subgenome_mapping.2$lib <- gsub("([[:print:]]+)(_paired.star.dir.)(napus|rapa|juncea)([[:punct:]]+)(subgenome_mapping_summary_main)", "\\1", subgenome_mapping.2$V3)
subgenome_mapping.2$gt <- c(rep("Da-Ae", 10), rep("Da-Ol-1", 10))

subgenome_mapping.2$lib
subgenome_mapping.2$ratio
subgenome_mapping.2$V2 

p.subgenome_mapping <- ggplot(data = subgenome_mapping.2)
p.subgenome_mapping <- p.subgenome_mapping + geom_bar(aes(x=V2, y = ratio, fill=V2), stat="identity")
p.subgenome_mapping <- p.subgenome_mapping + facet_grid(ref~lib)
p.subgenome_mapping <- p.subgenome_mapping + geom_text(data = subgenome_mapping.2, aes(x=factor(V2), y=ratio, label=factor(round(ratio, digits = 2))))
# geom_text(data=DEGs_number_between_gt.melt.parent,aes(x=DE,y=number*1.05,label=factor(number)))
p.subgenome_mapping 
ggsave(p.subgenome_mapping, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/mapping_napus_rapa_juncea/p.subgenome_mapping.png", width = 12, height=8)
# conlussion: the lowest but good partion of reads mapped to B subgenome, and the percentage is about the same between Da-Ae & Da_Ol-1. This suggest that that the genetic makeup between Da-Ae & Da-Ol-1 is not that different... 

# Ae_mapped <- read.delim("~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/data/Ae_mapping_chr_2", header = F)
# Ae_mapped$value <- gsub("([[:digit:]]+)(|)([[:print:]]+)", "\\1", Ae_mapped$V1)
# Ae_mapped
# 
# 
# Ol_mapped <- read.delim("~/Desktop/Brassica_project/KIAT_RNA_seq/SNP_calling/data/Ol_mapping_chr_2", header = F)
# Ol_mapped$value <- gsub("([[:digit:]]+)(|)([[:print:]]+)", "\\1", Ol_mapped$V1)
# Ol_mapped
# 
# sum(as.numeric(Ae_mapped[c(1,3,5,7,9,11,13,15,17,19), 2])) 
# sum(as.numeric(Ol_mapped[c(1,3,5,7,9,11,13,15,17,19), 2]))
# sum(as.numeric(Ae_mapped[c(1,3,5,7,9,11,13,15,17,19), 2])) + sum(as.numeric(Ol_mapped[c(1,3,5,7,9,11,13,15,17,19), 2]))
# # 383980934 A subgenome mapped reads 
# sum(as.numeric(Ae_mapped[c(22,24,26,28,30,32,34,36,38), 2])) + sum(as.numeric(Ol_mapped[c(22,24,26,28,30,32,34,36,38), 2])) 
# sum(as.numeric(Ae_mapped[c(22,24,26,28,30,32,34,36,38), 2]))
# sum(as.numeric(Ol_mapped[c(22,24,26,28,30,32,34,36,38), 2])) 
# # 786750509 C subgenome mapped reads 
# 
# Ae_mapped[c(1,3,5,7,9,11,13,15,17,19), 1]
# Ol_mapped[c(1,3,5,7,9,11,13,15,17,19), 1]
# Ae_mapped[c(22,24,26,28,30,32,34,36,38), 1]
# Ol_mapped[c(22,24,26,28,30,32,34,36,38), 1]

#### reads level mapping detail check ... 
reads_mapping <- read.csv("~/Desktop/Brassica_project/KIAT_RNA_seq/mapping_napus_rapa_juncea/reads_mapping_summary.csv", header = F)
head(reads_mapping)
reads_mapping <- reads_mapping[grep("^A$|^C$|^non_unique$|^unmapped$", reads_mapping$V2),]
reads_mapping_sum <- aggregate(reads_mapping$V1, list(reads_mapping$V3), sum)
names(reads_mapping_sum) <- c("V3", "total")
head(reads_mapping_sum)
reads_mapping.2 <- merge(reads_mapping, reads_mapping_sum, by="V3")
reads_mapping.2$ratio <- reads_mapping.2$V1/reads_mapping.2$total

reads_mapping.2

reads_mapping.2$ref <- gsub("([[:print:]]+)(juncea|rapa)(_)(6_Da-Ae_Asub.summary|6_Da-Ae_Csub.summary|All1_Gae_Da-Ol_Asub.summary|All1_Gae_Da-Ol_Csub.summary|6_Ae_summary|All1_Gae_Ol_summary)", "\\2", reads_mapping.2$V3)
reads_mapping.2$sub <- c(rep("A", 4), rep("B", 4), rep("A", 3), rep("B", 3), rep("A", 7))
reads_mapping.2$gt <- c(rep("Da-Ae", 8), rep("Da-Ol-1", 6), rep("Da-Ae", 4), rep("Da-Ol-1", 3))
reads_mapping.2$group <- paste(reads_mapping.2$ref, reads_mapping.2$sub, sep = "_")

p.reads_mapping <- ggplot(data = reads_mapping.2)
p.reads_mapping <- p.reads_mapping + geom_bar(aes(x=V2, y = ratio, fill=V2), stat="identity")
p.reads_mapping <- p.reads_mapping + facet_grid(group~gt)
p.reads_mapping <- p.reads_mapping + geom_text(data = reads_mapping.2, aes(x=factor(V2), y=ratio, label=factor(round(ratio, digits = 2))))
# geom_text(data=DEGs_number_between_gt.melt.parent,aes(x=DE,y=number*1.05,label=factor(number)))
p.reads_mapping 
ggsave(p.reads_mapping, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/mapping_napus_rapa_juncea/p.reads_mapping.png", width = 12, height=8)

# conclussion: 

# p.reads_mapping.2 <- ggplot(data = reads_mapping.2) 
# p.reads_mapping.2 <- p.reads_mapping.2 + geom_bar(aes(x="", y= ratio, fill=V2), stat="identity")
# p.reads_mapping.2 <- p.reads_mapping.2 + coord_polar("y") 
# p.reads_mapping.2 <- p.reads_mapping.2 + facet_grid(gt~group) + theme(axis.text.x=element_blank())
# p.reads_mapping.2 <- p.reads_mapping.2 + geom_text(data = reads_mapping.2, aes(x="", y=ratio, label=factor(round(ratio, digits = 2)), size=3), check_overlap = TRUE) 
# p.reads_mapping.2 
# ggsave(p.reads_mapping, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/mapping_napus_rapa_juncea/p.reads_mapping.png", width = 12, height=8)

```

# SNP filter inside linux or R 
```{r}
######## SNP filtering can be done inside R, but need bgzip & tabix to index the vcf file.... 
#### spent quite a long time to figure this out... 

### visulize unique & non-unique SNPs in IGV, done... 

### SNP filtering inside linux using script (see /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/2016_fall/SNP_calling/alignment_freebayes/vcf_file/log file in whitney)

### alignment score can also be taken into consideration for filtering before SNP calling, because there are multiple possible causes of false-positive SNPs: 1) sequencing error [Phred score] 2) alignment error [mapping score inside bam file] 3) incompleness of ref genome (ref: Wang et al. 2015; )  

### filter based on read depth (DP)? alternate allele observation count (AO)? QA (sum of quality of the alternate observation count)  

Ae_Ol_after_filter.main.unique.test <- Ae_Ol_after_filter.main.unique[Ae_Ol_after_filter.main.unique$CHROM=="chrA02",]

colnames(Ae_Ol_after_filter.main.unique.test)
total_depth_check <- Ae_Ol_after_filter.main.unique.test[,c("CHROM","POS", "Ae_tot.depth", "Ol_tot.depth")]
head(total_depth_check)
total_depth_check_long <- melt(total_depth_check, id.vars = c("CHROM", "POS"))
head(total_depth_check_long)

pl.total_depth <- ggplot(data = total_depth_check_long)
pl.total_depth <- pl.total_depth + geom_line(aes(x=POS, y = value, color=variable))
# pl.total_depth <- pl.total_depth + facet_facet(~variable)
# pl.SNP.main.2 <- pl.SNP.main.2 + labs(list(title = "", x = "chromosome ID", y = "number of SNPs"))
# pl.SNP.main.2 <- pl.SNP.main.2 + theme(legend.position = "none")
pl.total_depth

# based on the total depth visualization result, some regions are having high depth for both Ae and Ol, whereas some regions are only exibiting high depth for one of the genotype 

hist(Ae_Ol_after_filter.main.unique.test[Ae_Ol_after_filter.main.unique.test$Ae_tot.depth>10000,]$Ae_tot.depth)
# average depth 
# Ae 
sum(Ae_Ol_after_filter.main.unique$Ae_tot.depth)/nrow(Ae_Ol_after_filter.main.unique) # 185.4755 
# Ol 
sum(Ae_Ol_after_filter.main.unique$Ol_tot.depth)/nrow(Ae_Ol_after_filter.main.unique) # 207.3724 
# 

```

# SNP annotation 
```{r}
# snpEff # this version of annotated vcf file is not correct... because "ERROR_CHROMOSOME_NOT_FOUND" 
########### 1) 
# build my own database ############### 
# 1) download snpEff databse, unzip 
# mkdir data 
# cd data 
# mkdir napus 
# cd napus  
# use the gff3 file in 
# /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/Reference/B.napus/Brassica_napus.annotation_v5_modified_modified.gff3 
#########
# mv Brassica_napus.annotation_v5_modified_modified.gff3 genes.gff 
# wget ftp://bradata:zhl410ivf@brassicadb.org/Brassica_napus/Brassica_napus_v4.1.chromosomes.fa.gz
# mv Brassica_napus_v4.1.chromosomes.fa.gz sequences.fa.gz 
# gunzip *.gz
# cd ..
# mkdir genomes  
# cd genomes/ 
# wget ftp://bradata:zhl410ivf@brassicadb.org/Brassica_napus/Brassica_napus_v4.1.chromosomes.fa.gz
# mv Brassica_napus_v4.1.chromosomes.fa.gz sequences.fa.gz 
# gunzip sequences.fa.gz
# cd ..  (to snpEff home directory)
# vim snpEff.config #### 
# add "# Brassica napus genome  napus.genome: Brassica_napus" 
# java -jar snpEff.jar build -gff3 -v napus 

### SNP annotation  
# cd cd ~/bin/snpEff/
# (GATK version) java -jar snpEff.jar napus /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/result/GATK/raw_vcf/vcf_split_by_chr/GATK_SNP_biallelic.recode.vcf > /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/result/GATK/raw_vcf/vcf_split_by_chr/GATK_SNP_biallelic.recode.ann.vcf 

# (freebayes version) java -jar snpEff.jar napus /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/result/GATK/raw_vcf/vcf_split_by_chr/GATK_SNP_biallelic.recode.vcf > /Network/Servers/avalanche.plb.ucdavis.edu/Volumes/Mammoth/Users/ruijuanli/SNP_parent/result/GATK/raw_vcf/vcf_split_by_chr/GATK_SNP_biallelic.recode.ann.vcf 

``` 

# randomly pick SNPs that can be confirmed in F2 VS. not found in F2 
```{r}
library(Biostrings)
source("/Users/ruijuanli/Desktop/Brassica_project/KIAT_RNA_seq/analysis/function_BnRNAseq.R")

# get genome fasta sequence
napus <-readDNAStringSet("~/Desktop/Brassica_project/IGV_file/Brassica_napus_v4.1.chromosomes.fa")

# get GATK final vcf content 
load("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.Ae.Ol.GATK.Rdata") 

######### get SNP present in F2 #############
F2.Y <- read.csv("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/final/Final_F2_SNP_Sites.csv")  

test.2.F2.Y <- SNP.reform(F2.Y) # randomly pick 150 because otherwise designing primer will become a problem, since we want to avoid picking primer from the 20-30 bp flanking region of the candiate SNPs (sanger sequencing has bad quality for the 1st 20-30 nucleotides).  

seq.final.F2.Y <- get.fasta(test.2.F2.Y, napus) # get fasta seq for the 150 SNPs
vcf.F2.Y <- get.vcf(test.2.F2.Y, vcf.Ae.Ol.GATK) 
sum((vcf.F2.Y$CHROM %in% test.2.F2.Y$CHROM) & (vcf.F2.Y$POS %in% test.2.F2.Y$POS)) # check, correct!!! 
test.2.F2.Y

save(vcf.F2.Y, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.F2.Y")
writeXStringSet(seq.final.F2.Y, filepath = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/SNP_validation/seq.final.F2.fa") 

# get 2nd dataset because the number of SNPs are not enough 
test.2.F2.Y.2 <- SNP.reform(F2.Y) # randomly pick 150 because otherwise designing primer will become a problem, since we want to avoid picking primer from the 20-30 bp flanking region of the candiate SNPs (sanger sequencing has bad quality for the 1st 20-30 nucleotides).  

seq.final.F2.Y.2 <- get.fasta(test.2.F2.Y, napus) # get fasta seq for the 150 SNPs
vcf.F2.Y.2 <- get.vcf(test.2.F2.Y.2, vcf.Ae.Ol.GATK) 
sum((vcf.F2.Y.2$CHROM %in% test.2.F2.Y.2$CHROM) & (vcf.F2.Y.2$POS %in% test.2.F2.Y.2$POS)) # check, correct!!! 
test.2.F2.Y.2

save(vcf.F2.Y.2, file = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/vcf.F2.2.Y")
writeXStringSet(seq.final.F2.Y.2, filepath = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/SNP_validation/seq.final.F2.2.fa") 

######### get SNP not present in F2 ############
# get SNP data
all <- read.csv("~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/data/final/vcf.Ae.Ol.intersect.final.csv",row.names = 1) 

all$feature <- paste(all$CHROM, all$POS, sep = "_") 
F2.Y$feature <- paste(F2.Y$CHROM, F2.Y$POS, sep = "_")

sum(F2.Y$feature %in% all$feature)
sum(all$feature %in% F2.Y$feature)

F2.N <- all[(all$feature %in% F2.Y$feature) == "FALSE",]
dim(F2.N)
dim(all) 
dim(F2.Y)

dim(F2.Y) + dim(F2.N)

test.2.F2.N <- SNP.reform(F2.N) # randomly pick 96 
seq.final.F2.N <- get.fasta(test.2.F2.N, napus) # get fasta seq for the 96 SNPs
vcf.F2.N <- get.vcf(test.2.F2.N, vcf.Ae.Ol.GATK) # get vcf 
sum((vcf.F2.N$CHROM %in% test.2.F2.N$CHROM) & (vcf.F2.N$POS %in% test.2.F2.N$POS)) # check, correct!!! 

####### plot depth & quality for F2.Y & F2.N dataset  

########## present in F2 with good quality ##############
depth.F2.Y <- vcf.F2.Y[,c("Ae_approx.depth", "Ol_approx.depth")]
GQ.F2.Y <- vcf.F2.Y[,c("Ae_genotype.qual", "Ol_genotype.qual")]

# reformat 
library(reshape)
depth.F2.Y.long <- melt(depth.F2.Y)
depth.F2.Y.long$gt <- gsub("(Ae|Ol)(_)(approx.depth)","\\1",depth.F2.Y.long$variable)
depth.F2.Y.long$gt

GQ.F2.Y.long <- melt(GQ.F2.Y)
GQ.F2.Y.long$gt <- gsub("(Ae|Ol)(_)(genotype.qual)","\\1",GQ.F2.Y.long$variable) 
GQ.F2.Y.long$gt

# plot total depth by Ae_gt  
p.tot.depth.F2.Y <- ggplot(data = depth.F2.Y.long) 
p.tot.depth.F2.Y <- p.tot.depth.F2.Y + geom_histogram(aes(value, fill=gt)) # + scale_x_log10() 
p.tot.depth.F2.Y <- p.tot.depth.F2.Y + facet_wrap(~gt, ncol = 1)
p.tot.depth.F2.Y  
ggsave(p.tot.depth.F2.Y, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/p.tot.depth.F2.Y.png", width = 6, height = 8)

# plot GQ
p.GQ.F2.Y <- ggplot(data = GQ.F2.Y.long) 
p.GQ.F2.Y <- p.GQ.F2.Y + geom_histogram(aes(value, fill=gt))  
p.GQ.F2.Y <- p.GQ.F2.Y + facet_wrap(~gt, ncol = 1)
p.GQ.F2.Y 
ggsave(p.GQ.F2.Y, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/p.GQ.F2.Y.png", width = 6, height = 8) 

########### not present in F2 
depth.F2.N <- vcf.F2.N[,c("Ae_approx.depth", "Ol_approx.depth")]
GQ.F2.N <- vcf.F2.N[,c("Ae_genotype.qual", "Ol_genotype.qual")]

# reformat 
depth.F2.N.long <- melt(depth.F2.N)
depth.F2.N.long$gt <- gsub("(Ae|Ol)(_)(approx.depth)","\\1",depth.F2.N.long$variable)
depth.F2.N.long$gt

GQ.F2.N.long <- melt(GQ.F2.N)
GQ.F2.N.long$gt <- gsub("(Ae|Ol)(_)(genotype.qual)","\\1",GQ.F2.N.long$variable) 
GQ.F2.N.long$gt

# plot total depth by Ae_gt 
p.tot.depth.F2.N <- ggplot(data = depth.F2.N.long) 
p.tot.depth.F2.N <- p.tot.depth.F2.N + geom_histogram(aes(value, fill=gt)) # + scale_x_log10() 
p.tot.depth.F2.N <- p.tot.depth.F2.N + facet_wrap(~gt, ncol = 1)
p.tot.depth.F2.N  
ggsave(p.tot.depth.F2.N, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/p.tot.depth.F2.N.png", width = 6, height = 8)

# plot GQ
p.GQ.F2.N <- ggplot(data = GQ.F2.N.long) 
p.GQ.F2.N <- p.GQ.F2.N + geom_histogram(aes(value, fill=gt))  
p.GQ.F2.N <- p.GQ.F2.N + facet_wrap(~gt, ncol = 1)
p.GQ.F2.N 
ggsave(p.GQ.F2.N, filename = "~/Desktop/Brassica_project/KIAT_RNA_seq/parent_SNP/output/figure/p.GQ.F2.N.png", width = 6, height = 8)  

# design primer using primer-BLAST 
# use primer blast to design primer one by one: 
# rules to follow: 
# 1) SNPs must fall within the targeted region
# 2) avoid having SNP in the 1st 20nt of sanger sequencing, that's forward primer end position + 20 && reverse primer start postion - 20 
# 3) if the tempelate is not unique, don't use that SNP 

### To do: grep 300 SNPs, and include only SNPs in the seq name... 




```