GATK_VCF_version2.Rmd

---
title: "gatk version2"
author: "Shankar K Shakya"
date: "May 22, 2017"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.width = 18, fig.height = 12, echo = F, warning = F, message = F, verbose = F)
```

### GATK generated variant analysis

```{r, echo=TRUE}

rm(list = ls())
library(vcfR)
library(poppr)
library(ggplot2)
library(reshape2)
vcf.gatk <- read.vcfR("Pcinna144.vcf")

# 144 samples, 178110 variants

```

## Filtering samples by depth

Varaints that were sequenced at unusual depth, i.e less than 10 % quantile and greater than 80% quantile is converted to NA. Minimum depth of 4x is maintained per variant.

```{r, echo=TRUE}
dp <- extract.gt(vcf.gatk, element = "DP", as.numeric=TRUE)

quants <- apply(dp, MARGIN=2, quantile, probs=c(0.1, 0.8), na.rm=TRUE)
dp2 <- sweep(dp, MARGIN=2, FUN = "-", quants[1,])
dp[dp2 < 0] <- NA

dp2 <- sweep(dp, MARGIN=2, FUN = "-", quants[2,])
dp[dp2 > 0] <- NA

dp[dp < 4] <- NA

vcf.gatk@gt[,-1][ is.na(dp) == TRUE ] <- NA

```


## Omitting variants


```{r, echo=TRUE}

dp <- extract.gt(vcf.gatk, element = "DP", as.numeric=TRUE)

myMiss <- apply(dp, MARGIN = 1, function(x){ sum( is.na(x) ) } )
myMiss <- myMiss / ncol(dp)
vcf.gatk <- vcf.gatk[myMiss < 0.1, ]


```


## Omitting samples

```{r, echo=TRUE}

myMiss <- apply(dp, MARGIN = 2, function(x){ sum( is.na(x))} )
myMiss <- myMiss / nrow(dp)
vcf.gatk@gt <- vcf.gatk@gt[, c(TRUE, myMiss < 0.7)]

```



## Filtering by Minor allele frequency

```{r, echo=TRUE}

mymaf <- maf(vcf.gatk, element = 2)
mymaf <- mymaf[mymaf[,4] > 0.05, ]

vcf.gatk@fix[,3] <- paste(vcf.gatk@fix[,1], vcf.gatk@fix[,2], sep = "_")

fix <- vcf.gatk@fix[vcf.gatk@fix[,3] %in% rownames(mymaf), ]
true_ind <- which(vcf.gatk@fix[,3] %in% rownames(mymaf))
vcf.gatk@fix <- vcf.gatk@fix[true_ind, ]

vcf.gatk@gt <- vcf.gatk@gt[true_ind, ]
# 
#write.vcf(vcf.gatk, file = "MAF_Filtered.vcf.gatk6262.gz")

```

## Depth Plot

```{r}
library(vcfR)
library(ggplot2)

#vcf.gatk <- read.vcfR("MAF_Filtered.vcf.gatk6262.gz")

dp <- extract.gt(vcf.gatk, element = "DP", as.numeric=TRUE)
par(mar=c(12,4,4,2))
boxplot(dp, col=2:8, las=3)
title(ylab = "Depth (DP)")

```


## Converting VCFR object to genlight object


```{r}

library(adegenet)

vcf.gatk.gl <- vcfR2genlight(vcf.gatk)
mypca <- glPca(vcf.gatk.gl, parallel = FALSE, nf = 20)
#saveRDS(mypca, file = "vcf.gatk.nf20_pca.RData")

#mypca <- readRDS("vcf.gatk.nf20_pca.RData")

pcascores <- mypca$scores
rownames(pcascores) <- unlist(strsplit(rownames(pcascores), split = ".fq"))

