Script general

popg_general_01262021.R

+#######################################
+## Population genetics : general script
+## J. Plantade
+## Jan 2021 
+#######################################
+
+# packages
+library(ade4)
+library(pegas)
+library(adegenet)
+library(PopGenome)
+library(hierfstat)
+library(ggplot2)
+library(poppr)
+
+################
+### setting data
+
+# data import into a loci object
+loci <- read.vcf("data/CD4/CD4_panTro3.vcf_fixed.vcf", from = 1, to = 572)
+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
+# Bartlett test (homogeneity of variance)
+# H0 : Hobs = Hexp
+bartlett.test(list(div$Hexp, div$Hobs)) # S : population is structured
+# 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)
+}
+# -> possibility to to Bartlett test and plot Hobs for each pop
+
+#######
+### Fst
+
+# 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 names
+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 -> couple of population.
+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
+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
+
+# Fst values VS loci position
+plot(locinfo$POS, fst[,1], cex=0.2, xlab="Position", ylab="Fst value", main=colnames(fst)[1])
+abline(0.15, 0)
+
+#############
+## Clustering with DAPC 
+
+# 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
+
+# Compare inferred groups (inf) with actual groups (ori)
+# Find actual groups with function pop on the genind object
+pop(genind)
+tab <- table(pop(genind), grp$grp)
+# Display
+table.value(tab, col.labels = colnames(tab), row.labels = rownames(tab))
+
+# 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=1, clab=0, leg=TRUE, txt.leg = (paste("Cluster", 1:4)))
+
+##########
+## Tajimas'D with popGenome
+
+# list of populations
+populations <- list("paniscus","troglodytes_ellioti","troglodytes_schweinfurthii","troglodytes_troglodytes","troglodytes_verus")
+# create a GENOME class object with vcf file
+genome <- readData(path = "data/CD4/", format = "VCF")
+# function that create a list, which element is the ensemble of individual for each pop
+listing_ind <- function(npan, nellio, nschwein, ntroglo, nverus){
+  pan_vec <- c(rep("NA", npan))
+  ellio_vec <- c(rep("NA", nellio))
+  schwein_vec <- c(rep("NA", nschwein))
+  troglo_vec <- c(rep("NA", ntroglo))
+  verus_vec <- c(rep("NA", nverus))
+  for (i in 1:npan){
+    pan_vec[i] <- paste0("pan_", i)
+  }
+  for (i in 1:nellio){
+    ellio_vec[i] <- paste0("ellioti_", i)
+  }
+  for (i in 1:nschwein){
+    schwein_vec[i] <- paste0("schwein_", i)
+  }
+  for (i in 1:ntroglo){
+    troglo_vec[i] <- paste0("troglo_", i)
+  }
+  for (i in 1:nverus){
+    verus_vec[i] <- paste0("verus_", i)
+  }
+  populations <- list(pan_vec, ellio_vec, schwein_vec, troglo_vec, verus_vec)
+}
+# apply listing_ind with the rignt number of individuals for each pop
+# and save it in the object 'populations'
+populations <- listing_ind(11, 10, 6, 4, 5)
+# set populations in the genome object
+genome <- set.populations(genome, new.populations = populations, diploid = TRUE)
+# check
+genome@populations
+# neutrality tests
+genome <- neutrality.stats(genome) # error