inst/extdata/make_eems_plots_JLW.R
In JeffWeinell/misc.wrappers: Miscellaneous Wrapper Functions for Spatial Population Genomics
#####
## All of these functions are from the R package 'reemsplots2' by Desislava Petkova (github repository 'dipetkov/reemsplots2').
## Only very minor edits were been made by me (Jeff Weinell), for a few of the functions, but in most cases the functions are identical to those in reemsplots2.
#' Plot effective migration and diversity surfaces
#'
#' Given a vector of EEMS output directories, this function generates several figures to visualize EEMS results. It is a good idea to examine all these figures, which is why they are generated by default.
#' \describe{
#' \item{mrates01}{effective migration surface. This contour plot visualizes the estimated effective migration rates \code{m}, on the log10 scale after mean centering.}
#' \item{mrates02}{posterior probability contours \code{P(log(m) > 0) = p} and \code{P(log(m) < 0) = p} for the given probability level \code{p}. Since migration rates are visualized on the log10 scale after mean centering, 0 corresponds to the overall mean migration rate. This contour plot emphasizes regions with effective migration that is significantly higher/lower than the overall average.}
#' \item{qrates01}{effective diversity surface. This contour plot visualizes the estimated effective diversity rates \code{q}, on the log10 scale after mean centering.}
#' \item{qrates02}{posterior probability contours \code{P(log(q) > 0) = p} and \code{P(log(q) < 0) = p}. Similar to \code{mrates02} but applied to the effective diversity rates.}
#' \item{rdist01}{scatter plot of the observed vs the fitted between-deme component of genetic dissimilarity, where one point represents a pair of sampled demes.}
#' \item{rdist01}{scatter plot of the observed vs the fitted within-deme component of genetic dissimilarity, where one point represents a sampled deme.}
#' \item{rdist03}{scatter plot of observed genetic dissimilarities between demes vs observed geographic distances between demes.}
#' \item{pilogl01}{posterior probability trace}
#' }
#' @details
#' The function \code{make_eems_plots} will work given the results from a single EEMS run (one directory in \code{mcmcpath}) but it is better to run EEMS several times, randomly initializing the MCMC chain each time. In other words, simulate several realizations of the Markov chain and let each realization start from a different state in the parameter space (by using a different random seed).
##'
#' The \code{mrates} and \code{qrates} figures visualize (properties of) the effective migration and diversity rates across the habitat. The other figures can help to check that the MCMC sampler has converged (the trace plot \code{pilogl}) and that the EEMS model fits the data well (the scatter plots of genetic dissimilarities \code{rdist}).
#'
#' To describe the within-deme and between-deme components of genetic dissimilarity, let \code{D(a,b)} be the dissimilarity between one individual from deme \code{a} and another individual from deme \code{b}. Then the within-deme component for \code{a} and \code{b} is simply \code{D(a,a)} and \code{D(b, b)}, respectively. The between-deme component is \code{D(a,b) - [D(a,a) + D(b,b)] / 2} and it represents dissimilarity that is due to the spatial structure of the population and is not a consequence of the local diversity in the two demes.
#' @param mcmcpath A vector of EEMS output directories, for the same dataset. Warning: There is minimal checking that the directories all correspond to the same dataset.
#' @param longlat A logical value indicating whether the coordinates are given as pairs (longitude, latitude) or (latitude, longitude).
#' @param dpi Resolution of the contour raster. The default is 250.
#' @param add_grid A logical value indicating whether to add the population grid or not.
#' @param col_grid The color of the population grid. Defaults to gray.
#' @param add_outline A logical value indicating whether to add the habitat outline or not.
#' @param col_outline The color of the habitat outline. Defaults to white.
#' @param add_demes A logical value indicating whether to add the observed demes or not.
#' @param col_demes The color of the demes. Defaults to black.
#' @param eems_colors The EEMS color scheme as a vector of colors, ordered from low to high. Defaults to a DarkOrange to Blue divergent palette with six orange shades, white in the middle, six blue shades. Acknowledgement: The default color scheme is adapted from the \code{dichromat} package.
#' @param m_colscale A fixed range for log10-transformed migration rates. If the estimated rates fall outside the specified range, then the color scale is ignored. The default range is \code{[-2.5, +2.5]}.
#' @param q_colscale A fixed range for log10-transformed diversity rates. The default range is \code{-0.1, +0.1}.
#' @param add_abline Add the line \code{y = x} to the diagnostic scatter plots of observed vs fitted genetic dissimilarities.
#' @param prob_level A probability \code{p} to define the posterior probability contours \code{P(log(m) > 0) = p} and \code{P(log(m) < 0) = p}. Defaults to \code{0.9}.
#' @references Light A and Bartlein PJ (2004). The End of the Rainbow? Color Schemes for Improved Data Graphics. EOS Transactions of the American Geophysical Union, 85(40), 385.
#' @examples
#' # Use the provided example or supply the path to your own EEMS run.
#' mcmcpath <- system.file("extdata", "EEMS-barrier", package = "reemsplots2")
#'
#' # Generate contour plots of migration and diversity rates
#' # as well as several diagnostic plots
#' plots <- make_eems_plots(mcmcpath, longlat = TRUE)
#' names(plots)
#' \donttest{
#' # Save the various plots
#' library("ggplot2")
#' plotpath <- file.path(path.expand("~"), "EEMS-barrier")
#' ggsave(paste0(plotpath, "-mrates01.png"), plots$mrates01, dpi = 600,
#' width = 6, height = 4)
#' ggsave(paste0(plotpath, "-mrates02.png"), plots$mrates02, dpi = 600,
#' width = 6, height = 4)
#' ggsave(paste0(plotpath, "-qrates01.png"), plots$qrates01, dpi = 600,
#' width = 6, height = 4)
#' ggsave(paste0(plotpath, "-qrates02.png"), plots$qrates02, dpi = 600,
#' width = 6, height = 4)
#' ggsave(paste0(plotpath, "-rdist01.pdf"), plots$rdist01,
#' width = 6.5, height = 6)
#' ggsave(paste0(plotpath, "-rdist02.pdf"), plots$rdist02,
#' width = 6.5, height = 6)
#' ggsave(paste0(plotpath, "-rdist03.pdf"), plots$rdist03,
#' width = 6.5, height = 6)
#' ggsave(paste0(plotpath, "-pilogl01.pdf"), plots$pilogl01,
#' width = 7, height = 5)
#' }
#' @seealso \code{\link{plot_population_grid}}, \code{\link{plot_resid_heatmap}}, \code{\link{plot_voronoi_tiles}}
#' @export
make_eems_plots <- function(mcmcpath, longlat = TRUE, dpi = 250,add_grid = FALSE, col_grid = "#BBBBBB",add_demes = FALSE, col_demes = "#000000",add_outline = FALSE, col_outline = "#FFFFFF",eems_colors = NULL, prob_level = 0.9,m_colscale = NULL, q_colscale = NULL,add_abline = FALSE) {
# Returns NULL if everything is okay
check_mcmcpath_contents(mcmcpath=mcmcpath)
func_params <- list(add_grid = add_grid, add_demes = add_demes, add_outline = add_outline,eems_colors = eems_colors, prob_level = prob_level,col_grid = col_grid, col_demes = col_demes,col_outline = col_outline,m_colscale = m_colscale, q_colscale = q_colscale,add_abline = add_abline)
# List of parameters defining how plots will be generated
plot_params <- check_plot_params(pars=func_params)
dimns <- read_dimns(mcmcpath=mcmcpath[1], longlat=longlat, nmrks = dpi)
plots <- list()
### Generate the first two of the four EEMS maps
p <- eems_contours(mcmcpath=mcmcpath, dimns=dimns, longlat=longlat, plot_params=plot_params, is_mrates = TRUE)
plots$mrates01 <- p[[1]]
plots$mrates02 <- p[[2]]
### Generate the third and fourth of the four EEMS maps
p <- eems_contours(mcmcpath=mcmcpath, dimns=dimns, longlat=longlat, plot_params=plot_params, is_mrates = FALSE)
plots$qrates01 <- p[[1]]
plots$qrates02 <- p[[2]]
### Generate the first three EEMS graphs (not maps)
# This will fail if fewer than two observed demes
obs_demes <- read_matrix(file.path(mcmcpath[1], "rdistoDemes.txt"), ncol = 3)
sizes <- obs_demes[, 3]
if (!sum(sizes > 1) < 2) {
dissimilarities <- pairwise_dist(mcmcpath=mcmcpath, longlat=longlat, plot_params=plot_params)
p <- plot_pairwise_dissimilarities_(dissimilarities=dissimilarities, add_abline=add_abline)
plots$rdist01 <- p[[1]]
plots$rdist02 <- p[[2]]
plots$rdist03 <- p[[3]]
}
### Generate the last EEMS graph (not a map)
# This will succeed even with fewer than two observed demes
plots$pilogl01 <- plot_log_posterior(mcmcpath=mcmcpath)
plots
}
eems_contours <- function(mcmcpath, dimns, longlat, plot_params, is_mrates) {
if (is_mrates){
message("Generate effective migration surface ","(posterior mean of m rates). ","See plots$mrates01 and plots$mrates02.")
}
else{
message("Generate effective diversity surface ","(posterior mean of q rates). ","See plots$qrates01 and plots$qrates02.")
}
zrates <- rep(0, dimns$nmrks)
pr_gt0 <- rep(0, dimns$nmrks)
pr_lt0 <- rep(0, dimns$nmrks)
niters <- 0
# Loop over each directory in mcmcpath to average the contour plots
for (path in mcmcpath) {
voronoi <- read_voronoi(mcmcpath=path, longlat=longlat, is_mrates=is_mrates)
# The function tiles2contours is coded in c++
rslt <- tiles2contours(tiles=voronoi$tiles, rates=voronoi$rates, seeds=cbind(voronoi$xseed, voronoi$yseed), marks=dimns$marks, distm=dimns$dist_metric)
zrates <- zrates + rslt$zrates
niters <- niters + rslt$niters
pr_gt0 <- pr_gt0 + rslt$pr_gt0
pr_lt0 <- pr_lt0 + rslt$pr_lt0
}
zrates <- zrates / niters
pr_gt0 <- pr_gt0 / niters
pr_lt0 <- pr_lt0 / niters
p1 <- filled_eems_contour(dimns=dimns, zmean=zrates, plot_params=plot_params, is_mrates=is_mrates)
p2 <- filled_prob_contour(dimns=dimns, probs=(pr_gt0 - pr_lt0), plot_params=plot_params, is_mrates=is_mrates)
result=list(p1, p2)
result
}
plot_log_posterior <- function(mcmcpath) {
message("Generate posterior probability trace. ","See plots$pilog01.")
rleid <- function(x) {
r <- rle(x)
rep(seq_along(r$lengths), r$lengths)
}
pl_df <- NULL
for (path in mcmcpath) {
pl <- read_matrix(file.path(path, "mcmcpilogl.txt"))
pl_df <- dplyr::bind_rows(pl_df, dplyr::as_data_frame(pl) %>% dplyr::mutate(path))
}
pl_df <- pl_df %>% setNames(c("pi", "logl", "path")) %>% dplyr::mutate(mcmcpath = factor(data.table::rleid(path))) %>% dplyr::group_by(mcmcpath) %>% dplyr::mutate(iter = dplyr::row_number(), pilogl = pi + logl)
ggplot2::ggplot(pl_df, ggplot2::aes(x = iter, y = pilogl, color = mcmcpath)) + ggplot2::geom_path() + ggplot2::labs(x = "MCMC iteration (after burn-in and thinning)",y = "log posterior",title = "Have the MCMC chains converged?",subtitle = "If not, restart EEMS and/or increase numMCMCIter") + ggplot2::theme_bw() + ggplot2::theme(panel.grid.minor = ggplot2::element_blank(),panel.grid.major.x = ggplot2::element_blank())
}
plot_pairwise_dissimilarities_ <- function(dissimilarities, add_abline) {
message("Generate average dissimilarities within and between demes. ","See plots$rdist01, plots$rdist02 and plots$rdist03.")
p1 <- ggplot2::ggplot(dissimilarities$between %>% filter(size > 1),ggplot2::aes(fitted, obsrvd)) + ggplot2::geom_point(shape = 1) + ggplot2::theme_minimal() + ggplot2::labs(x = expression(paste("Fitted dissimilarity between demes ", Delta[alpha * beta], " - (", Delta[alpha * alpha], " + ", Delta[beta * beta], ") / 2")), y = expression(paste("Observed dissimilarity between demes ",D[alpha * beta], " - (",D[alpha * alpha], " + ",D[beta * beta], ") / 2")), title = expression(paste("Dissimilarities between pairs of ", "sampled demes (", alpha, ", ", beta, ")")),subtitle = paste("Singleton demes, if any, are excluded from this", "plot (but not from EEMS)"))
p2 <- ggplot2::ggplot(dissimilarities$within %>% filter(size > 1), ggplot2::aes(fitted, obsrvd)) + ggplot2::geom_point(shape = 1) + ggplot2::theme_minimal() + ggplot2::labs(x = expression(paste("Fitted dissimilarity within demes ", Delta[alpha * alpha])), y = expression(paste("Observed dissimilarity within demes ",D[alpha * alpha])), title = expression(paste("Dissimilarities within sampled ","demes ", alpha)), subtitle = paste("Singleton demes, if any, are excluded from ","this plot (but not from EEMS)"))
p3 <- ggplot2::ggplot(dissimilarities$ibd %>% filter(size > 1), ggplot2::aes(fitted, obsrvd)) + ggplot2::geom_point(shape = 1) + ggplot2::theme_minimal() + ggplot2::labs(x = "Great circle distance between demes (km)",y = expression(paste("Observed dissimilarity between demes ",D[alpha * beta], " - (",D[alpha * alpha], " + ",D[beta * beta], ") / 2")),title = expression(paste("Dissimilarities between pairs of ","sampled demes (", alpha, ", ", beta, ")")),subtitle = paste("Singleton demes, if any, are excluded from this","plot (but not from EEMS)"))
if (add_abline) {
p1 <- p1 + ggplot2::geom_smooth(method = "lm", se = FALSE)
p2 <- p2 + ggplot2::geom_smooth(method = "lm", se = FALSE)
}
result=list(p1, p2, p3)
result
}
check_mcmcpath_contents <- function(mcmcpath) {
for (path in mcmcpath) {
for (file in c("rdistJtDobsJ.txt", "rdistJtDhatJ.txt", "rdistoDemes.txt",
"mcmcmtiles.txt", "mcmcmrates.txt", "mcmcxcoord.txt",
"mcmcycoord.txt", "mcmcqtiles.txt", "mcmcqrates.txt",
"mcmcwcoord.txt", "mcmczcoord.txt", "mcmcpilogl.txt",
"outer.txt", "demes.txt", "edges.txt", "ipmap.txt",
"eemsrun.txt")) {
if (!file.exists(file.path(path, file)))
stop("Each EEMS output folder should include ", file)
}
}
}
read_vector <- function(file) {
stopifnot(file.exists(file))
scan(file, what = numeric(), quiet = TRUE)
}
read_matrix <- function(file, ncol = 2) {
stopifnot(file.exists(file))
matrix(scan(file, what = numeric(), quiet = TRUE),
ncol = ncol, byrow = TRUE)
}
read_dimns <- function(mcmcpath, longlat, nmrks = 100) {
outer <- read_matrix(file.path(mcmcpath, "outer.txt"))
ipmap <- read_vector(file.path(mcmcpath, "ipmap.txt"))
demes <- read_matrix(file.path(mcmcpath, "demes.txt"))
edges <- read_matrix(file.path(mcmcpath, "edges.txt"))
dist_metric <- get_dist_metric(mcmcpath)
if (!longlat) {
outer <- outer[, c(2, 1)]
demes <- demes[, c(2, 1)]
}
xlim <- range(outer[, 1])
ylim <- range(outer[, 2])
aspect <- (diff(ylim) / diff(xlim)) / cos(mean(ylim) * pi / 180)
aspect <- abs(aspect)
if (aspect > 1) {
nxmrks <- nmrks
nymrks <- round(nxmrks * aspect)
} else {
nymrks <- nmrks
nxmrks <- round(nymrks / aspect)
}
# Construct a rectangular "raster" of equally spaced pixels/marks
xmrks <- seq(from = xlim[1], to = xlim[2], length = nxmrks)
ymrks <- seq(from = ylim[1], to = ylim[2], length = nymrks)
marks <- cbind(rep(xmrks, times = nymrks), rep(ymrks, each = nxmrks))
# Exclude pixels that fall outside the habitat outline
outer_poly <-sp::SpatialPolygons(list(Polygons(list(Polygon(outer, hole = FALSE)), "1")))
marks <- sp::SpatialPoints(marks)[outer_poly, ]
marks <- marks@coords
### set column names to x and y
outer <- dplyr::as_data_frame(outer) %>% setNames(c("x", "y"))
###
ipmap <- dplyr::data_frame(id = ipmap) %>% dplyr::count(id)
### set column names to x and y
demes <- dplyr::as_data_frame(demes) %>% setNames(c("x", "y")) %>% dplyr::mutate(id = dplyr::row_number()) %>% dplyr::left_join(ipmap) %>% dplyr::arrange(id) %>% dplyr::mutate(n = dplyr::if_else(is.na(n), 0L, n))
edges <- dplyr::bind_cols(demes[edges[, 1], ] %>% dplyr::select(x, y),demes[edges[, 2], ] %>% dplyr::select(x, y)) %>% setNames(c("x", "y", "xend", "yend"))
list(marks = marks, nmrks = nrow(marks), xlim = xlim, ylim = ylim,outer = outer, demes = demes, edges = edges, dist_metric = dist_metric)
}
decompose_distances <- function(diffs, sizes = NULL) {
# Diffs can have NAs on the main diagonal; these elements correspond to demes
# with a single observation. For such deme a, no dissimilarities between
# two distinct individuals are observed. I approximate diffs(a,a) with the
# average diffs(b,b) computed across demes b with multiple samples.
if (!is.null(sizes))
diag(diffs)[sizes < 2] <- mean(diag(diffs)[sizes >= 2])
within <- diag(diffs)
selfsim <- matrix(within, nrow(diffs), ncol(diffs))
between <- diffs - (selfsim + t(selfsim)) / 2
between <- between[upper.tri(between, diag = FALSE)]
list(within = within, between = between)
}
check_plot_params <- function(pars) {
if (is.logical(pars$add_grid)) {
pars$add_grid <- pars$add_grid[1]
} else {
pars$add_grid <- FALSE
}
if (is_color(pars$col_grid)) pars$col_grid <- pars$col_grid[1]
else pars$col_grid <- "#BBBBBB"
if (is.logical(pars$add_outline)) pars$add_outline <- pars$add_outline[1]
else pars$add_outline <- FALSE
if (is_color(pars$col_outline)) pars$col_outline <- pars$col_outline[1]
else pars$col_outline <- "#EEEEEE"
if (is.logical(pars$add_demes)) pars$add_demes <- pars$add_demes[1]
else pars$add_demes <- FALSE
if (is_color(pars$col_demes)) pars$col_demes <- pars$col_demes[1]
else pars$col_demes <- "#000000"
if (is.logical(pars$add_seeds)) pars$add_seeds <- pars$add_seeds[1]
else pars$add_seeds <- TRUE
if (is.numeric(pars$m_colscale)) pars$m_colscale <- pars$m_colscale
else pars$m_colscale <- c(-2.5, 2.5)
if (is.numeric(pars$q_colscale)) pars$q_colscale <- pars$q_colscale
else pars$q_colscale <- c(-0.1, 0.1)
if (length(pars$eems_colors) < 2 || any(!is_color(pars$eems_colors)))
pars$eems_colors <- default_eems_colors()
if (is.null(pars$prob_level)) prob_level <- 0.9
else prob_level <- pars$prob_level
prob_level <- prob_level[prob_level > 0.5 & prob_level < 1]
if (length(prob_level) != 1) prob_level <- 0.9
pars$prob_level <- prob_level
pars
}
geo_distm <- function(coord, longlat, plot_params) {
if (!longlat) coord <- coord[, c(2, 1)]
dist <- sp::spDists(coord, longlat = TRUE)
dist <- dist[upper.tri(dist, diag = FALSE)]
dist
}
pairwise_dist <- function(mcmcpath, longlat, plot_params) {
# List of observed demes, with number of samples taken collected Each row
# specifies: x coordinate, y coordinate, n samples
obs_demes <- read_matrix(file.path(mcmcpath[1], "rdistoDemes.txt"), ncol = 3)
sizes <- obs_demes[, 3]
if (sum(sizes > 1) < 2) {
message("There should be at least two observed demes ",
"to plot pairwise dissimilarities")
return(NULL)
}
npops <- nrow(obs_demes)
demes <- seq(npops)
diffs_obs <- matrix(0, npops, npops)
diffs_hat <- matrix(0, npops, npops)
for (path in mcmcpath) {
tempi <- read_matrix(file.path(path, "rdistoDemes.txt"), ncol = 3)
if (sum(dim(obs_demes) != dim(tempi)) || sum(obs_demes != tempi)) {
message("EEMS results for at least two different population grids. ",
"Plot pairwise dissimilarity for each grid separately.")
return(list(between = data_frame(), within = data_frame(),
ibd = data_frame()))
}
diffs_obs <- diffs_obs +
as.matrix(read.table(file.path(path, "rdistJtDobsJ.txt")))
diffs_hat <- diffs_hat +
as.matrix(read.table(file.path(path, "rdistJtDhatJ.txt")))
}
diffs_obs <- diffs_obs / length(mcmcpath)
diffs_hat <- diffs_hat / length(mcmcpath)
alpha <- matrix(demes, nrow = npops, ncol = npops)
beta <- t(alpha)
tempi <- matrix(sizes, npops, npops)
smaller_deme <- pmin(tempi, t(tempi))
smaller_deme <- smaller_deme[upper.tri(smaller_deme, diag = FALSE)]
alpha <- alpha[upper.tri(alpha, diag = FALSE)]
beta <- beta[upper.tri(beta, diag = FALSE)]
# Under pure isolation by distance, we expect the genetic dissimilarities
# between demes increase with the geographic distance separating them
dist <- geo_distm(obs_demes[, 1:2], longlat, plot_params)
bw_obs <- decompose_distances(diffs_obs, sizes)
bw_hat <- decompose_distances(diffs_hat)
b_component <- data_frame(alpha_x = obs_demes[, 1][alpha],
alpha_y = obs_demes[, 2][alpha],
beta_x = obs_demes[, 1][beta],
beta_y = obs_demes[, 2][beta],
fitted = bw_hat$between,
obsrvd = bw_obs$between,
size = smaller_deme)
w_component <- data_frame(alpha_x = obs_demes[, 1][demes],
alpha_y = obs_demes[, 2][demes],
fitted = bw_hat$within,
obsrvd = bw_obs$within,
size = sizes)
g_component <- data_frame(alpha_x = obs_demes[, 1][alpha],
alpha_y = obs_demes[, 2][alpha],
beta_x = obs_demes[, 1][beta],
beta_y = obs_demes[, 2][beta],
fitted = dist,
obsrvd = bw_obs$between,
size = smaller_deme)
list(between = b_component, within = w_component, ibd = g_component)
}
#######
## Add these functions to plot_eems.R
load_required_packages <- function(packages) {
for (package in packages) {
if (!requireNamespace(package, quietly = TRUE))
stop("The ", package, " package is required. ",
"Please install it first.")
else
message("Loading ", package, ".")
}
}
tiles2contours <- function(tiles, rates, seeds, marks, distm) {
.Call('_reemsplots2_tiles2contours', PACKAGE = 'reemsplots2', tiles, rates, seeds, marks, distm)
}
read_voronoi <- function(mcmcpath, longlat, is_mrates) {
if (is_mrates) {
rates <- read_vector(file.path(mcmcpath, "mcmcmrates.txt"))
tiles <- read_vector(file.path(mcmcpath, "mcmcmtiles.txt"))
xseed <- read_vector(file.path(mcmcpath, "mcmcxcoord.txt"))
yseed <- read_vector(file.path(mcmcpath, "mcmcycoord.txt"))
} else {
rates <- read_vector(file.path(mcmcpath, "mcmcqrates.txt"))
tiles <- read_vector(file.path(mcmcpath, "mcmcqtiles.txt"))
xseed <- read_vector(file.path(mcmcpath, "mcmcwcoord.txt"))
yseed <- read_vector(file.path(mcmcpath, "mcmczcoord.txt"))
}
if (!longlat) {
tempi <- xseed
xseed <- yseed
yseed <- tempi
}
list(rates = log10(rates), tiles = tiles, xseed = xseed, yseed = yseed)
}
# Get the distance metric from `eemsrun.txt`; use `euclidean` by default.
get_dist_metric <- function(mcmcpath) {
dist_metric <- "euclidean"
lines <- tolower(readLines(file.path(mcmcpath, "eemsrun.txt")))
for (line in lines) {
if (grepl("\\s*distance\\s*=\\s*", line))
dist_metric <- gsub("\\s*distance\\s*=\\s*(\\w+)", "\\1", line)
}
if (dist_metric != "euclidean" && dist_metric != "greatcirc")
stop("eemsrun.txt should specify `euclidean` or `greatcirc` distance.")
dist_metric
}
default_eems_colors <- function() {
# Use the default DarkOrange to Blue color scheme, which combines
# two color schemes from the `dichromat` package. These are based
# on a collection of color schemes for scientific graphics:
# See http://geog.uoregon.edu/datagraphics/color_scales.htm
# To reproduce the default eems colors,
# let oranges be dichromat::colorschemes$BluetoDarkOrange.12[12:7]
# and blues be dichromat::colorschemes$BrowntoBlue.12[7:12]
c("#994000", "#CC5800", "#FF8F33", "#FFAD66", "#FFCA99", "#FFE6CC",
"#FBFBFB",
"#CCFDFF", "#99F8FF", "#66F0FF", "#33E4FF", "#00AACC", "#007A99")
}
# Check that a string, or a vector of strings, represents a hex color
is_color <- function(x) {
if (is.null(x)) return(FALSE)
sapply(x, function(x) {
tryCatch(is.matrix(col2rgb(x)), error = function(e) FALSE)
})
}
theme_void <- function() {
ggplot2::theme(line = ggplot2::element_blank(), rect = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(), axis.title = ggplot2::element_blank(),
legend.text = ggplot2::element_text(size = ggplot2::rel(0.8)),
legend.title = ggplot2::element_text(hjust = 0),
legend.text.align = 1)
}
filled_contour_rates <- function(z, dimns) {
w <- cbind(dimns$marks, z)
colnames(w) <- c("x", "y", "z")
ggplot2::ggplot(dplyr::as_data_frame(w), ggplot2::aes(x = x, y = y)) +
ggplot2::geom_tile(ggplot2::aes(fill = z)) +
ggplot2::scale_x_continuous(limits = dimns$xlim) +
ggplot2::scale_y_continuous(limits = dimns$ylim) +
ggplot2::coord_quickmap() +
theme_void() +
ggplot2::theme(legend.text.align = 1)
}
filled_contour_graph <- function(p, dimns, plot_params) {
if (plot_params$add_grid) {
p <- p + ggplot2::geom_segment(data = dimns$edges,
ggplot2::aes(x = x, y = y, xend = xend, yend = yend),
color = plot_params$col_grid)
}
if (plot_params$add_demes) {
p <- p + ggplot2::geom_point(data = dimns$demes %>% dplyr::filter(n > 0), ggplot2::aes(x = x, y = y, size = n), shape = 1,color = plot_params$col_demes) +
ggplot2::scale_size_continuous(guide = FALSE)
}
if (plot_params$add_outline) {
p <- p + ggplot2::geom_path(data = dimns$outer, ggplot2::aes(x = x, y = y),
color = plot_params$col_outline)
}
p
}
filled_eems_contour <- function(dimns, zmean, plot_params, is_mrates) {
if (is_mrates) {
title <- "log(m)"
colscale <- plot_params$m_colscale
} else {
title <- "log(q)"
colscale <- plot_params$q_colscale
}
limits <- range(zmean, colscale, na.rm = TRUE, finite = TRUE)
p <- filled_contour_rates(zmean, dimns)
p <- filled_contour_graph(p, dimns, plot_params) +
ggplot2::scale_fill_gradientn(colors = plot_params$eems_colors,
limits = limits, name = title)
p
}
filled_prob_contour <- function(dimns, probs, plot_params, is_mrates) {
probs <- (probs + 1) / 2
probs[probs < 0] <- 0
probs[probs > 1] <- 1
if (is_mrates) r <- "m" else r <- "q"
breaks <- c(1 - plot_params$prob_level, plot_params$prob_level)
labels <- c(paste0("P{log(", r, ") < 0} = ", plot_params$prob_level),
paste0("P{log(", r, ") > 0} = ", plot_params$prob_level))
p <- filled_contour_rates(probs, dimns)
p <- filled_contour_graph(p, dimns, plot_params) +
ggplot2::geom_contour(ggplot2::aes(z = z), breaks = breaks, color = "white") +
ggplot2::scale_fill_gradientn(colors = default_eems_colors(),
limits = c(0, 1), name = "",
breaks = breaks, labels = labels)
p
}
JeffWeinell/misc.wrappers documentation built on Sept. 20, 2023, 12:42 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#####
## All of these functions are from the R package 'reemsplots2' by Desislava Petkova (github repository 'dipetkov/reemsplots2').
## Only very minor edits were been made by me (Jeff Weinell), for a few of the functions, but in most cases the functions are identical to those in reemsplots2.
#' Plot effective migration and diversity surfaces
#'
#' Given a vector of EEMS output directories, this function generates several figures to visualize EEMS results. It is a good idea to examine all these figures, which is why they are generated by default.
#' \describe{
#' \item{mrates01}{effective migration surface. This contour plot visualizes the estimated effective migration rates \code{m}, on the log10 scale after mean centering.}
#' \item{mrates02}{posterior probability contours \code{P(log(m) > 0) = p} and \code{P(log(m) < 0) = p} for the given probability level \code{p}. Since migration rates are visualized on the log10 scale after mean centering, 0 corresponds to the overall mean migration rate. This contour plot emphasizes regions with effective migration that is significantly higher/lower than the overall average.}
#' \item{qrates01}{effective diversity surface. This contour plot visualizes the estimated effective diversity rates \code{q}, on the log10 scale after mean centering.}
#' \item{qrates02}{posterior probability contours \code{P(log(q) > 0) = p} and \code{P(log(q) < 0) = p}. Similar to \code{mrates02} but applied to the effective diversity rates.}
#' \item{rdist01}{scatter plot of the observed vs the fitted between-deme component of genetic dissimilarity, where one point represents a pair of sampled demes.}
#' \item{rdist01}{scatter plot of the observed vs the fitted within-deme component of genetic dissimilarity, where one point represents a sampled deme.}
#' \item{rdist03}{scatter plot of observed genetic dissimilarities between demes vs observed geographic distances between demes.}
#' \item{pilogl01}{posterior probability trace}
#' }
#' @details
#' The function \code{make_eems_plots} will work given the results from a single EEMS run (one directory in \code{mcmcpath}) but it is better to run EEMS several times, randomly initializing the MCMC chain each time. In other words, simulate several realizations of the Markov chain and let each realization start from a different state in the parameter space (by using a different random seed).
##'
#' The \code{mrates} and \code{qrates} figures visualize (properties of) the effective migration and diversity rates across the habitat. The other figures can help to check that the MCMC sampler has converged (the trace plot \code{pilogl}) and that the EEMS model fits the data well (the scatter plots of genetic dissimilarities \code{rdist}).
#'
#' To describe the within-deme and between-deme components of genetic dissimilarity, let \code{D(a,b)} be the dissimilarity between one individual from deme \code{a} and another individual from deme \code{b}. Then the within-deme component for \code{a} and \code{b} is simply \code{D(a,a)} and \code{D(b, b)}, respectively. The between-deme component is \code{D(a,b) - [D(a,a) + D(b,b)] / 2} and it represents dissimilarity that is due to the spatial structure of the population and is not a consequence of the local diversity in the two demes.
#' @param mcmcpath A vector of EEMS output directories, for the same dataset. Warning: There is minimal checking that the directories all correspond to the same dataset.
#' @param longlat A logical value indicating whether the coordinates are given as pairs (longitude, latitude) or (latitude, longitude).
#' @param dpi Resolution of the contour raster. The default is 250.
#' @param add_grid A logical value indicating whether to add the population grid or not.
#' @param col_grid The color of the population grid. Defaults to gray.
#' @param add_outline A logical value indicating whether to add the habitat outline or not.
#' @param col_outline The color of the habitat outline. Defaults to white.
#' @param add_demes A logical value indicating whether to add the observed demes or not.
#' @param col_demes The color of the demes. Defaults to black.
#' @param eems_colors The EEMS color scheme as a vector of colors, ordered from low to high. Defaults to a DarkOrange to Blue divergent palette with six orange shades, white in the middle, six blue shades. Acknowledgement: The default color scheme is adapted from the \code{dichromat} package.
#' @param m_colscale A fixed range for log10-transformed migration rates. If the estimated rates fall outside the specified range, then the color scale is ignored. The default range is \code{[-2.5, +2.5]}.
#' @param q_colscale A fixed range for log10-transformed diversity rates. The default range is \code{-0.1, +0.1}.
#' @param add_abline Add the line \code{y = x} to the diagnostic scatter plots of observed vs fitted genetic dissimilarities.
#' @param prob_level A probability \code{p} to define the posterior probability contours \code{P(log(m) > 0) = p} and \code{P(log(m) < 0) = p}. Defaults to \code{0.9}.
#' @references Light A and Bartlein PJ (2004). The End of the Rainbow? Color Schemes for Improved Data Graphics. EOS Transactions of the American Geophysical Union, 85(40), 385.
#' @examples
#' # Use the provided example or supply the path to your own EEMS run.
#' mcmcpath <- system.file("extdata", "EEMS-barrier", package = "reemsplots2")
#'
#' # Generate contour plots of migration and diversity rates
#' # as well as several diagnostic plots
#' plots <- make_eems_plots(mcmcpath, longlat = TRUE)
#' names(plots)
#' \donttest{
#' # Save the various plots
#' library("ggplot2")
#' plotpath <- file.path(path.expand("~"), "EEMS-barrier")
#' ggsave(paste0(plotpath, "-mrates01.png"), plots$mrates01, dpi = 600,
#' width = 6, height = 4)
#' ggsave(paste0(plotpath, "-mrates02.png"), plots$mrates02, dpi = 600,
#' width = 6, height = 4)
#' ggsave(paste0(plotpath, "-qrates01.png"), plots$qrates01, dpi = 600,
#' width = 6, height = 4)
#' ggsave(paste0(plotpath, "-qrates02.png"), plots$qrates02, dpi = 600,
#' width = 6, height = 4)
#' ggsave(paste0(plotpath, "-rdist01.pdf"), plots$rdist01,
#' width = 6.5, height = 6)
#' ggsave(paste0(plotpath, "-rdist02.pdf"), plots$rdist02,
#' width = 6.5, height = 6)
#' ggsave(paste0(plotpath, "-rdist03.pdf"), plots$rdist03,
#' width = 6.5, height = 6)
#' ggsave(paste0(plotpath, "-pilogl01.pdf"), plots$pilogl01,
#' width = 7, height = 5)
#' }
#' @seealso \code{\link{plot_population_grid}}, \code{\link{plot_resid_heatmap}}, \code{\link{plot_voronoi_tiles}}
#' @export
make_eems_plots <- function(mcmcpath, longlat = TRUE, dpi = 250,add_grid = FALSE, col_grid = "#BBBBBB",add_demes = FALSE, col_demes = "#000000",add_outline = FALSE, col_outline = "#FFFFFF",eems_colors = NULL, prob_level = 0.9,m_colscale = NULL, q_colscale = NULL,add_abline = FALSE) {
# Returns NULL if everything is okay
check_mcmcpath_contents(mcmcpath=mcmcpath)
func_params <- list(add_grid = add_grid, add_demes = add_demes, add_outline = add_outline,eems_colors = eems_colors, prob_level = prob_level,col_grid = col_grid, col_demes = col_demes,col_outline = col_outline,m_colscale = m_colscale, q_colscale = q_colscale,add_abline = add_abline)
# List of parameters defining how plots will be generated
plot_params <- check_plot_params(pars=func_params)
dimns <- read_dimns(mcmcpath=mcmcpath[1], longlat=longlat, nmrks = dpi)
plots <- list()
### Generate the first two of the four EEMS maps
p <- eems_contours(mcmcpath=mcmcpath, dimns=dimns, longlat=longlat, plot_params=plot_params, is_mrates = TRUE)
plots$mrates01 <- p[[1]]
plots$mrates02 <- p[[2]]
### Generate the third and fourth of the four EEMS maps
p <- eems_contours(mcmcpath=mcmcpath, dimns=dimns, longlat=longlat, plot_params=plot_params, is_mrates = FALSE)
plots$qrates01 <- p[[1]]
plots$qrates02 <- p[[2]]
### Generate the first three EEMS graphs (not maps)
# This will fail if fewer than two observed demes
obs_demes <- read_matrix(file.path(mcmcpath[1], "rdistoDemes.txt"), ncol = 3)
sizes <- obs_demes[, 3]
if (!sum(sizes > 1) < 2) {
dissimilarities <- pairwise_dist(mcmcpath=mcmcpath, longlat=longlat, plot_params=plot_params)
p <- plot_pairwise_dissimilarities_(dissimilarities=dissimilarities, add_abline=add_abline)
plots$rdist01 <- p[[1]]
plots$rdist02 <- p[[2]]
plots$rdist03 <- p[[3]]
}
### Generate the last EEMS graph (not a map)
# This will succeed even with fewer than two observed demes
plots$pilogl01 <- plot_log_posterior(mcmcpath=mcmcpath)
plots
}
eems_contours <- function(mcmcpath, dimns, longlat, plot_params, is_mrates) {
if (is_mrates){
message("Generate effective migration surface ","(posterior mean of m rates). ","See plots$mrates01 and plots$mrates02.")
}
else{
message("Generate effective diversity surface ","(posterior mean of q rates). ","See plots$qrates01 and plots$qrates02.")
}
zrates <- rep(0, dimns$nmrks)
pr_gt0 <- rep(0, dimns$nmrks)
pr_lt0 <- rep(0, dimns$nmrks)
niters <- 0
# Loop over each directory in mcmcpath to average the contour plots
for (path in mcmcpath) {
voronoi <- read_voronoi(mcmcpath=path, longlat=longlat, is_mrates=is_mrates)
# The function tiles2contours is coded in c++
rslt <- tiles2contours(tiles=voronoi$tiles, rates=voronoi$rates, seeds=cbind(voronoi$xseed, voronoi$yseed), marks=dimns$marks, distm=dimns$dist_metric)
zrates <- zrates + rslt$zrates
niters <- niters + rslt$niters
pr_gt0 <- pr_gt0 + rslt$pr_gt0
pr_lt0 <- pr_lt0 + rslt$pr_lt0
}
zrates <- zrates / niters
pr_gt0 <- pr_gt0 / niters
pr_lt0 <- pr_lt0 / niters
p1 <- filled_eems_contour(dimns=dimns, zmean=zrates, plot_params=plot_params, is_mrates=is_mrates)
p2 <- filled_prob_contour(dimns=dimns, probs=(pr_gt0 - pr_lt0), plot_params=plot_params, is_mrates=is_mrates)
result=list(p1, p2)
result
}
plot_log_posterior <- function(mcmcpath) {
message("Generate posterior probability trace. ","See plots$pilog01.")
rleid <- function(x) {
r <- rle(x)
rep(seq_along(r$lengths), r$lengths)
}
pl_df <- NULL
for (path in mcmcpath) {
pl <- read_matrix(file.path(path, "mcmcpilogl.txt"))
pl_df <- dplyr::bind_rows(pl_df, dplyr::as_data_frame(pl) %>% dplyr::mutate(path))
}
pl_df <- pl_df %>% setNames(c("pi", "logl", "path")) %>% dplyr::mutate(mcmcpath = factor(data.table::rleid(path))) %>% dplyr::group_by(mcmcpath) %>% dplyr::mutate(iter = dplyr::row_number(), pilogl = pi + logl)
ggplot2::ggplot(pl_df, ggplot2::aes(x = iter, y = pilogl, color = mcmcpath)) + ggplot2::geom_path() + ggplot2::labs(x = "MCMC iteration (after burn-in and thinning)",y = "log posterior",title = "Have the MCMC chains converged?",subtitle = "If not, restart EEMS and/or increase numMCMCIter") + ggplot2::theme_bw() + ggplot2::theme(panel.grid.minor = ggplot2::element_blank(),panel.grid.major.x = ggplot2::element_blank())
}
plot_pairwise_dissimilarities_ <- function(dissimilarities, add_abline) {
message("Generate average dissimilarities within and between demes. ","See plots$rdist01, plots$rdist02 and plots$rdist03.")
p1 <- ggplot2::ggplot(dissimilarities$between %>% filter(size > 1),ggplot2::aes(fitted, obsrvd)) + ggplot2::geom_point(shape = 1) + ggplot2::theme_minimal() + ggplot2::labs(x = expression(paste("Fitted dissimilarity between demes ", Delta[alpha * beta], " - (", Delta[alpha * alpha], " + ", Delta[beta * beta], ") / 2")), y = expression(paste("Observed dissimilarity between demes ",D[alpha * beta], " - (",D[alpha * alpha], " + ",D[beta * beta], ") / 2")), title = expression(paste("Dissimilarities between pairs of ", "sampled demes (", alpha, ", ", beta, ")")),subtitle = paste("Singleton demes, if any, are excluded from this", "plot (but not from EEMS)"))
p2 <- ggplot2::ggplot(dissimilarities$within %>% filter(size > 1), ggplot2::aes(fitted, obsrvd)) + ggplot2::geom_point(shape = 1) + ggplot2::theme_minimal() + ggplot2::labs(x = expression(paste("Fitted dissimilarity within demes ", Delta[alpha * alpha])), y = expression(paste("Observed dissimilarity within demes ",D[alpha * alpha])), title = expression(paste("Dissimilarities within sampled ","demes ", alpha)), subtitle = paste("Singleton demes, if any, are excluded from ","this plot (but not from EEMS)"))
p3 <- ggplot2::ggplot(dissimilarities$ibd %>% filter(size > 1), ggplot2::aes(fitted, obsrvd)) + ggplot2::geom_point(shape = 1) + ggplot2::theme_minimal() + ggplot2::labs(x = "Great circle distance between demes (km)",y = expression(paste("Observed dissimilarity between demes ",D[alpha * beta], " - (",D[alpha * alpha], " + ",D[beta * beta], ") / 2")),title = expression(paste("Dissimilarities between pairs of ","sampled demes (", alpha, ", ", beta, ")")),subtitle = paste("Singleton demes, if any, are excluded from this","plot (but not from EEMS)"))
if (add_abline) {
p1 <- p1 + ggplot2::geom_smooth(method = "lm", se = FALSE)
p2 <- p2 + ggplot2::geom_smooth(method = "lm", se = FALSE)
}
result=list(p1, p2, p3)
result
}
check_mcmcpath_contents <- function(mcmcpath) {
for (path in mcmcpath) {
for (file in c("rdistJtDobsJ.txt", "rdistJtDhatJ.txt", "rdistoDemes.txt",
"mcmcmtiles.txt", "mcmcmrates.txt", "mcmcxcoord.txt",
"mcmcycoord.txt", "mcmcqtiles.txt", "mcmcqrates.txt",
"mcmcwcoord.txt", "mcmczcoord.txt", "mcmcpilogl.txt",
"outer.txt", "demes.txt", "edges.txt", "ipmap.txt",
"eemsrun.txt")) {
if (!file.exists(file.path(path, file)))
stop("Each EEMS output folder should include ", file)
}
}
}
read_vector <- function(file) {
stopifnot(file.exists(file))
scan(file, what = numeric(), quiet = TRUE)
}
read_matrix <- function(file, ncol = 2) {
stopifnot(file.exists(file))
matrix(scan(file, what = numeric(), quiet = TRUE),
ncol = ncol, byrow = TRUE)
}
read_dimns <- function(mcmcpath, longlat, nmrks = 100) {
outer <- read_matrix(file.path(mcmcpath, "outer.txt"))
ipmap <- read_vector(file.path(mcmcpath, "ipmap.txt"))
demes <- read_matrix(file.path(mcmcpath, "demes.txt"))
edges <- read_matrix(file.path(mcmcpath, "edges.txt"))
dist_metric <- get_dist_metric(mcmcpath)
if (!longlat) {
outer <- outer[, c(2, 1)]
demes <- demes[, c(2, 1)]
}
xlim <- range(outer[, 1])
ylim <- range(outer[, 2])
aspect <- (diff(ylim) / diff(xlim)) / cos(mean(ylim) * pi / 180)
aspect <- abs(aspect)
if (aspect > 1) {
nxmrks <- nmrks
nymrks <- round(nxmrks * aspect)
} else {
nymrks <- nmrks
nxmrks <- round(nymrks / aspect)
}
# Construct a rectangular "raster" of equally spaced pixels/marks
xmrks <- seq(from = xlim[1], to = xlim[2], length = nxmrks)
ymrks <- seq(from = ylim[1], to = ylim[2], length = nymrks)
marks <- cbind(rep(xmrks, times = nymrks), rep(ymrks, each = nxmrks))
# Exclude pixels that fall outside the habitat outline
outer_poly <-sp::SpatialPolygons(list(Polygons(list(Polygon(outer, hole = FALSE)), "1")))
marks <- sp::SpatialPoints(marks)[outer_poly, ]
marks <- marks@coords
### set column names to x and y
outer <- dplyr::as_data_frame(outer) %>% setNames(c("x", "y"))
###
ipmap <- dplyr::data_frame(id = ipmap) %>% dplyr::count(id)
### set column names to x and y
demes <- dplyr::as_data_frame(demes) %>% setNames(c("x", "y")) %>% dplyr::mutate(id = dplyr::row_number()) %>% dplyr::left_join(ipmap) %>% dplyr::arrange(id) %>% dplyr::mutate(n = dplyr::if_else(is.na(n), 0L, n))
edges <- dplyr::bind_cols(demes[edges[, 1], ] %>% dplyr::select(x, y),demes[edges[, 2], ] %>% dplyr::select(x, y)) %>% setNames(c("x", "y", "xend", "yend"))
list(marks = marks, nmrks = nrow(marks), xlim = xlim, ylim = ylim,outer = outer, demes = demes, edges = edges, dist_metric = dist_metric)
}
decompose_distances <- function(diffs, sizes = NULL) {
# Diffs can have NAs on the main diagonal; these elements correspond to demes
# with a single observation. For such deme a, no dissimilarities between
# two distinct individuals are observed. I approximate diffs(a,a) with the
# average diffs(b,b) computed across demes b with multiple samples.
if (!is.null(sizes))
diag(diffs)[sizes < 2] <- mean(diag(diffs)[sizes >= 2])
within <- diag(diffs)
selfsim <- matrix(within, nrow(diffs), ncol(diffs))
between <- diffs - (selfsim + t(selfsim)) / 2
between <- between[upper.tri(between, diag = FALSE)]
list(within = within, between = between)
}
check_plot_params <- function(pars) {
if (is.logical(pars$add_grid)) {
pars$add_grid <- pars$add_grid[1]
} else {
pars$add_grid <- FALSE
}
if (is_color(pars$col_grid)) pars$col_grid <- pars$col_grid[1]
else pars$col_grid <- "#BBBBBB"
if (is.logical(pars$add_outline)) pars$add_outline <- pars$add_outline[1]
else pars$add_outline <- FALSE
if (is_color(pars$col_outline)) pars$col_outline <- pars$col_outline[1]
else pars$col_outline <- "#EEEEEE"
if (is.logical(pars$add_demes)) pars$add_demes <- pars$add_demes[1]
else pars$add_demes <- FALSE
if (is_color(pars$col_demes)) pars$col_demes <- pars$col_demes[1]
else pars$col_demes <- "#000000"
if (is.logical(pars$add_seeds)) pars$add_seeds <- pars$add_seeds[1]
else pars$add_seeds <- TRUE
if (is.numeric(pars$m_colscale)) pars$m_colscale <- pars$m_colscale
else pars$m_colscale <- c(-2.5, 2.5)
if (is.numeric(pars$q_colscale)) pars$q_colscale <- pars$q_colscale
else pars$q_colscale <- c(-0.1, 0.1)
if (length(pars$eems_colors) < 2 || any(!is_color(pars$eems_colors)))
pars$eems_colors <- default_eems_colors()
if (is.null(pars$prob_level)) prob_level <- 0.9
else prob_level <- pars$prob_level
prob_level <- prob_level[prob_level > 0.5 & prob_level < 1]
if (length(prob_level) != 1) prob_level <- 0.9
pars$prob_level <- prob_level
pars
}
geo_distm <- function(coord, longlat, plot_params) {
if (!longlat) coord <- coord[, c(2, 1)]
dist <- sp::spDists(coord, longlat = TRUE)
dist <- dist[upper.tri(dist, diag = FALSE)]
dist
}
pairwise_dist <- function(mcmcpath, longlat, plot_params) {
# List of observed demes, with number of samples taken collected Each row
# specifies: x coordinate, y coordinate, n samples
obs_demes <- read_matrix(file.path(mcmcpath[1], "rdistoDemes.txt"), ncol = 3)
sizes <- obs_demes[, 3]
if (sum(sizes > 1) < 2) {
message("There should be at least two observed demes ",
"to plot pairwise dissimilarities")
return(NULL)
}
npops <- nrow(obs_demes)
demes <- seq(npops)
diffs_obs <- matrix(0, npops, npops)
diffs_hat <- matrix(0, npops, npops)
for (path in mcmcpath) {
tempi <- read_matrix(file.path(path, "rdistoDemes.txt"), ncol = 3)
if (sum(dim(obs_demes) != dim(tempi)) || sum(obs_demes != tempi)) {
message("EEMS results for at least two different population grids. ",
"Plot pairwise dissimilarity for each grid separately.")
return(list(between = data_frame(), within = data_frame(),
ibd = data_frame()))
}
diffs_obs <- diffs_obs +
as.matrix(read.table(file.path(path, "rdistJtDobsJ.txt")))
diffs_hat <- diffs_hat +
as.matrix(read.table(file.path(path, "rdistJtDhatJ.txt")))
}
diffs_obs <- diffs_obs / length(mcmcpath)
diffs_hat <- diffs_hat / length(mcmcpath)
alpha <- matrix(demes, nrow = npops, ncol = npops)
beta <- t(alpha)
tempi <- matrix(sizes, npops, npops)
smaller_deme <- pmin(tempi, t(tempi))
smaller_deme <- smaller_deme[upper.tri(smaller_deme, diag = FALSE)]
alpha <- alpha[upper.tri(alpha, diag = FALSE)]
beta <- beta[upper.tri(beta, diag = FALSE)]
# Under pure isolation by distance, we expect the genetic dissimilarities
# between demes increase with the geographic distance separating them
dist <- geo_distm(obs_demes[, 1:2], longlat, plot_params)
bw_obs <- decompose_distances(diffs_obs, sizes)
bw_hat <- decompose_distances(diffs_hat)
b_component <- data_frame(alpha_x = obs_demes[, 1][alpha],
alpha_y = obs_demes[, 2][alpha],
beta_x = obs_demes[, 1][beta],
beta_y = obs_demes[, 2][beta],
fitted = bw_hat$between,
obsrvd = bw_obs$between,
size = smaller_deme)
w_component <- data_frame(alpha_x = obs_demes[, 1][demes],
alpha_y = obs_demes[, 2][demes],
fitted = bw_hat$within,
obsrvd = bw_obs$within,
size = sizes)
g_component <- data_frame(alpha_x = obs_demes[, 1][alpha],
alpha_y = obs_demes[, 2][alpha],
beta_x = obs_demes[, 1][beta],
beta_y = obs_demes[, 2][beta],
fitted = dist,
obsrvd = bw_obs$between,
size = smaller_deme)
list(between = b_component, within = w_component, ibd = g_component)
}
#######
## Add these functions to plot_eems.R
load_required_packages <- function(packages) {
for (package in packages) {
if (!requireNamespace(package, quietly = TRUE))
stop("The ", package, " package is required. ",
"Please install it first.")
else
message("Loading ", package, ".")
}
}
tiles2contours <- function(tiles, rates, seeds, marks, distm) {
.Call('_reemsplots2_tiles2contours', PACKAGE = 'reemsplots2', tiles, rates, seeds, marks, distm)
}
read_voronoi <- function(mcmcpath, longlat, is_mrates) {
if (is_mrates) {
rates <- read_vector(file.path(mcmcpath, "mcmcmrates.txt"))
tiles <- read_vector(file.path(mcmcpath, "mcmcmtiles.txt"))
xseed <- read_vector(file.path(mcmcpath, "mcmcxcoord.txt"))
yseed <- read_vector(file.path(mcmcpath, "mcmcycoord.txt"))
} else {
rates <- read_vector(file.path(mcmcpath, "mcmcqrates.txt"))
tiles <- read_vector(file.path(mcmcpath, "mcmcqtiles.txt"))
xseed <- read_vector(file.path(mcmcpath, "mcmcwcoord.txt"))
yseed <- read_vector(file.path(mcmcpath, "mcmczcoord.txt"))
}
if (!longlat) {
tempi <- xseed
xseed <- yseed
yseed <- tempi
}
list(rates = log10(rates), tiles = tiles, xseed = xseed, yseed = yseed)
}
# Get the distance metric from `eemsrun.txt`; use `euclidean` by default.
get_dist_metric <- function(mcmcpath) {
dist_metric <- "euclidean"
lines <- tolower(readLines(file.path(mcmcpath, "eemsrun.txt")))
for (line in lines) {
if (grepl("\\s*distance\\s*=\\s*", line))
dist_metric <- gsub("\\s*distance\\s*=\\s*(\\w+)", "\\1", line)
}
if (dist_metric != "euclidean" && dist_metric != "greatcirc")
stop("eemsrun.txt should specify `euclidean` or `greatcirc` distance.")
dist_metric
}
default_eems_colors <- function() {
# Use the default DarkOrange to Blue color scheme, which combines
# two color schemes from the `dichromat` package. These are based
# on a collection of color schemes for scientific graphics:
# See http://geog.uoregon.edu/datagraphics/color_scales.htm
# To reproduce the default eems colors,
# let oranges be dichromat::colorschemes$BluetoDarkOrange.12[12:7]
# and blues be dichromat::colorschemes$BrowntoBlue.12[7:12]
c("#994000", "#CC5800", "#FF8F33", "#FFAD66", "#FFCA99", "#FFE6CC",
"#FBFBFB",
"#CCFDFF", "#99F8FF", "#66F0FF", "#33E4FF", "#00AACC", "#007A99")
}
# Check that a string, or a vector of strings, represents a hex color
is_color <- function(x) {
if (is.null(x)) return(FALSE)
sapply(x, function(x) {
tryCatch(is.matrix(col2rgb(x)), error = function(e) FALSE)
})
}
theme_void <- function() {
ggplot2::theme(line = ggplot2::element_blank(), rect = ggplot2::element_blank(),
axis.text = ggplot2::element_blank(), axis.title = ggplot2::element_blank(),
legend.text = ggplot2::element_text(size = ggplot2::rel(0.8)),
legend.title = ggplot2::element_text(hjust = 0),
legend.text.align = 1)
}
filled_contour_rates <- function(z, dimns) {
w <- cbind(dimns$marks, z)
colnames(w) <- c("x", "y", "z")
ggplot2::ggplot(dplyr::as_data_frame(w), ggplot2::aes(x = x, y = y)) +
ggplot2::geom_tile(ggplot2::aes(fill = z)) +
ggplot2::scale_x_continuous(limits = dimns$xlim) +
ggplot2::scale_y_continuous(limits = dimns$ylim) +
ggplot2::coord_quickmap() +
theme_void() +
ggplot2::theme(legend.text.align = 1)
}
filled_contour_graph <- function(p, dimns, plot_params) {
if (plot_params$add_grid) {
p <- p + ggplot2::geom_segment(data = dimns$edges,
ggplot2::aes(x = x, y = y, xend = xend, yend = yend),
color = plot_params$col_grid)
}
if (plot_params$add_demes) {
p <- p + ggplot2::geom_point(data = dimns$demes %>% dplyr::filter(n > 0), ggplot2::aes(x = x, y = y, size = n), shape = 1,color = plot_params$col_demes) +
ggplot2::scale_size_continuous(guide = FALSE)
}
if (plot_params$add_outline) {
p <- p + ggplot2::geom_path(data = dimns$outer, ggplot2::aes(x = x, y = y),
color = plot_params$col_outline)
}
p
}
filled_eems_contour <- function(dimns, zmean, plot_params, is_mrates) {
if (is_mrates) {
title <- "log(m)"
colscale <- plot_params$m_colscale
} else {
title <- "log(q)"
colscale <- plot_params$q_colscale
}
limits <- range(zmean, colscale, na.rm = TRUE, finite = TRUE)
p <- filled_contour_rates(zmean, dimns)
p <- filled_contour_graph(p, dimns, plot_params) +
ggplot2::scale_fill_gradientn(colors = plot_params$eems_colors,
limits = limits, name = title)
p
}
filled_prob_contour <- function(dimns, probs, plot_params, is_mrates) {
probs <- (probs + 1) / 2
probs[probs < 0] <- 0
probs[probs > 1] <- 1
if (is_mrates) r <- "m" else r <- "q"
breaks <- c(1 - plot_params$prob_level, plot_params$prob_level)
labels <- c(paste0("P{log(", r, ") < 0} = ", plot_params$prob_level),
paste0("P{log(", r, ") > 0} = ", plot_params$prob_level))
p <- filled_contour_rates(probs, dimns)
p <- filled_contour_graph(p, dimns, plot_params) +
ggplot2::geom_contour(ggplot2::aes(z = z), breaks = breaks, color = "white") +
ggplot2::scale_fill_gradientn(colors = default_eems_colors(),
limits = c(0, 1), name = "",
breaks = breaks, labels = labels)
p
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.