Pinf_diversity_MX.Rmd

---
title: "P.infestans diversity 2015/16 in Mexico"
author: "Shankar K Shakya"
date: "August 22, 2017"
output: html_document
---

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


## Summary

This document includes R scripts for reproducibility.


## Reading data

First step was to read the data and filter locus for missing data (D13 locus was removed). Only strictly diploid isolates were analyzed. 
# Per cent missing data by population x locus

```{r}

rm(list = ls())

library(poppr, quietly = TRUE)
P.inf <- read.genalex("MX_pop.csv", ploidy = 4)

source("multiploid2diploid.R")
diploid_pinf <- multiploid2diploid(P.inf, to = 2)
info_table(diploid_pinf, type = "missing", plot = T)

diploid_pinf <- missingno(diploid_pinf, type = "loci", cutoff = 0.07)

splitStrata(diploid_pinf) <- ~Region/Year
setPop(diploid_pinf) <- ~Region

```

## Table 2. Diversity summary by region 

```{r}

# Corrected simpson diversity index
uSimp <- function(x){
  lambda <- vegan::diversity(x, "simpson")
  x <- drop(as.matrix(x))
  if (length(dim(x)) > 1){
    N <- rowSums(x)
  } else {
    N <- sum(x)
  }
  return((N/(N-1))*lambda)
}

unbiased <- poppr(diploid_pinf, uSimp = uSimp)

unbiased <- unbiased[, c(1,2,3,4,9,10,11)]
colnames(unbiased) <- c("Region", "N", "MLG", "eMLG", "E5", "Corrected lambda", "Hexp")
unbiased <- unbiased[c(6,1,4,3,2,5), ]
rownames(unbiased) <- NULL

library(knitr)
kable(unbiased, digits = 2,
      align = c("c"),
      caption = "Table 2. P. infestans sampling locations in Mexico and population statistics for diploid genotypes based on eleven SSR loci.  Total number of samples (N), observed multilocus genotype (MLG), expected multilocus genotype (eMLG), eveness (E5), corrrected Simpson's diversity index (Corrected lambda) and expected heterozygosity (Hexp)")

```


## Table 3. rbarD and Hardy-Weinberg equilibrium

```{r}

ia_Pinf <- poppr(diploid_pinf, clonecorrect = TRUE, strata = ~Region, sample = 999, quiet = TRUE)
ia_Pinf <- ia_Pinf[,c(1:3,13:14)]  
colnames(ia_Pinf) <- c("Region", "N", "MLG", "rbarD", "P-value")

library(pegas)

diploid_pinf_cc <- clonecorrect(diploid_pinf)
hwe_per_pop <- seppop(diploid_pinf_cc) %>% lapply(hw.test, B = 1000)
per_pop_mat <- sapply(hwe_per_pop, "[", i = TRUE, j = 3)
alpha <- 0.05
per_pop_mat[per_pop_mat > alpha] <- 1


library(lattice)

hwe_plot <- levelplot(per_pop_mat, xlab = "Locus", ylab = "Population")

num_loci_hwe <- per_pop_mat
num_loci_hwe[num_loci_hwe != 1] <- 0
num_loci_hwe <- colSums(num_loci_hwe)


ia_Pinf <- ia_Pinf[-7, ]
ia_Pinf <- cbind(ia_Pinf, num_loci_hwe) 
ia_Pinf <- ia_Pinf[c(6,1,4,3,2,5), ]
rownames(ia_Pinf) <- NULL
colnames(ia_Pinf[6]) <- "Number of loci under HWE"
 
kable(ia_Pinf, 
      digits = 2,
      align = c("c"),
      caption = "Table 3. Standardized index of association (rbarD), a measure of linkage disequilibrium among loci, and number of loci under Hardy-Weinberg equilibrium on diploid clone corrected dataset of P.infestans based on eleven SSR markers")


```


## Table 4. Pairwise FST 

