population genomics general script

popg_general.R

+#######################################
+## Population genetics : general script
+## J. Plantade
+## Jan 2021 
+#######################################
+
+# packages
+library(ade4)
+library(pegas)
+library(adegenet)
+library(PopGenome)
+library(hierfstat)
+library(ggplot2)
+library(poppr)
+library(rmarkdown)
+library(knitr)
+library(stringr)
+library(dplyr)
+
+################
+### setting data
+
+# data import into a loci object
+loci <- read.vcf("data_fixed_999/999_CD4_panTro3.vcf_fixed.vcf")
+nloci <- length(loci)
+# convert loci object into a genind object
+genind <- loci2genind(loci)
+# add population factor
+pop(genind)<-rep(c("paniscus","troglodytes_ellioti",
+                   "troglodytes_schweinfurthii","troglodytes_troglodytes","troglodytes_verus"), 
+                 c(11,10,6,4, 5))
+# convert genind object into a hierfstat object
+# hierf <-  genind2hierfstat(genind)
+
+##################
+### heterozygosity
+
+div <- summary(genind)
+locinfo<-VCFloci("data/CD4_panTro3.vcf_fixed.vcf", what = "all", chunk.size = 1e9, quiet = FALSE)
+# ggplot2 : Hobs-Hexp VS position
+g1 <- ggplot() + geom_point(aes(x=locinfo$POS, y=div$Hobs-div$Hexp), size=0.3) + xlab("Position") + ylab("Hobs - Hexp") + ggtitle("Difference between Hobs and Hexp for each position") + theme_light()
+g1
+
+# test Hobs and Hexp normality : Kolmogorov-Smirnov test for samples > 50
+hist(div$Hobs)
+ks.test(x = div$Hobs, y = "pnorm")
+hist(div$Hexp)
+ks.test(x = div$Hobs, y = "pnorm")
+# WARNING : if distributions are not normal we cannot use the Bartlett test (not robust enough) -> use the Levene test
+
+# BARTLETT TEST (homogeneity of variance between two samples : obs VS exp)
+# H0 : Var(Hobs) = Var(Hexp)
+bartlett.test(list(div$Hexp, div$Hobs)) # S : population is structured
+# theory : https://eric.univ-lyon2.fr/~ricco/cours/cours/Comp_Pop_Tests_Parametriques.pdf 
+
+# extract H values and bind them
+Hobs <- div$Hobs
+Hexp <- div$Hexp
+H <- cbind(Hobs, Hexp)
+H <- as.data.frame(H)
+# tidy the data : df -> one column with values, one column with the origin of values (observed or expected)
+library(tidyr)
+longer_H <- H %>%
+  pivot_longer("Hobs":"Hexp", names_to = "origin", values_to = "H_value")
+# LEVENE TEST
+library(lawstat)
+test <- levene.test(y = longer_H$H_value, group = longer_H$origin)
+
+# create one genind object for each population
+subspecies <- c("paniscus","troglodytes_ellioti","troglodytes_schweinfurthii","troglodytes_troglodytes","troglodytes_verus")
+pop <- c("pan", "ellioti", "schwein", "troglo", "verus")
+npop <- length(pop)
+for (i in 1:length(pop)){
+  tmp <- popsub(genind, sublist = c(subspecies[i]))
+  assign(pop[i], tmp)
+}
+# extracting heterozygosity for each population
+list_genind <- c(pan, ellioti, schwein, troglo, verus)
+for (i in 1:length(list_genind)){
+  tmp <- summary(list_genind[[i]])
+  assign(paste0(pop[i], "_div"), tmp)
+}
+# plot Hobs - Hexp for each population (for geom_rect see report_pan.Rmd)
+gg_pan <- ggplot() + 
+  geom_rect(aes(xmin = start, ymin = ymin, xmax = end, ymax = ymax), fill = 'grey') + 
+  geom_point(aes(x=locinfo$POS, y=pan_div$Hobs-pan_div$Hexp), size = 0.1) + 
+  labs(x = "position", y = "Hobs-Hexp", title = "P. paniscus") + 
+  theme_classic()
+gg_ellioti <- ggplot() + 
+  geom_rect(aes(xmin = start, ymin = ymin, xmax = end, ymax = ymax), fill = 'grey') + 
+  geom_point(aes(x=locinfo$POS, y=ellioti_div$Hobs-ellioti_div$Hexp), size = 0.1) + 
+  labs(x = "position", y = "Hobs-Hexp", title = "P. t. ellioti") + 
+  theme_classic()
+gg_schwein <- ggplot() + 
+  geom_rect(aes(xmin = start, ymin = ymin, xmax = end, ymax = ymax), fill = 'grey') + 
+  geom_point(aes(x=locinfo$POS, y=schwein_div$Hobs-schwein_div$Hexp), size = 0.1) + 
+  labs(x = "position", y = "Hobs-Hexp", title = "P. t. schweinfurthii") + 
+  theme_classic()
+gg_troglo <- ggplot() + 
+  geom_rect(aes(xmin = start, ymin = ymin, xmax = end, ymax = ymax), fill = 'grey') + 
+  geom_point(aes(x=locinfo$POS, y=troglo_div$Hobs-troglo_div$Hexp), size = 0.1) + 
+  labs(x = "position", y = "Hobs-Hexp", title = "P. t. troglodytes") + 
+  theme_classic()
+gg_verus <- ggplot() + 
+  geom_rect(aes(xmin = start, ymin = ymin, xmax = end, ymax = ymax), fill = 'grey') + 
+  geom_point(aes(x=locinfo$POS, y=verus_div$Hobs-verus_div$Hexp), size = 0.1) + 
+  labs(x = "position", y = "Hobs-Hexp", title = "P. t. verus") + 
+  theme_classic()
+
+#######
+### Fst
+# calculation of F statistics and extraction of Fst values
+# setting short names for each population
+pop <- c("pan", "ellioti", "schwein", "troglo", "verus")
+# Matrix X with colnames and rownames = pop, create the name of pairs
+npop <- 5
+M <- matrix(NA, nrow = npop, ncol = npop)
+M <- as.data.frame(M)            
+rownames(M) <- pop            
+colnames(M) <- pop
+for (i in 1:npop){
+  start <- colnames(M)[i]
+  for (j in i+1:npop-i){
+    end <- rownames(M)[j]
+    M[j, i] <- paste0(start, "_", end)
+  }
+}
+# creates a list of pair's name
+pair_name <- c(rep("NA", (npop^2-npop)/2))
+k <- 1
+for (i in 1:(npop-1)){
+  for (j in i+1:(npop-i)){
+    pair_name[k] <- M[j, i]
+    k <- k+1
+  }
+}
+# genind subset -> pair of populations.
+pair_list <- NULL
+for (i in 1:(npop-1)){ # in matrix M we go from col 2 to 5
+  for (j in (i+1):(npop)){ # in matrix M we go from row 3 to 6
+    tmp <- popsub(genind, sublist = c(subspecies[i], subspecies[j]))
+    tmp2 <- assign(M[j, i], tmp)
+    pair_list <- c(pair_list, tmp2) # list of all pop pairs
+  }
+}
+# calculation of F statistics and extraction of Fst values for each pair
+fst <- matrix(NA, nrow = nloci, ncol = length(pair_list))
+for (i in 1:length(pair_list)){
+  tmp <- Fst(as.loci(pair_list[[i]]), na.alleles = "")
+  fst[,i] <- tmp[,2]
+} 
+fst <- as.data.frame(fst)
+colnames(fst) <- pair_name
+# replace all NaN values by proper NA values in fst data frame
+fst[fst == "NaN"] <- NA
+# create a list containing every plot of fst values
+gg <- lapply(1:length(pair_name), function(i) {
+  pair <- pair_name[i]
+  title <- str_replace(pair, pattern = "_", replacement = "/")
+  ggplot() + geom_point(aes(x=locinfo$POS, y=fst[,i]), size=0.3) + xlab("position") + ylab("Fst value") + ggtitle(title) + theme_light()
+})
+#######################
+## Clustering with DAPC 
+
+loci <- read.vcf("data_fixed_999/999_CD4_panTro3.vcf_fixed.vcf")
+genind <- loci2genind(loci)
+### Looking for 5 clusters = species/subspecies
+pop(genind)<-rep(c("paniscus","troglodytes_ellioti",
+                   "troglodytes_schweinfurthii","troglodytes_troglodytes","troglodytes_verus"), 
+                 c(11,10,6,4, 5))
+
+# Identify the number of cluster
+grp <- find.clusters(genind, max.n.clust = 10) #PCs to retain = 35; nb of clusters = 4
+
+# Kstat : BIC for each value of k, here k=4 has the lower BIC
+grp$Kstat
+# stat : gives the selected nb of clusters (selected k) and its BIC value
+grp$stat
+# grp : gives the id of the group to which each individual is assigned
+head(grp$grp)
+# size : number of individuals assigned to each group
+grp$size
+
+# extract the group partition
+group_df <- as.data.frame(grp$grp)
+# add the infectious status for each individual
+inf_status <- rep(c("nonSIV", "SIV","nonSIV"), c(21,10, 5))
+size <- rep(0, 36)
+tab <- cbind(group_df, pop(genind), inf_status, size)
+colnames(tab) <- c("group", "population", "infectious_status", "size")
+# set point size in function of the size of the group
+grp_size <- grp$size
+for (i in 1:length(tab$group)){
+  for (j in 1:6){
+    if (tab[i, 1]==j){
+      tab[i,4] <- grp_size[j]
+    }
+  }
+}
+# plot the groups with population, infectious status and group size
+gg <- ggplot(data = tab, aes(x = group, y = population, colour = infectious_status, size = size)) + geom_point()
+gg
+
+# decompose the find.cluster function
+kmean <- kmeans(x = genind@tab, centers = 4)
+kmeanBIC <- function (clust) {
+  m = ncol(clust$centers)
+  n = length(clust$cluster)
+  k = nrow(clust$centers)
+  D = clust$tot.withinss
+  return(BIC = D + log(n)*m*k)
+}
+# save the BIC values in a df
+k <- c(seq(1,10,1))
+BIC <- c(rep("NA", 10))
+bic_value <- cbind(k, BIC)
+for (i in 1:10){
+  # cluster the data with k=i
+  clust <- kmeans(x = genind@tab, centers = i)
+  # compute the BIC value for k=i
+  bic <- kmeanBIC(clust)
+  bic_value[i,2] <- bic
+} # return wird values
+
+# describe clusters using dapc1
+dapc1 <- dapc(genind, grp$grp) # PCs=10, Linear discriminant=5
+# visualization of clusters
+myCol <- c("darkblue", "purple", "green", "orange")
+scatter(dapc1, scree.da = FALSE, bg="white", pch=20, cell=0, cstar=0, col=myCol, solid=0.4, cex=2, clab=0, leg=TRUE, txt.leg = (paste("Cluster", 1:4)))
+
+### looking for infected/non infected populations
+# reload the genind object
+# data import into a loci object
+
+# then run the dapc again
+
+##########
+## Tajimas'D with popGenome
+
+# create a GENOME class object with vcf file
+genome <- readData(path = "popgenome_data/", format = "VCF")
+
+### Marie's method : setting pop considering that individuals are diploid
+# extract the names of individuals in the genome object
+pops <- get.individuals(genome)[[1]]
+# subset populations
+# grep : search for matches of a pattern and return the indices that yielded a match
+# by putting the indices in [] we extract the individuals names
+paniscus <- pops[grep("Pan_paniscus",pops)]
+ellioti <- pops[grep("Pan_troglodytes_ellioti",pops)]
+schwein <- pops[grep("Pan_troglodytes_schweinfurthii",pops)]
+troglo <- pops[grep("Pan_troglodytes_troglodytes",pops)]
+verus <- pops[grep("Pan_troglodytes_verus",pops)]
+# we set the populations with a list
+genome  <- set.populations(genome, list(paniscus, troglo, ellioti, schwein, verus))
+# check : two copies for each individual because they are diploid
+genome@populations
+###
+
+# set the outgroup (it changes the D values)
+genome <- set.outgroup(genome, paniscus)
+# check
+genome@outgroup
+
+# neutrality tests
+genome <- neutrality.stats(genome)
+# get the Tajima's D values
+genome@Tajima.D
+# pop :
+# 1 > paniscus
+# 2 > troglo
+# 3 > ellioti
+# 4 > schwein
+# 5 > verus
+
+# figure of D values for one gene
+library(tidyr)
+taj <- genome@Tajima.D
+taj <- as.data.frame(t(taj))
+taj <- cbind(taj, pop = c("paniscus", "t. troglodytes", "t. ellioti", "t. schweinwurthii", "t. verus"))
+colnames(taj)[1] <- c("d_value")
+taj <- taj %>% drop_na()
+
+gg <- ggplot(data = taj, mapping = aes(x = d_value, y = pop)) + geom_point() + labs(title = "Tajima's D value for each P. troglodytes subspecies", x = "D value", y = "populations")
+gg