ts viz

Author: microsud

.Rhistory

@@ -1,7 +1,7 @@
 Package: microbiomeutilities
 Type: Package
-Title: microbiomeutilities: An R package for utilities to guide in-depth marker gene amplicon data analysis
-Version: 1.00.08
+Title: microbiomeutilities: Utilities for Microbiome Analytics
+Version: 1.00.09
 Authors@R: c(
     person("Sudarshan", "Shetty", email = "[email protected]", role = c("aut", "cre")),
     person("Leo", "Lahti", role = "aut"))

DESCRIPTION

@@ -14,12 +14,15 @@ export(plasticity)
 export(plot_abund_prev)
 export(plot_alpha_diversities)
 export(plot_alpha_rcurve)
+export(plot_area)
 export(plot_diversity_stats)
 export(plot_listed_taxa)
 export(plot_ordination_utils)
 export(plot_ordiplot_core)
+export(plot_paired_abundances)
 export(plot_read_distribution)
 export(plot_select_taxa)
+export(plot_spaghetti)
 export(plot_taxa_boxplot)
 export(plot_taxa_composition)
 export(plot_taxa_cv)
@@ -55,3 +58,4 @@ importFrom(stats,median)
 importFrom(stats,quantile)
 importFrom(stats,sd)
 importFrom(tibble,rownames_to_column)
+importFrom(tidyr,pivot_longer)

NAMESPACE

@@ -25,3 +25,33 @@
 #' }
 #'
 "zackular2014"
+
+#' Test data 2
+#'
+#' Data from a Stansfield J, Dozmorov M. "16s rRNA sequencing data 
+#' from the Human Microbiome Project 2". Randomly choosen 13 particpants
+#' with multiple timepoints for rectum samples.
+#'
+#' @docType data
+#'
+#' @usage data("hmp2")
+#'
+#' @format An object of class \code{"phyloseq"}.
+#'
+#' @keywords datasets
+#' @references
+#' \itemize{
+#' \item{}{Stansfield J, Dozmorov M (2019). HMP2Data: 16s rRNA sequencing data 
+#' from the Human Microbiome Project 2. R package version 1.1.0,
+#' \url{	https://bioconductor.org/packages/HMP2Data/}}
+#' }
+#' @examples
+#' \dontrun{
+#' library(microbiomeutilities)
+#' data("hmp2")
+#' pseq <- hmp2
+#' print(hmp2)
+#' }
+#'
+"hmp2"
+

R/data.R

@@ -0,0 +1,65 @@
+#' @title Area plot
+#' @description Create an area plot for longitudinal samples 
+#' with \code{\link{ggplot2}} package.
+#' @param x \code{\link{phyloseq-class}} object.
+#' @param xvar Column name to plot on x-axis.
+#' @param facet.by Column with variable that has multiple measurements.
+#' @param fill.colors brewer.pal(6,"Paired"). Specify colors.
+#' @param abund.thres = 0.001 check \code{\link{microbiome}} package aggregate_rare function.
+#' @param prev.thres = 0.1 check \code{\link{microbiome}} package aggregate_rare function.
+#' @param level Taxonomic level. OTU/ASV level not supported. 
+#' @param ncol wrap, specify number of columns.
+#' @param nrow wrap, specify number of rows.
+#' @return  \code{\link{ggplot}} object.
+#' @export 
+#' @examples 
+#' \dontrun{
+#' library(microbiomeutilities)
+#' data("hmp2")
+#' ps <- hmp2
+#' ps.rel <- microbiome::transform(ps, "compositional") 
+#' p <- plot_area(ps.rel, xvar="visit_number", 
+#'               level = "Phylum",
+#'               facet.by = "subject_id",
+#'               fill.colors=brewer.pal(6,"Paired"))
+#'}           
+#' 
+#' @keywords visualization
+#' 
+plot_area <- function(x,
+                      xvar=NULL,
+                      level = NULL,
+                      facet.by = NULL,
+                      fill.colors=brewer.pal(6,"Paired"),
+                      abund.thres=0.001,
+                      prev.thres=0.5,
+                      ncol=5,
+                      nrow=5){
+  
+  if(is.null(xvar)){
+    stop("xvar cannot be empty")
+  }
+  if(is.null(facet.by)){
+    stop("facet.by cannot be empty")
+  }
+  x.lev<-xdf<-comp.plt<-facet_var <- NULL
+  x.lev <- aggregate_rare(x, 
+                          level=level,
+                          detection = abund.thres, 
+                          prevalence = prev.thres)
+  
+  xdf <- phy_to_ldf(x.lev, transform.counts = NULL)
+  xdf$facet_var <- xdf[,facet.by]
+  comp.plt <- ggplot(xdf)
+  comp.plt <- comp.plt + 
+    geom_area(aes_string(x = xvar, 
+                         y = "Abundance", 
+                         fill = level)) +
+    facet_wrap(~facet_var, 
+               scales = "free",
+               ncol = ncol,
+               nrow = nrow) +
+   scale_fill_manual("Taxa", values = fill.colors)
+  return(comp.plt + theme_biome_utils() + ylab("Relative abundance"))
+  
+}