# rownames(pcascores)[grep("port", rownames(pcascores))] <- "Portugal"
# rownames(pcascores)[grep("ntaiw", rownames(pcascores))] <- "N.Taiwan"
# rownames(pcascores)[grep("staiw", rownames(pcascores))] <- "S.Taiwan"
# rownames(pcascores)[grep("centtaiw", rownames(pcascores))] <- "Cent.Taiwan"
# rownames(pcascores)[grep("nsw", rownames(pcascores))] <- "NSW"
# rownames(pcascores)[grep("safrica", rownames(pcascores))] <- "SAfrica"
# rownames(pcascores)[grep("nvietnam", rownames(pcascores))] <- "N.Vietnam"
# rownames(pcascores)[grep("waus", rownames(pcascores))] <- "W.Aus"
# rownames(pcascores)[grep("france", rownames(pcascores))] <- "France"
# rownames(pcascores)[grep("italy", rownames(pcascores))] <- "Italy"
# rownames(pcascores)[grep("domrep", rownames(pcascores))] <- "DomRep"
# rownames(pcascores)[grep("chile", rownames(pcascores))] <- "Chile"
# rownames(pcascores)[grep("newguin", rownames(pcascores))] <- "PNG"
# rownames(pcascores)[grep("tasaus", rownames(pcascores))] <- "Tasmania"
# rownames(pcascores)[grep("queensaus", rownames(pcascores))] <- "Queensland"
# rownames(pcascores)[grep("netaiw", rownames(pcascores))] <- "NE.Taiwan"

# rownames(pcascores)[grep("Port", rownames(pcascores))] <- "Portugal"
# rownames(pcascores)[grep("Taiw", rownames(pcascores))] <- "Taiwan"
# rownames(pcascores)[grep("SAfrica", rownames(pcascores))] <- "SAfrica"
# rownames(pcascores)[grep("France", rownames(pcascores))] <- "France"
# rownames(pcascores)[grep("Italy", rownames(pcascores))] <- "Italy"
# rownames(pcascores)[grep("DomRep", rownames(pcascores))] <- "DomRep"
# rownames(pcascores)[grep("Chile", rownames(pcascores))] <- "Chile"
# rownames(pcascores)[grep("NewGuin", rownames(pcascores))] <- "PNG"
# rownames(pcascores)[grep(c("QueensAus|TasAus|WAus|NSWalesAus"), rownames(pcascores))] <- "Australia"
# rownames(pcascores)[grep("NVietnam", rownames(pcascores))] <- "Vietnam"
# 
# 
# df <- mypca$scores
# df <- data.frame(Country = c(rownames(pcascores)), df)
# 
# library(ggfortify)
# 
# rownames(pcascores)[grep("Portugal|France|Italy", rownames(pcascores))] <- "Europe"
# rownames(pcascores)[grep("Taiwan|Vietnam", rownames(pcascores))] <- "Asia"
# rownames(pcascores)[grep("SAfrica", rownames(pcascores))] <- "Africa"
# rownames(pcascores)[grep("Chile" , rownames(pcascores))] <- "South America"
# rownames(pcascores)[grep("Australia" , rownames(pcascores))] <- "Australia"
# rownames(pcascores)[grep("PNG" , rownames(pcascores))] <- "PNG"
# rownames(pcascores)[grep("DomRep" , rownames(pcascores))] <- "DomRep"
# 
# df <- data.frame(Continent = c(rownames(pcascores)), df)
# 
# newdf <- df[,c(3,4)]
# 
# 
# head(pcascores)
# rownames(pcascores)

pop <- read.csv("Pcinna_pop.csv", header = TRUE)
pop <- pop[pop$Isolate %in% rownames(pcascores), ]
pop <- pop[match(rownames(pcascores), pop$Isolate), ]
#newpca_pop <- cbind(pop, pcascores)


```


## Discriminant analysis of principal components


```{r, echo=TRUE}

library(poppr)

pop(vcf.gatk.gl) <- pop$Country
mydapc <- dapc(vcf.gatk.gl, pop = pop(vcf.gatk.gl), parallel = F,  n.pca = 50, n.da = 5, glPca = mypca )

scatter.dapc(mydapc, clabel = 0.75, pch=15:20, scree.pca = TRUE, scree.da = FALSE, 
        posi.pca = "topright", posi.leg = "topleft", legend = F, 
        cleg = 0.9, inset.solid = 1, xax = 1, yax = 2, cex.lab = 1, cex = 2, solid = 1, cstar = 0)



newpop <- popsub(vcf.gatk.gl , blacklist = c("Taiwan", "Vietnam"))

newdapc <- dapc(newpop, pop = pop(newpop), n.pca = 50, n.da = 5, parallel = FALSE)