```{r, message=FALSE, warning=FALSE}

diploid_pinf_cc <- clonecorrect(diploid_pinf)

library(strataG)
gen2gtype <- function (x) {
  
  gen.mat <- genind2df(x, usepop = TRUE, oneColPerAll = TRUE)
  strata <- x@strata
  rownames(strata) <- rownames(gen.mat)
  x@strata <- strata
  gen.mat[gen.mat == "NA"] <- NA
  has.pop <- !is.null(x@pop)
  df2gtypes(x = gen.mat, ploidy = x@ploidy[1], id.col = NULL, 
            strata.col = if (has.pop) 
              1
            else NULL, loc.col = if (has.pop) 
              2
            else 1, 
            schemes = x@strata,  other = other(x))
}

diploid_pinf_cc.g <- gen2gtype(diploid_pinf_cc)

pairwise_fst <- pairwiseTest(diploid_pinf_cc.g, stats = "fst", nrep = 10000)

fst_mat <- pairwise_fst$pair.mat$Fst
fst_mat[upper.tri(fst_mat)] <- t(fst_mat)[upper.tri(fst_mat)]
fst_mat[is.na(fst_mat)] <- 0


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


kable(fst_mat, digits = 2,
      caption = "Table 4. Pairwise FST values (lower triangle) for diploid P.infestans genotypes from regions of Mexico using eleven SSR markers. Upper triangle represents the associated P-value for pairwise comparisons based on 10,000 permutation replicates.")


# 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 FST") + 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

```


## Table 5. Analysis of molecular variance

```{r}

repeat_length <- c(2, 2, 2, 3, 3, 2, 2, 2, 2, 2, 2)
bruvodist <- bruvo.dist(clonecorrect(diploid_pinf), replen = repeat_length)

amv <- poppr.amova(diploid_pinf, 
                   hier = ~Region, 
                   clonecorrect = TRUE, 
                   dist = bruvodist,
                   within = F)


myamv <- amv$results
rownames(myamv) <- c("Between regions", "Within regions", "Total")
myamv <- cbind(myamv, amv$componentsofcovariance$`%`)

colnames(myamv) <- c("df", "SS", "MSS", "% variance")

knitr::kable(myamv, digits = 2,
             caption = "Table 5. Analysis of Molecular Variance for clone corrected diploid P. infestans based on Bruvo’ s distance ")


```




## Fig. 1. Sampling states and sites in Mexico.

```{r, warning=FALSE}
library(ggplot2)
library(raster)
library(gridExtra)
library(grid)
library(maptools)

mx <- readShapePoly("mexstates/mexstates.shp")
mystates <- mx[grep("Mexico|Tlaxcala|Michoacan", mx@data$ADMIN_NAME),]

states <- coordinates(mystates)
states <- data.frame(states)
states$label <- mystates@data$ADMIN_NAME

site <- read.csv("Dgeo_MX.csv", header = TRUE)
colnames(site) <- c("label", "X2", "X1")
states <- rbind(states, site[c(3,1,2)])
states <- states[4:9, ]
colnames(states)[3] <- "Sites"

pol <- data.frame(xmin=-103.7455,xmax=-97.61348 ,ymin=17.92083 ,ymax=20.39417)
mystates <- fortify(mystates)
colnames(mystates)[6] <- "States"
mystates$States[mystates$States == 18] <- "Michoacan"
mystates$States[mystates$States == 19] <- "Central Mexico"
mystates$States[mystates$States == 20] <- "Tlaxcala"

mycol <- rainbow(6)
mycol[2] <- "#000000"

mx_map <- ggplot()+
  geom_polygon(data=mx, aes(long,lat,group=group),colour="grey10",fill="#fff7bc") +
  coord_equal()+
  theme_bw()+labs(x="Longitude",y = "Latitude") + 
  geom_rect(data = pol, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), alpha=0, colour="red", size = 1, linetype = 1) +
  theme(axis.text = element_text(size =10, face = "bold", family = "Microsoft Sans Serif", colour = "black")) +
  theme(axis.title = element_text(size = 10, face = "bold", family = "Microsoft Sans Serif", colour = "black"))

#ggsave("insetmx_map.tiff", plot = last_plot(), height = 3, width = 4, dpi = 600)

mx_map

statemap <- ggplot() + 
  geom_polygon(data=mystates, aes(x = long, y = lat, group=group, fill = States)) +
  coord_equal() + 
  theme_bw() + labs(x="Longitude", y = "Latitude") +
  geom_point(data = states, aes(x=X1, y=X2, color = Sites), size = 1.9) + 
  scale_color_manual(values = mycol) +
  theme(axis.text = element_text(size = 10, face = "bold", family = "Microsoft Sans Serif", colour = "black")) +
  theme(axis.title = element_text(size = 10, face = "bold", family = "Microsoft Sans Serif")) +
  theme(text = element_text(size =10, face = "bold", family = "Microsoft Sans Serif", colour = "black")) +
  theme(plot.margin=unit(c(0,0,0,0),"mm")) 

statemap

#ggsave("statemap.tiff", plot = last_plot(), width = 7, height = 3.5, dpi = 600)



```