R/plot_area.R

@@ -51,12 +51,17 @@ plot_listed_taxa <- function(x,
                              ncol = NULL,
                              nrow = NULL) {
   x.rel <- x.prun <- x.df <- p.box <- p.vio <- p.strp <- NULL
-
+  
+  if(is.na(group) | is.null(group)) {
+    stop(" 'group' argument cannot be empty")
+  }
+  
   x.rel <- microbiome::transform(x, "compositional")
 
   x.prun <- prune_taxa(select.taxa, x.rel)
 
   x.df <- phy_to_ldf(x.prun, transform.counts = NULL)
+  
 
   if (!is.null(group.order)) {
     x.df[, group] <- factor(x.df[, group],

R/plot_listed_taxa.R

@@ -0,0 +1,218 @@
+#' @title A paired-boxplot for user specified list of taxa
+#' @description User specified taxa are plotted.
+#' @details Useful for instances where user is interested only in some taxa and thier change after
+#'          an intervention. This can also be used at higher taxonomic
+#'          levels, output from phyloseq::tax_glom or microbiome::aggregate_taxa.
+#'
+#' @param x \code{\link{phyloseq-class}} object.
+#' @param select.taxa a character list of taxa to be plotted. eg. select.taxa <- c("OTU-370251", "OTU-311173", "OTU-341024").
+#' @param group Grouping variable to compare. x axis, eg. before-after, t1-t2.
+#' @param group.colors Colors for plotting groups.
+#' @param dot.opacity For ggplot alpha to determine opacity for points.
+#' @param dot.size For ggplot point size.
+#' @param jitter.width Value to avoid over plotting by moving points.
+#' @param add.box Logical. If boxplot to the added. Default=TRUE
+#' @param box.opacity For ggplot alpha to determine opacity for box.
+#' @param add.violin Logical. If half violin to the added. Default=TRUE
+#' @param violin.opacity If add.violin=TRUE, opacity for violin.
+#' @param group.order Default is NULL. a list specifying order of x-axis.
+#' @param ncol If 2 or more taxa to plot, specify number of columns.
+#' @param nrow If 2 or more taxa to plot, specify number of rows.
+#' @param line Variable to use for lines. E.g. "subject" before-after
+#' @param line.down Line Color for change when negative. Decreased abundance.
+#' @param line.stable Line Color for no change.
+#' @param line.up Line Color for change when positive. Increased abundance.
+#' @param line.na.value "grey50" for no/missing observations.
+#' @param line.guide "none" to plot guide for line.
+#' @param line.opacity Line opacity.
+#' @param line.size Size of line to plot.
+#'
+#' @return  \code{\link{ggplot}} object. This can be further modified using ggpubr.
+#' @importFrom tidyr pivot_longer
+#' @export
+#' @examples
+#' \dontrun{
+#' library(microbiome)
+#' library(microbiomeutilities)
+#' library(gghalves)
+#' library(tidyr)
+#' data(peerj32) # Source: https://peerj.com/articles/32/
+#' pseq <- peerj32$phyloseq # Ren
+#' pseq.rel <- microbiome::transform(pseq, "compositional")
+#' select.taxa <- c("Akkermansia", "Dialister")
+#' group.colors <- c("brown3", "steelblue", "grey70")
+#' p <- plot_paired_abundances(pseq.rel,
+#'   select.taxa = select.taxa,
+#'   group = "time",
+#'   group.colors = group.colors,
+#'   dot.opacity = 0.25,
+#'   dot.size = 2,
+#'   group.order = NULL,
+#'   line = "subject"
+#' )
+#' p
+#' }
+#' @keywords visualization
+
+plot_paired_abundances <- function(x,
+                                   select.taxa = NULL,
+                                   group = NULL,
+                                   group.colors = NULL,
+                                   dot.opacity = 0.25,
+                                   dot.size = 2,
+                                   add.box = FALSE,
+                                   box.opacity = 0.25,
+                                   group.order = NULL,
+                                   add.violin = TRUE,
+                                   violin.opacity = 0.25,
+                                   ncol = NULL,
+                                   nrow = NULL,
+                                   line = NULL,
+                                   line.down = "#7209b7",
+                                   line.stable = "#8d99ae",
+                                   line.up = "#14213d",
+                                   line.na.value = "grey50",
+                                   line.guide = "legend",
+                                   line.opacity = 0.25,
+                                   line.size = 1,
+                                   jitter.width = 0) {
+  Abundance <- change <- change.order <- change.sign <- linevar <- NULL
+  # check arguments
+  if (is.na(line) | is.null(line)) {
+    stop(" 'line' argument cannot be empty")
+  }
+
+  if (is.na(group) | is.null(group)) {
+    stop(" 'group' argument cannot be empty")
+  }
+  #x <- pseq.rel
+  xmeta <- meta(x)
+  for (i in select.taxa) {
+    df.tx <- as.data.frame(abundances(x)[i, ])
+    colnames(df.tx) <- i
+    xmeta <- cbind(xmeta, df.tx)
+  }
+
+  if(length(unique(xmeta[, group])) > 2){
+    stop("Only two group comparison e.g. before n after")
+  }
+  
+  # check if factor else convert to factor
+  if (!is.factor(xmeta[, group])) {
+    xmeta$group <- factor(as.character(xmeta[, group]))
+  }
+
+  # convert to wide format
+  xmeta_lf <- xmeta %>%
+    pivot_longer(
+      cols = all_of(select.taxa),
+      names_to = "taxa",
+      values_to = "Abundance"
+    )
+
+  xmeta_lf$linevar <- factor(xmeta_lf[[line]])
+
+  x.grp <- sym(group)
+
+  df2 <- suppressWarnings(xmeta_lf %>%
+    arrange(taxa,linevar, !!x.grp) %>%
+    group_by(taxa,linevar) %>%
+    summarise(change = diff(Abundance)))
+
+  xmeta_lf_2 <- suppressWarnings(xmeta_lf %>%
+                     arrange(taxa,linevar, !!x.grp) %>%
+                     group_by(taxa,linevar))
+  
+  xmeta_lf_2 <- xmeta_lf_2 %>% 
+    left_join(df2)
+  
+  #xmeta_lf$change <- df2$change[match(xmeta_lf$linevar, df2$linevar)]
+  
+  xmeta_lf_2$change.sign <- sign(xmeta_lf_2$change)
+  
+  if (!is.null(group.order)) {
+    xmeta_lf_2[, group] <- factor(xmeta_lf_2[, group],
+                                  levels = group.order
+    )
+  }
+  xmeta_lf_2 <- xmeta_lf_2 %>%
+    mutate(change.order = ifelse(change.sign == 1, "Up",
+      ifelse(change.sign == -1, "Down",
+        "Stable"
+      )
+    ))
+  # start plotting
+  p <- ggplot(
+    data = xmeta_lf_2,
+    aes_string(
+      x = "group",
+      y = "Abundance",
+      fill = "group"
+    )
+  ) +
+    geom_point(aes_string(
+      x = "group",
+      fill = "group"
+    ),
+    position = position_jitter(width = jitter.width),
+    size = dot.size,
+    alpha = dot.opacity,
+    shape = 21
+    ) +
+    geom_line(aes(
+      group = linevar,
+      color = change.order
+    ),
+    size = line.size,
+    alpha = line.opacity
+    #lty = line.type
+    )
+
+  p <- p + scale_fill_manual(values = group.colors)
+
+  if (add.box == TRUE) {
+    p <- p + geom_boxplot(
+      width = 0.2,
+      outlier.shape = NA,
+      alpha = box.opacity
+    )
+  }
+
+#geom_half_violin(
+#  data = iris %>% filter(Species=="versicolor"), 
+#  aes(x = Species, y = Sepal.Length, fill = Species), side="r")
+  if (add.violin == TRUE) {
+    p <- p + geom_half_violin(data = xmeta_lf_2 %>% 
+                                filter(group==unique(xmeta_lf_2$group)[1]),
+      position = position_nudge(x = -0.15, y = 0),
+      alpha = violin.opacity,
+      side = "l"
+    ) +
+      geom_half_violin(data = xmeta_lf_2 %>% 
+                         filter(group==unique(xmeta_lf_2$group)[2]),
+                       position = position_nudge(x = 0.15, y = 0),
+                       alpha = violin.opacity,
+                       side = "r"
+      )
+  }
+
+  if (length(select.taxa) >= 2) {
+    p <- p + facet_wrap(~taxa,
+      scales = "free",
+      ncol = ncol,
+      nrow = nrow
+    )
+  }
+  p <- p + scale_color_manual(
+    values = c(
+      Up = line.up,
+      Down = line.down,
+      Stable = line.stable
+    ),
+    guide = line.guide
+  )
+
+  p <- p + theme_biome_utils() + xlab(group)
+
+  return(p)
+}