R/gl.map.structure.r

In dartR.popgen: Analysing 'SNP' and 'Silicodart' Data Generated by Genome-Wide Restriction Fragment Analysis

#' @name gl.map.structure
#' @title Maps a STRUCTURE plot using a genlight object
#' @description
#' This function takes the output of plotstructure (the q matrix) and maps the
#' q-matrix across using the population centers from the genlight object that
#' was used to run the structure analysis via \code{\link{gl.run.structure}})
#' and plots the typical structure bar plots on a spatial map, providing a
#' barplot for each subpopulation. Therefore it requires coordinates from a
#'  genlight object. This kind of plots should support the interpretation of the
#'   spatial structure of a population, but in principle is not different from
#'   \code{\link{gl.plot.structure}}
#'
#' @param qmat Q-matrix from a structure run followed by a clumpp run object
#' [from \code{\link{gl.run.structure}} and \code{\link{gl.plot.structure}}]
#'  [required].
#' @param x Name of the genlight object containing the coordinates in the
#'  \code{\@other$latlon} slot to calculate the population centers [required].
#' @param K The number for K to be plotted [required].
#' @param provider Provider	passed to leaflet. Check \link[leaflet]{providers}
#' for a list of possible backgrounds [default "Esri.NatGeoWorldMap"].
#' @param scalex Scaling factor to determine the size of the bars in x direction
#' [default 1].
#' @param scaley Scaling factor to determine the size of the bars in y direction
#'  [default 1].
#' @param movepops A two-dimensional data frame that allows to move the center of
#' the barplots manually in case they overlap. Often if populations are
#' horizontally close to each other. This needs to be a data.frame of the
#' dimensions [rows=number of populations, columns = 2 (lon/lat)]. For each
#' population you have to specify the x and y (lon and lat) units you want to
#' move the center of the plot, (see example for details) [default NULL].
#' @param pop.labels Switch for population labels below the parplots
#' [default TRUE].
#' @param pop.labels.cex Size of population labels [default 12].
#' @return An interactive map that shows the structure plots broken down by
#' population.
#' @author Bernd Gruber (Post to \url{https://groups.google.com/d/forum/dartr})
#' @details
#' Creates a mapped version of structure plots. For possible background maps
#' check as specified via the provider:
#' \url{http://leaflet-extras.github.io/leaflet-providers/preview/index.html}.
#' You may need to adjust scalex and scaley values [default 1], as the size
#' depends on the scale of the map and the position of the populations.
#' @return returns the map and a list of the qmat split into sorted matrices per
#'  population. This can be used to create your own map.
#' @examples
#' # examples need structure to be installed on the system (see above)
#' \dontrun{
#' bc <- bandicoot.gl[,1:100]
#' sr <- gl.run.structure(bc, k.range = 2:5, num.k.rep = 3, exec = './structure.exe')
#' ev <- gl.evanno(sr)
#' ev
#' qmat <- gl.plot.structure(sr, k=2:4)#' #head(qmat)
#' gl.map.structure(qmat, bc,K=3)
#' gl.map.structure(qmat, bc,K=4)
#' # move population 4 (out of 5) 0.5 degrees to the right and populations 1
#' # 0.3 degree to the north of the map.
#' mp <- data.frame(lon=c(0,0,0,0.5,0), lat=c(-0.3,0,0,0,0))
#' gl.map.structure(qmat, bc,K=4, movepops=mp)
#' }
#' @export
#' @seealso \code{\link{gl.run.structure}},  \code{clumpp},
#' \code{\link{gl.plot.structure}}
#' @references
#' \itemize{
#' \item Pritchard, J.K., Stephens, M., Donnelly, P. (2000) Inference of
#' population structure using multilocus genotype data. Genetics 155, 945-959.
#' \item Archer, F. I., Adams, P. E. and Schneiders, B. B. (2016) strataG: An R
#'  package for manipulating, summarizing and analysing population genetic data.
#'   Mol Ecol Resour. doi:10.1111/1755-0998.12559
#' \item Evanno, G., Regnaut, S., and J. Goudet. 2005. Detecting the number of
#' clusters of individuals using the software STRUCTURE: a simulation study.
#' Molecular Ecology 14:2611-2620.
#' \item Mattias Jakobsson and Noah A. Rosenberg. 2007. CLUMPP: a cluster
#' matching and permutation program for dealing with label switching and
#' multimodality in analysis of population structure. Bioinformatics
#' 23(14):1801-1806. Available at
#' \href{http://web.stanford.edu/group/rosenberglab/clumppDownload.html}{clumpp}
#' }

