R/LocalPlot.R

In AdmixPoly: Global and Local Admixture Inference in Polyploids

#' Plot local ancestry dosages
#'
#' @param Ancestry
#' List of local ancestry dosages for an individual of size `C` (number of chromosomes)
#' with each element being a matrix with ancestry dosages of dimension
#' `M` (number of markers) x `K` (number of groups)
#' @param GeneticMap
#' Dataframe with three columns: `Chromosome`, `Marker`, and `Distance`,
#' with `M` rows. Distance must be in centiMorgan
#' @param DisplayNAs
#' Should missing data be displayed as gaps (FALSE by default)
#' 
#' @return
#' A local admixture plot as a ggplot object
#'
#' @description
#' This function generates local admixture plot
#' 
#' @details
#' The function `LocalPlot()` based on [ggplot2::ggplot()] performs local 
#' admixture plot using local ancestry dosage estimated from the [AdmixLocal()] 
#' function and a genetic map.
#' 
#' The plot displays local ancestry dosages for one individual, with one facet 
#' for each of the chromosomes. The x-axis represents the genetic
#' distance and the y-axis represents ancestry dosages, ranging from 0 to 
#' the ploidy level of the individual.
#' 
#' When the individual has missing genotypic data as compared to the genetic map,
#' the gaps can be displayed as blanks by setting `DisplayNAs=TRUE`
#' 
#'
#' @import tibble
#' @import dplyr
#' @import tidyr
#' @import ggplot2
#' @importFrom magrittr %>%
#' @importFrom purrr imap_dfr
#'
#' @seealso
#' * [SimulatePop()] to simulate a polyploid admixed population.
#' * [AdmixGlobal()] to perform global (genome-wide) admixture inference.
#' * [AdmixLocal()] to perform local admixture inference.
#'
#' @export
#'
#' @examples
#' ## Simulate Simulate a polyploid admixed population
#' DataSim <- SimulatePop(K=3L, N=10L, P=6L, M=50L, C=5L, L=10L, Seed=123, NbThreads=1)
#'
#' ## Perform global admixture inference
#' ResAdmixGlobal <- AdmixGlobal(Geno=DataSim$Geno, K=3, Verbose=FALSE, NbThreads=1)
#'
#' ## Perform local admixture inference for one individual
#' ResAdmixLocal <- AdmixLocal(Geno=DataSim$Geno, ResAdmixGlobal, "Ind4", 6L,
#'                             DataSim$GeneticMap, Verbose=FALSE, NbThreads=1)
#'                             
#' ## Local admixture barplot 
#' LocalPlot(ResAdmixLocal$Posterior, DataSim$GeneticMap)
LocalPlot <- function(Ancestry, GeneticMap, DisplayNAs = FALSE){
  ## Checks
  stopifnot("Ancestry must be a list" =
              is.list(Ancestry))
  stopifnot("Each element of Ancestry must be numeric matrix with row and column names" =
              all(sapply(Ancestry,function(i)
                is.matrix(i)&&is.numeric(i)&&!is.null(rownames(i))&&!is.null(colnames(i)))))
  stopifnot("GeneticMap must be a dataframe with three named columns
             (Marker, Chromosome and Distance)" =
              is.data.frame(GeneticMap) &&
              all(c("Marker","Chromosome","Distance")%in%colnames(GeneticMap)))
  stopifnot("Ancestry and GeneticMap must have the same markers" =
              all(unlist(sapply(Ancestry,rownames))%in%GeneticMap$Marker))
  stopifnot("Ancestry and GeneticMap must have the same chromosomes" =
              all(names(Ancestry)%in%GeneticMap$Chromosome))
  stopifnot("DisplayNAs must be a boolean" =
              is.logical(DisplayNAs))

  ## Define global variable for tidy operations
  Marker <- NULL
  Chromosome <- NULL
  Distance <- NULL
  Group <- NULL
  Dosage <- NULL
  ymin <- NULL
  ymax <- NULL
  . <- NULL

  ## Format dataframe
  Ancestry_df <- imap_dfr(Ancestry,~as_tibble(.x) %>%
                            mutate(Marker=rownames(.x), Chromosome=.y)) %>%
    {
      if (DisplayNAs) {
        full_join(x = ., y = GeneticMap, by = c("Chromosome", "Marker"))
      } else {
        left_join(x = ., y = GeneticMap, by = c("Chromosome", "Marker"))
      }
    } %>%
    # left_join(GeneticMap,by=c("Chromosome","Marker")) %>%
    pivot_longer(!Marker&!Chromosome&!Distance,names_to = "Group",values_to = "Dosage") %>%
    group_by(Marker) %>%
    reframe(Group=Group,Dosage=Dosage,Distance=Distance,Chromosome=Chromosome,
            ymin=c(0,cumsum(Dosage)[-length(Dosage)]),
            ymax=cumsum(Dosage),
            .groups = "keep")

  ## Local admixture plot
  p <- ggplot(Ancestry_df,aes(x=Distance, ymin=ymin, ymax=ymax, fill=Group)) +
    facet_wrap(~Chromosome,ncol = 1) +
    scale_y_continuous(expand=expansion(mult = c(0,0))) +
    scale_x_continuous(expand=expansion(mult = c(0,0))) +
    scale_fill_brewer(palette = "Set1") +
    xlab("Distance") +
    ylab("Ancestry dosage") +
    theme_bw() +
    geom_ribbon(alpha=0.7) +
    theme(panel.grid = element_blank())

  ## Outputs
  return(p)
}