library("BPEC") ## To load the Brown Frog data through BPEC: data("MacrocnemisRawSeqs", package = "BPEC") data("MacrocnemisCoordsLocs", package = "BPEC") coordsLocs <- MacrocnemisCoordsLocs rawSeqs <- MacrocnemisRawSeqs ## Instead, to load the Brown Frog data through the source files directly: rawSeqs <- bpec.loadSeq("haplotypes.nex") coordsLocs <- bpec.loadCoords("coordsLocsFile.txt", header = TRUE) ## To run the analysis: set.seed(108) bpecout <- bpec.mcmc(rawSeqs, coordsLocs, maxMig = 3, iter = 1000000, ds = 3, postSamples = 1000, dims = 8) ## To produce the visualization plots: ################################ Figure 6 of the paper################################# ## Unlist the sequences: rawSeqsList <- unlist(rawSeqs) rawSeqsList <- matrix(rawSeqsList, nrow = length(rawSeqs[[1]]), ncol = length(rawSeqs)) rawSeqsList <- t(rawSeqsList) ## Identify any nucleotide sites which are identical for all sequences: delCol <- NULL for (i in 1:ncol(rawSeqsList)) { if (length(unique(rawSeqsList[, i])) == 1) { delCol <- c(delCol, i) } } rm(rawSeqsList) ## Remove such nucleotide sites: if (is.null(delCol) == FALSE) { for (i in 1:length(rawSeqs)) { rawSeqs[[i]] <- rawSeqs[[i]][-delCol] } } ## Camerani sequences are characterised by a 'C' nucleotide in the 3rd ## SNP site of the sequences: camerani <- numeric(length(rawSeqs)) for (i in 1:length(rawSeqs)) { if(rawSeqs[[i]][3] == "c"){ camerani[i] <- 1 } } ## In each location, count how many sequences within are camerani, ## divide by total number of sequences in each location. The ## geographical information is available per sequence, so the first ## instance of each location will contain the camerani percentage, ## whereas any duplicate locations will be stored as '-1' (but with ## their sequence information included in the first instance). cameraniperloc <- camerani done <- numeric(length(rawSeqs)) for (i in 1:length(rawSeqs)) { nlocs <- 1 for (j in 1:length(rawSeqs)) { if (done[j] == 0 && i!=j && coordsLocs[i, 1] == coordsLocs[j, 1] && coordsLocs[i, 2] == coordsLocs[j, 2]) { cameraniperloc[i] <- cameraniperloc[i] + cameraniperloc[j] nlocs <- nlocs + 1 done[j] <- 1 cameraniperloc[j] <- -1 } } if (nlocs > 1) { cameraniperloc[i] <- round(cameraniperloc[i]/nlocs, 2) } done[i] <- 1 } ## Plot all unique locations, where the colour is given by the ## percentage of camerani: savex <- coordsLocs[cameraniperloc > -0.5, 2] savey <- coordsLocs[cameraniperloc > -0.5, 1] savecol <- rgb(cameraniperloc[cameraniperloc > -0.5], 0, 0, max = 1) style1 <- paste("feature:all|element:labels|visibility:off", "&style=feature:road|element:all|visibility:off", "&style=feature:administrative|element:geometry|visibility:off", "&style=feature:transit.line|element:geometry|visibility:off", "&style=feature:water|element:geometry.fill|color:0xffffff", sep = "") totalMap <- ggmap::get_googlemap(c((min(savex) + max(savex))/2, (min(savey) + max(savey))/2), zoom = 7, maptype = "terrain", color = "bw", style = style1, messaging = FALSE) contourMap <- ggmap::ggmap(totalMap) savexy <- data.frame(longitude = savex, latitude = savey) contourMap <- contourMap + ggplot2::geom_point(data = savexy, ggplot2::aes(x = longitude,y = latitude), size = 1.4, color = savecol) contourMap <- contourMap + labs(x = "longitude", y = "latitude") plot(contourMap) ################################ Figure 7 of the paper################################# pdf("BPECcontourPlot.pdf") par(mfrow = c(1, 2)) bpec.contourPlot(bpecout, GoogleEarth = 0, mapType = "google", colorCode = c(7, 5, 6, 3, 2), mapCentre = NULL, zoom = 7) bpec.contourPlot(bpecout, GoogleEarth = 0, mapType = "osm", colorCode = c(7, 5, 6, 3, 2), mapCentre = NULL, zoom = 7) dev.off() ################################ Figure 9 of the paper################################# pdf("BPECcovariatesPlot.pdf") par(mfrow = c(2, 3)) bpec.covariatesPlot(bpecout) dev.off() ################################ Figure 10 of the paper################################# pdf("BPECtreePlot.pdf") par(mfrow = c(1, 1), ask = TRUE) bpec.tree <- bpec.treePlot(bpecout, colorCode = c(7, 5, 6, 3, 2)) bpec.geo <- bpec.geoTree(bpecout, file = "GoogleEarthTree.kml") dev.off()