gl.map.structure <- function(qmat,
                             x,
                             K,
                             provider = "Esri.NatGeoWorldMap",
                             scalex = 1,
                             scaley = 1,
                             movepops = NULL,
                             pop.labels = TRUE,
                             pop.labels.cex = 12) {
  eq.k <- sapply(qmat, function(x) {
    ncol(x) - 4 == K
  })

  if (sum(eq.k) == 0) {
    stop(error(paste(
      "No entries for K =", K, "found in 'qmat'.\n"
    )))
  }

  qmat <- as.data.frame(qmat[eq.k][[1]])

  # ff <- qmat[, 4:(ncol(qmat))]
  ff <- qmat[, which(grepl("cluster", colnames(qmat)))]

  df <- x@other$latlon
  centers <-
    apply(df, 2, function(xx) {
      tapply(xx, pop(x), mean, na.rm = TRUE)
    })

  if (!is.null(movepops)) {
    if (nrow(movepops) != nrow(centers)) {
      stop(
        error(
          "The provided movepops data.frame has not the corret number of rows, please check. It needs to have the same numbers of rows as the number populations in your genlight object."
        )
      )
    }
    centers[, 1] <- centers[, 1] + movepops[, 1]
    centers[, 2] <- centers[, 2] + movepops[, 2]
  }


  sc <-
    match(rownames(centers), levels(factor(qmat$orig.pop)))
  if (any(is.na(sc))) {
    message(
      error(
        "Population names (coordinates) in the genlight object do not match population in your q-matrix. Please check both."
      )
    )
  }
  centers <- centers[sc, ]
  cx <- centers[, "lon"]
  cy <- centers[, "lat"]
  sx <- abs(diff(range(centers[, "lon"]))) / (100) * scalex
  sy <- 20 * sx * scaley
  #
  qmat$orig.pop <- factor(qmat$orig.pop)
  npops <- length(levels(qmat$orig.pop))
  ll <- data.frame(cbind(as.numeric(qmat$orig.pop), ff))
  zz <- do.call(order, unname(as.list(ll)))
  bb <- qmat[zz, ]
  bb$orig.pop <- factor(bb$orig.pop)
  # ff <- bb[, 4:(ncol(bb))]

  ff <- bb[, which(grepl("cluster", colnames(bb)))]

  out <- list()
  m1 <- leaflet::leaflet() %>%
    leaflet::addProviderTiles(provider = provider)
  for (p in 1:npops) {
    qmi <- ff[bb$orig.pop == levels(bb$orig.pop)[p], ]
    out[[p]] <- bb[bb$orig.pop == levels(bb$orig.pop)[p], ]
    names(out)[p] <- levels(bb$orig.pop)[p]
    qmi1 <- cbind(rep(0, nrow(qmi)), qmi)
    for (xx in 1:nrow(qmi1)) {
      qmi1[xx, ] <- cumsum(as.numeric(qmi1[xx, ]))
    }


    for (ii in 1:nrow(qmi1)) {
      for (i in 1:(ncol(qmi1) - 1)) {
        oo <- (ii - nrow(qmi) / 2) * sx

        m1 <- m1 %>%
          leaflet::addRectangles(
            cx[p] + oo,
            cy[p] + qmi1[ii, i] * sy,
            cx[p] + oo + sx,
            cy[p] + qmi1[ii, i + 1] * sy,
            opacity = 0,
            color = rainbow(ncol(ff))[i],
            fillOpacity = 0.8
          )
      }
    }
  }
  if (pop.labels) {
    m1 <- m1 %>%
      leaflet::addLabelOnlyMarkers(
        lng = centers[, "lon"],
        lat = centers[, "lat"] - sy * 0.1,
        label = rownames(centers),
        labelOptions = leaflet::labelOptions(
          noHide = T,
          direction = "center",
          textOnly = T,
          textsize = paste0(pop.labels.cex, "px")
        )
      )
  }


  print(m1)
  # mapshot(m1, file='./Rplot.png', remove_controls = TRUE)
  return(list(qmats=out,map=m1))
  # %>% addLegend(labels=paste('Group',1:ncol(ff)), colors=rainbow(ncol(ff)),position ='topright' )

  # if (save) mapshot(m1, file='./Rplot.png', remove_controls = TRUE)
}