## Fig 2. Neighbor-joining tree

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

source("https://bioconductor.org/biocLite.R")
biocLite("ggtree", type = "source")

library(ggtree)
library(ggrepel)

repeat_length <- c(2, 2, 2, 3, 3, 2, 2, 2, 2, 2, 2)

diploid_pinf <- clonecorrect(diploid_pinf)
tree_nj <- bruvo.boot(diploid_pinf, replen = repeat_length, 
                   sample = 1000, cutoff = 70, tree = "nj")

tree_nj <- readRDS("tree_nj.RData")
tree_nj$node.labels[tree_nj$node.labels > 70]  <- "*" 

mycol <- rainbow(6)
mycol[2] <- "#000000"
mycol <- mycol %>% setNames(popNames(diploid_pinf))

x <- apeBoot(tree_nj, tree_nj$node.label)
strata <- dplyr::data_frame(taxa = indNames(diploid_pinf)) %>% cbind(strata(diploid_pinf)) 
#tree_nj$tip.label <- as.character(strata$Region) 
# 
# plot(tree_nj,
#             tip.color = mycol, show.node.label = TRUE, use.edge.length = F,
#             edge.width = 2,
#             align.tip.label = T,
#             font = 2, adj = 1, label.offset = 1,
#             main = "NJ tree.", cex = 0.65, node.pos = 1, type = "unrooted")


test_tree <- ggtree(x, layout = "circular", branch.length = "none") +
  geom_text(aes(label = bootstrap), color = "red")
  
  
final_tree <- test_tree %<+% strata +
              geom_tippoint(aes(color = Region), alpha = 1, size = 1) +
              scale_color_manual(values = mycol) +
              theme(legend.text = element_text(size = 10)) +
              theme(legend.position = "right") + 
              theme(text = element_text(size =10, face = "bold", family = "Microsoft Sans Serif", colour = "black"))

```


## Fig 3 & 4. Delta K plot and STRUCTURE plot

```{r}

library(strataG)
# Sys.setenv(PATH = paste("C:/Program Files (x86)/Structure2.3.4/bin", Sys.getenv("PATH"), sep = ";"))
Sys.setenv(PATH = paste("C:/Users/Shankar/Downloads/CLUMPP_Windows.1.1.2/", Sys.getenv("PATH"), sep = ";"))

# gen2gtype <- function (x) {
#   
#   gen.mat <- genind2df(x, usepop = TRUE, oneColPerAll = TRUE)
#   strata <- x@strata
#   rownames(strata) <- rownames(gen.mat)
#   x@strata <- strata
#   gen.mat[gen.mat == "NA"] <- NA
#   has.pop <- !is.null(x@pop)
#   df2gtypes(x = gen.mat, ploidy = x@ploidy[1], id.col = NULL, 
#             strata.col = if (has.pop) 
#               1
#             else NULL, loc.col = if (has.pop) 
#               2
#             else 1, 
#             schemes = x@strata,  other = other(x))
# }
# 
# diploid_pinf.g <- gen2gtype(diploid_pinf)
# 
# sr <- structureRun(diploid_pinf.g, k.range  = 1:6, burnin = 20000, num.k.rep = 15, noadmix = F, numreps = 100000)
#save(sr, file = "diploid_pinf_str_byregion_exceptD13locus_admix.RData")

str_data <- load("diploid_pinf_str_byregion_exceptD13locus_admix.RData")
sr <- get(str_data)

source("deltaK_plot.R")
evno <- deltaK_plot(sr, plot = TRUE)