scatter.dapc(newdapc, clabel = 0.75, pch=15:19, scree.pca = TRUE, scree.da = FALSE, 
        posi.pca = "topright", posi.leg = "topleft", legend = F, 
        cleg = 0.9, inset.solid = 1, xax = 2, yax = 1, cex.lab = 1, cex = 1.5, solid = 1, cstar = 0)



newpop2 <- popsub(vcf.gatk.gl , blacklist = c("Chile", "Taiwan", "Vietnam"))

newdapc2 <- dapc(newpop2, pop = pop(newpop2), n.pca = 25, n.da = 5, parallel = FALSE)

scatter.dapc(newdapc2, clabel = 0.75, pch=15:19, scree.pca = TRUE, scree.da = FALSE, 
        posi.pca = "topright", posi.leg = "topleft", legend = F, 
        cleg = 0.9, inset.solid = 1, xax = 2, yax = 1, cex.lab = 1, cex = 1.5, solid = 1, cstar = 0)




newpop3 <- popsub(vcf.gatk.gl , blacklist = c("PNG", "Chile", "Taiwan", "Vietnam"))

newdapc3 <- dapc(newpop3, pop = pop(newpop3), n.pca = 20, n.da = 4, parallel = FALSE)

scatter.dapc(newdapc3, clabel = 0.75, pch=15:19, scree.pca = TRUE, scree.da = FALSE, 
        posi.pca = "topright", posi.leg = "topleft", legend = F, 
        cleg = 0.9, inset.solid = 1, xax = 2, yax = 1, cex.lab = 1, cex = 1.5, solid = 1, cstar = 0)



```



## Index of association

```{r, echo=TRUE}

library(poppr)

pop(vcf.gatk.gl) <- pop$Country
ia_pop <- seppop(vcf.gatk.gl ) %>% lapply(samp.ia) 

par(mfrow = c(2,5))
lapply(names(ia_pop), function(x) hist(ia_pop[[x]], main=x))

mean <- lapply(ia_pop, function(x) mean(na.omit(x)))

```

## Index of association boxplot and barplot using variants with no missing data

```{r}

vcf_ia <- vcf.gatk
fix <- vcf_ia@fix
gt <- vcf_ia@gt
gt_vcf.ia <- extract.gt(vcf_ia)

vcf_ia@fix <- fix[rowSums(!is.na(gt_vcf.ia)) == 127, ]
vcf_ia@gt <- gt[fix[,3] %in% vcf_ia@fix[,3], ]

vcf_ia_gl <- vcfR2genlight(vcf_ia)
pop(vcf_ia_gl) <- pop$Country


rf.ia <- seppop(vcf_ia_gl) %>% # separate each population
  c(Total = vcf_ia_gl) %>%     # add the total population
  lapply(samp.ia, threads = 0, n.snp = 100L, reps = 100L) %>%
  data.frame %>%    # convert list to data frame w/ 1000 rows
  tidyr::gather(POP, value)

rf.ia$POP <- factor(rf.ia$POP)

ggia <- ggplot(rf.ia, aes(x = POP, y = value)) +
  geom_boxplot() +
  theme_bw() +
  theme(panel.grid.major.x = element_blank()) +
  theme(text = element_text(size = 18)) +
  ggtitle(expression(paste(bar(r)[d], " per population sampled over xxx SNPs")))
ggia


ia_pop <- seppop(vcf_ia_gl ) %>% lapply(samp.ia) 

par(mfrow = c(2,5))
lapply(names(ia_pop), function(x) hist(ia_pop[[x]], main=x))


# for (i in 1:length(vcf_list)){
# 
#   write.vcf(vcf_list[[i]], file = paste(names(vcf_list)[[i]], "6262.vcf.gz", sep = "." ))
# 
# }

#vcf_list



```

## Minimum spanning network

```{r,echo=TRUE}

library(RColorBrewer)

pop(vcf.gatk.gl) <- pop$Country
# myCol <- brewer.pal(nPop(vcf.gatk.gl ), "Dark2") %>% setNames(popNames(vcf.gatk.gl ))
msn <- poppr.msn(vcf.gatk.gl , distmat = bitwise.dist(vcf.gatk.gl ), palette = rainbow)
set.seed(99)
plot_poppr_msn(vcf.gatk.gl , msn, inds = "nepal")

