Add split.by parameter to PlotFootprint; #523

R/visualization.R

@@ -232,6 +232,9 @@ globalVariables(
 #' @param features A vector of features to plot
 #' @param assay Name of assay to use
 #' @param group.by A grouping variable
+#' @param split.by A metadata variable to split the plot by. For example,
+#' grouping by "celltype" and splitting by "batch" will create separate plots
+#' for each celltype and batch.
 #' @param idents Set of identities to include in the plot
 #' @param show.expected Plot the expected Tn5 integration frequency below the
 #' main footprint plot
@@ -256,6 +259,7 @@ PlotFootprint <- function(
   features,
   assay = NULL,
   group.by = NULL,
+  split.by = NULL,
   idents = NULL,
   label = TRUE,
   repel = TRUE,
@@ -265,17 +269,33 @@ PlotFootprint <- function(
   label.idents = NULL
 ) {
   assay <- SetIfNull(x = assay, y = DefaultAssay(object = object))
+  splitby_str <- "__signac_tmp__"
   if (!inherits(x = object[[assay]], what = "ChromatinAssay")) {
     stop("The requested assay is not a ChromatinAssay.")
   }
+  if (!is.null(x = split.by)) {
+    if (is.null(x = group.by)) {
+      grouping.var <- Idents(object = object)
+    } else {
+      grouping.var <- object[[group.by]][, 1]
+    }
+    # combine split.by and group.by information
+    combined.var <- paste0(object[[split.by]][, 1], splitby_str, grouping.var)
+    object$grouping_tmp <- combined.var
+    group.by <- "grouping_tmp"
+  }
   # TODO add option to show variance among cells
   plot.data <- GetFootprintData(
     object = object,
     features = features,
     assay = assay,
     group.by = group.by,
+    split.by = split.by,
     idents = idents
   )
+  if (length(x = plot.data) == 1) {
+    stop("Footprinting data not found")
+  }
   motif.sizes <- GetMotifSize(
     object = object,
     features = features,
@@ -311,6 +331,13 @@ PlotFootprint <- function(
       stop("Unknown normalization method requested")
     }
   }
+
+  # split back into group.by and split.by
+  if (!is.null(x = split.by)) {
+    splitvar <- strsplit(x = obs[['group']], split = splitby_str)
+    obs[['group']] <- sapply(X = splitvar, FUN = `[[`, 2)
+    obs[['split']] <- sapply(X = splitvar, FUN =`[[`, 1)
+  }
 
   # find flanking accessibility for each group and each feature
   flanks <- obs[obs$flanks, ]
@@ -369,6 +396,9 @@ PlotFootprint <- function(
       ggtitle(label = features[[i]]) +
       ylim(c(axis.min, axis.max)) +
       guides(color = guide_legend(override.aes = list(size = 1)))
+    if (!is.null(x = split.by)) {
+      p <- p + facet_wrap(facets = ~split)
+    }
     if (label) {
       if (repel) {
         if (!requireNamespace(package = "ggrepel", quietly = TRUE)) {
@@ -385,26 +415,30 @@ PlotFootprint <- function(
       }
     }
     if (show.expected) {
-      df <- expect[expect$feature == features[[i]], ]
-      p1 <- ggplot(
-        data = df,
-        mapping = aes(x = position, y = norm.value)
-      ) +
-        geom_line(size = 0.2) +
-        xlab("Distance from motif") +
-        ylab(label = "Expected\nTn5 enrichment") +
-        theme_classic()
-
-      # remove x-axis labels from top plot
-      p <- p + theme(
-        axis.title.x = element_blank(),
-        axis.text.x = element_blank(),
-        axis.line.x.bottom = element_blank(),
-        axis.ticks.x.bottom = element_blank()
-      )
-      p <- p + p1 + plot_layout(ncol = 1, heights = c(3, 1))
-      plotlist[[i]] <- p
+      if (!is.null(x = split.by)) {
+        warning("Cannot plot expected enrichment with split.by")
+      } else {
+        df <- expect[expect$feature == features[[i]], ]
+        p1 <- ggplot(
+          data = df,
+          mapping = aes(x = position, y = norm.value)
+        ) +
+          geom_line(size = 0.2) +
+          xlab("Distance from motif") +
+          ylab(label = "Expected\nTn5 enrichment") +
+          theme_classic()
+
+        # remove x-axis labels from top plot
+        p <- p + theme(
+          axis.title.x = element_blank(),
+          axis.text.x = element_blank(),
+          axis.line.x.bottom = element_blank(),
+          axis.ticks.x.bottom = element_blank()
+        )
+        p <- p + p1 + plot_layout(ncol = 1, heights = c(3, 1))
+      }
     }
+    plotlist[[i]] <- p
   }
   plots <- wrap_plots(plotlist)
   return(plots)