# qmat2 <- clumpp(sr, k = 2, sim.stat = "g.prime", repeats = 1000)
# qmat3 <- clumpp(sr, k = 3, sim.stat = "g.prime", repeats = 1000)
# qmat4 <- clumpp(sr, k = 4, sim.stat = "g.prime", repeats = 1000)
# qmat5 <- clumpp(sr, k = 5, sim.stat = "g.prime", repeats = 1000)
# 
# qmat_list <- list(qmat2, qmat3, qmat4, qmat5)
# saveRDS(qmat_list, "qmat_list.RData")

edited_qmat_list <- readRDS("qmat_list.RData")

for (i in 1:length(edited_qmat_list)){
  temp_mat <- edited_qmat_list[[i]]
  key <- c("Michoacan" = "MICHO", "Chapingo" = "CHA", "Toluca" = "TOLU", "San Gerónimo" = "SG", "Juchitepec" = "JFH", "Tlaxcala" = "TLAX")
  temp_mat$orig.pop <- factor(key[temp_mat$orig.pop], ordered = TRUE, levels = key)
  edited_qmat_list[[i]] <- temp_mat
} 

source("structurePlot.R")
plots <- lapply(edited_qmat_list, structure.Plot, horiz = FALSE)

library(cowplot)
plot_grid(plotlist = plots, nrow = 4, ncol = 1, hjust = 1) +
theme(panel.border = element_rect(colour = "blue", fill=NA, size=3)) +
theme(plot.margin=unit(c(0,0,0,0),"mm"))


```

## Fig 5. Discriminant Analysis of Principal components

```{r}

set.seed(999)
P.infx <- xvalDapc(tab(diploid_pinf, NA.method = "mean"), pop(diploid_pinf),
                   n.pca = 60:70,  n.rep = 100, 
                   parallel = "multicore", ncpus = 4L)

library(RColorBrewer)
myCol <- brewer.pal(nPop(diploid_pinf), "Dark2") %>% setNames(popNames(diploid_pinf))

diploid_pinf_no_MCH <- popsub(diploid_pinf, blacklist = "Michoacan")


set.seed(999)
P.infx1 <- xvalDapc(tab(diploid_pinf_no_MCH, NA.method = "mean"), pop(diploid_pinf_no_MCH),
                   n.pca = 60:70,  n.rep = 100,
                   parallel = "multicore", ncpus = 4L)

source("scatterdapc.R")

par(mfrow = c(1,2))
par(mar = c(1,1,1,3))

scatterdapc(P.infx$DAPC, col=myCol, clabel = 0.75, pch=15:19, scree.pca = TRUE, scree.da = FALSE, 
        posi.pca = "topright", posi.leg = "topright", legend = F, 
        cleg = 0.9, inset.solid = 1, xax = 1, yax = 2, cex.lab = 1, cex = 1, solid = 1, cstar = 0)
title(main = "A", col = "black", font.main = 2, adj = 0, line = -1)

scatterdapc(P.infx1$DAPC, col=myCol, clabel = 0.75, pch=15:19, scree.pca = TRUE, scree.da = FALSE, 
        posi.pca = "topright", posi.leg = "topright", legend = F, 
        cleg = 0.9, inset.solid = 1, xax = 1, yax = 2, cex.lab = 1, cex = 1, solid = 1, cstar = 0)
title(main = "B", col = "black", font.main =2, adj = 0, line = -1)



```

## Fig 6. Isolation by distance

```{r}

lat_lon <- read.csv("Dgeo_MX.csv")
colnames(lat_lon) <- c("name", "lat", "lon")
source("GeoDistanceInMetresMatrix.R")

geodist_mat <- GeoDistanceInMetresMatrix(lat_lon) / 1000
geodist_mat <- as.dist(geodist_mat)


new_fst_mat <- readRDS("new_fst_mat.RData")
fstdist_mat <- as.dist(new_fst_mat)

plot(geodist_mat, fstdist_mat, xlab = "Geographical distance (km)", ylab = "Genetic distance (FST)",
     pch = 19, bg = "black", cex.lab = 1.5, cex.axis = 1.5, font.lab = 2, 
           font = 2, family = "Microsoft Sans Serif")
abline(lm(fstdist_mat~geodist_mat))

summary(lm(fstdist_mat~geodist_mat))


```