# library(pals)
# pop(vcf.gatk.gl ) <- pop$Area
# myCol <- alphabet(n = 25) %>% setNames(popNames(vcf.gatk.gl ))
# 
# 
# 
# myCol <- glasbey(n = 10) %>% setNames(popNames(vcf.gatk.gl ))
# msn <- poppr.msn(vcf.gatk.gl , distmat = bitwise.dist(vcf.gatk.gl ), palette = myCol)
# set.seed(99)
# plot_poppr_msn(vcf.gatk.gl , msn, inds = "nepal")




```




## Allele frequency distribution

```{r, eval=FALSE, include=FALSE}

myFreq <- glMean(vcf.gatk.gl )
myFreq <- c(myFreq, 1-myFreq)
hist(myFreq, proba=TRUE, col="darkseagreen3", xlab="Allele frequencies",
main="Distribution of allele frequencies", nclass=20)
# temp <- density(myFreq, bw=.05)
# lines(temp$x, temp$y*2,lwd=3)



```



## Population differentiation

```{r, echo=TRUE}

source("genetic_diff.R")

myDiff <- genetic_diff(vcf.gatk, pop$Country)

#dpf <- melt(myDiff[ , c(1, 24)] , na.rm=TRUE)
dpf <- melt(myDiff[ , c("CHROM", "Gst", "Gprimest")] , na.rm=TRUE)
p <- ggplot(dpf, aes(x=variable, y=value)) + geom_violin(fill="#8dd3c7")
p <- p + xlab("Population differentiation statistics")
p <- p + ylab("Value")
p <- p + theme_bw()
p


vcf_list <- vector("list", length(unique(pop$Country)))
names(vcf_list) <- unique(pop$Country)

grep_list <-  list("port", "taiw", c("queensaus|tasaus|waus|nswalesaus"), "safrica", "nvietnam", "newguin", "chile", "italy", "france", "domrep")

#grep_list <- list(c("port|france|italy"), c("taiw|nvietnam"), c("queensaus|tasaus|waus|nswalesaus"), "safrica", "newguin", "chile", "domrep")


for (i in (1:length(vcf_list))) {
  pop <- names(vcf_list[i])
  vcf_pop <- vcf.gatk
  gt <- vcf_pop@gt
  cols <- grep(grep_list[[i]], colnames(gt), ignore.case = TRUE)
    
  vcf_pop@gt <- gt[,c(1,cols)]
  
  vcf_list[[i]] <- vcf_pop


}

#saveRDS(vcf_list, file = "vcf_list.RData", compress = FALSE)


#vcf_list <- readRDS("vcf_list.RData")

Pairs <- function(vcf_list) {
  pop <- names(vcf_list)
  pairs <- t(combn(pop,2))
  colnames(pairs) <- c("pop1", "pop2")
  as.data.frame(pairs, stringsAsFactors = FALSE)
}

pop.pairs <- Pairs(vcf_list)

pop.pairs$Pop_Combination <- apply(pop.pairs[, c("pop1", "pop2")], 1, paste, collapse = "-")

out <- vector("list", nrow(pop.pairs))
names(out) <- pop.pairs$Pop_Combination

temp.vcf <- vcf_list[[1]]

source("genetic_diff.R")
for (i in 1:nrow(pop.pairs)) {
  
  pair <- unlist(pop.pairs[i, ])
  
  pop1 <- pair[1]
  names(pop1) <- NULL
  vcf1 <- vcf_list[[grep(pop1, names(vcf_list))]]
  samp_num1 <- ncol(vcf1@gt) -1
  pop_len1 <- rep(pop1, samp_num1)
  
  
  pop2 <- pair[2]
  names(pop2) <- NULL
  vcf2 <- vcf_list[[grep(pop2, names(vcf_list))]]
  samp_num2 <- ncol(vcf2@gt) -1
  pop_len2 <- rep(pop2, samp_num2)
  
  pop <- c(pop_len1, pop_len2)
  
  
  temp.vcf@gt <- cbind(vcf1@gt, vcf2@gt[, -1])
  
  myDiff <- genetic_diff(temp.vcf, pops = factor(pop))
  out[[i]] <- myDiff
  
  
}


pairwise_diff_mat <- function(out) {

pairwise_gst_list <- lapply(out, function(x) colMeans(x[11], na.rm = TRUE))
pop.pairs$Gstprime <- unlist(pairwise_gst_list)
pop.pairs

pop <- pop.pairs[, 3:4]

mat <- matrix(NA, nrow = length(vcf_list), ncol = length(vcf_list))
rownames(mat) <- names(vcf_list)
colnames(mat) <- names(vcf_list)
mat[lower.tri(mat)] <- pop$Gstprime
mat[upper.tri(mat)] <- t(mat)[upper.tri(mat)]
mat[is.na(mat)] <- 0

return(mat)

}

pairwise_mat <- pairwise_diff_mat(out)


library(reshape2)
pairwise_mat[upper.tri(pairwise_mat)] <- NA
lower_tri <- melt(pairwise_mat, na.rm = TRUE)

library(ggplot2)
ggheatmap <- ggplot(lower_tri, aes(x = Var1, y = Var2, fill = value)) + geom_tile(color = "white") +
  scale_fill_gradient(low = "green", high = "red" , space = "Lab", name="Pairwise GSTprime") + theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, vjust = 1, size = 12, hjust = 1)) + coord_fixed() +
  labs(x = "Population", y = "Population") +
  theme(axis.text.y = element_text(size = 12)) +
  theme(axis.title = element_text(size = 12)) +
  geom_text(aes(label = round(value, 2)))

ggheatmap



library(ape)
myNj <- ape::nj(as.dist(pairwise_mat))
myNj$edge.length[myNj$edge.length < 0] <- 0

plot(myNj, edge.width = 2, main = "NJ tree based on pairwise Gst")
#add.scale.bar()

plot(hclust(dist(pairwise_mat)))




```


## Number of heterozygous positions

```{r Heterozygosity, echo=TRUE}

gt <- extract.gt(vcf.gatk)
hets <- is.het(gt, na_is_false = TRUE)

sum_het <- colSums(hets)

sum_het.pop <- as.data.frame(sum_het)
sum_het.pop <- cbind(Isolate = rownames(sum_het.pop), Country = newpca_pop$Country, Heterozygosity = sum_het.pop)
rownames(sum_het.pop) <- NULL

ggplot(sum_het.pop, aes(x = Country, y = sum_het, fill = Country)) + geom_boxplot(outlier.shape = NA) + geom_jitter() +
  ggtitle("Number of heterozygous sites") +
  labs(x = "Country", y = "Number of heterozygous position") +
  theme(axis.text.x = element_text(size =12, angle = 60, hjust = 1)) +
  theme(axis.text.y = element_text(size =12)) + theme(legend.position = "none") +
  theme(plot.title = element_text(size = 30)) +
  theme(axis.title = element_text(size = 25))

```


## Neighbor joining tree

```{r, eval=FALSE, include=FALSE}

library(ggtree)
library(ggrepel)
library(poppr)
library(dplyr)

tree_nj <- aboot(vcf.gatk.gl, tree = "nj", distance = bitwise.dist, sample = 100, showtree = FALSE)

# library(pals)
# myCol <- glasbey(n = 10) %>% setNames(popNames(vcf.gatk.gl))
# 
# x <- apeBoot(tree_nj, tree_nj$node.label)
# strata <- cbind(dplyr::data_frame(taxa = indNames(vcf.gatk.gl)), country = pop$Country) 
# 
# test_tree <- ggtree(x, layout = "unrooted", branch.length = "branch.length") +
#   geom_label_repel(aes(label = bootstrap),
#   nudge_x = -0.01, nudge_y = 0) +
#   scale_size(range = c(2, 4)) +
#   theme(plot.title = element_text(size = 30))
# 
# 
# final_tree <- test_tree %<+% strata +
#               geom_tippoint(aes(color = country), alpha = 1, size = 3) +
#               theme(text = element_text(size = 18)) +
#               scale_color_manual(values = myCol) +
#               theme(legend.text = element_text(size = 20)) +
#               theme(legend.position = "right") 
#           
# final_tree

library(phangorn)
tree_nj <- midpoint(tree_nj)

countryInfo <- split(tree_nj$tip.label, pop$Country)
tree2 <- groupOTU(tree_nj, countryInfo)
tree2$tip.label <- as.character(pop$Area)
tree2$tip.label <- paste(pop$Area, pop$MT, sep = "_")

ggtree(tree2, aes(color=group), layout='circular') + geom_tiplab(size=4, aes(angle=angle)) 



#ggtree(tree2, aes(color=group), layout='rectangular') + geom_tiplab(size=4)

```