Merge pull request #48 from FertigLab/krishnan_pkgdown

Update pkgdown files, reorganize data files

.gitignore

@@ -10,4 +10,3 @@ R/.Rhistory
 scenic/*
 
 scratch/
-inst/

DESCRIPTION

@@ -9,7 +9,7 @@ Authors@R: c(
     person("Elana", "Fertig", role = "ctb", email = "[email protected]", comment = c(ORCID = "0000-0003-3204-342X")),
     person("Jennifer", "Elisseeff", role = "ctb", email = "[email protected]", comment = c(ORCID = "0000-0002-5066-1996"))
     )
-Description: Domino2 is a package developed to analyze cell signaling through ligand - receptor - transcription factor networks in scRNAseq data.
+Description: Domino2 is a package developed to analyze cell signaling through ligand - receptor - transcription factor networks in scRNAseq data. It takes as input information transcriptomic data, requiring counts, z-scored counts, and cluster labels, as well as information on transcription factor activation (such as from SCENIC) and a database of ligand and receptor pairings (such as from cellphoneDB). This package creates an object storing ligand - receptor - transcription factor linkages by cluster and provides several methods for exploring, summarizing, and visualizing the analysis.
 BugReports: https://github.com/FertigLab/domino_development/issues
 Depends:
     R(>= 3.6.2),

NAMESPACE

@@ -33,7 +33,6 @@ import(plyr)
 import(stats)
 importClassesFrom(Matrix,dgCMatrix)
 importFrom(Matrix,rowSums)
-importFrom(NMF,aheatmap)
 importFrom(ggpubr,ggscatter)
 importFrom(igraph,E)
 importFrom(igraph,V)

NEWS.md

@@ -1,4 +1,14 @@
-# Domino2 v0.2.1
+# domino2 v0.2.2
+
+## Linkage functions
+- Addition of new class to summarize linkages in objects
+- Addition of helper functions to count linkages and compare between objects
+- Plotting function for differential linkages
+
+## Package structure
+- Adjustments made to meet BioConductor standards
+
+# domino2 v0.2.1
 
 ## Updates to domino object construction
 - Uniform formats for inputs of receptor-ligand interaction databases, tf activity features, and regulon gene lists for operability with alternative databases and tf activity inference methods

R/class_definitions.R

@@ -2,7 +2,7 @@
 #' @importClassesFrom Matrix dgCMatrix
 #'
 NULL
-#' The Domino Class
+#' The domino Class
 #' 
 #' The domino class contains all information necessary to calculate receptor-ligand
 #' signaling. It contains z-scored expression, cell cluster labels, feature values,
@@ -43,7 +43,7 @@ domino <- methods::setClass(
     misc = list("build"=FALSE)
   )
 )
-#' The Domino linkage summary class
+#' The domino linkage summary class
 #' 
 #' The linkage summary class contains linkages established in multiple domino
 #' objects through gene regulatory network inference and reference to receptor-
@@ -59,8 +59,15 @@ domino <- methods::setClass(
 #' @rdname linkage_summary-class
 #' @exportClass linkage_summary
 #' 
-linkage_summary <- setClass(Class="linkage_summary", slots=c(subject_names="factor", subject_meta="data.frame",
-  subject_linkages="list"))
+linkage_summary <- setClass(
+  Class = "linkage_summary",
+  slots = c(
+    subject_names = "factor",
+    subject_meta = "data.frame",
+    subject_linkages = "list"
+  )
+)
+
 #' Print domino object
 #' 
 #' Prints a summary of a domino object

R/convenience_fxns.R

@@ -16,6 +16,7 @@ NULL
 #' @param dom Domino object to rename clusters in
 #' @param clust_conv Named vector of conversions from old to new clusters. Values are taken as new clusters IDs and names as old cluster IDs.
 #' @return A domino object with clusters renamed in all applicable slots.
+#' @keywords internal
 #' @export 
 #' 
 rename_clusters <- function(dom, clust_conv) {
@@ -100,7 +101,9 @@ collate_network_items <- function(dom, clusters = NULL, return = NULL) {
 #' @param conversion_table A data.frame with column names corresponding to gene symbol types (mm.ens, hs.ens, mgi, hgnc)
 #' and rows corresponding to the gene symbols themselves
 #' @return Data frame of genes with original and corresponding converted symbols
+#' @keywords internal
 #' @export
+#' 
 table_convert_genes <- function(genes, from, to, conversion_table) {
   # Check inputs:
   stopifnot(`Genes must be a vector of characters` = (is(genes, "character") & is(genes, "vector")))

R/data.R

@@ -1,26 +1,12 @@
-#' Peripheral Blood Mononuclear Cell 3K (PBMC3K) data set
-#' 
-#' DESC ...
-#' 
-#' @format 'pbmc'
-#' A Seurat object with 13,714 features in 2700 cells
-#' \describe{
-#'   \item{orig.ident}
-#'   \item{nCount_RNA}
-#'   \item{nFeature_RNA}
-#'   \item{percent.mt}
-#'   \item{RNA_snn_res.0.5, seurat_clusters}
-#'   \item{cell_types}
-#' }
-
-#' pySCENIC AUC Result
-#' 
-#' DESC
-#' 
-#' @format 'auc_test'
-#' Data frame of transcription factor activity scores for 3 transcription factors in 300 cells.
-#' \describe{
-#'    \item{columns}{cells}
-#'    \item{rows}{transcription factors}
-#' }
-
+#' Example domino object from pbmc3k data
+#'
+#' A small example single cell data set from 10x Genomics that has been through preliminary analysis, 
+#' including normalization, scaling, PCA, UMAP, clustering, and cluster annotation and then used in an example
+#' domino2 analysis.
+#'
+#' @format ## `pbmc_dom`
+#' A domino object based on the 10x Genomics pbmc3k data set that serves as an example in tutorials
+#' and can also be used to explore the object structure and functionality.
+#' 
+#' @usage data(pbmc_dom)
+"pbmc_dom"

R/import_fxns.R

@@ -41,8 +41,10 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
     1)
   stopifnot(`Alternate conversion argument (not recommended) must be TRUE or FALSE` = is(alternate_convert,
     "logical"))
-  stopifnot(`If using alternate conversion table (not recommended), table must be provided as data.frame` = (alternate_convert &
-    is(alternate_convert_table, "data.frame")))
+  if(alternate_convert & is.null(alternate_convert_table)) {
+      stop("If using alternate conversion table (not recommended), a table must be provided")
+  }
+
   # Read in files if needed:
   if (is(genes, "character")) {
     genes <- read.csv(genes, stringsAsFactors = FALSE)
@@ -302,9 +304,11 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
     "data.frame") & c("gene_A", "gene_B", "type_A", "type_B") %in% colnames(rl_map)))
   stopifnot(`features must be either a file path or a named matrix with cells as columns and features as rows` = ((is(features,
     "character") & length(features) == 1) | (is(features, "matrix") & !is.null(rownames(features)) &
+    !is.null(colnames(features))) | (is(features, "data.frame") & !is.null(rownames(features)) &
     !is.null(colnames(features)))))
-  stopifnot(`Either a Seurat object OR z scores and clusters must be provided` = (is(ser, "Seurat") |
-    (is(features, "matrix") & !is.null(rownames(features)) & !is.null(colnames(features)) & is(clusters,
+  stopifnot(`Either a Seurat object OR counts, z scores, and clusters must be provided` = (is(ser, "Seurat") |
+    (!is.null(counts) & !is.null(rownames(counts)) & !is.null(colnames(counts)) &
+    is(z_scores, "matrix") & !is.null(rownames(z_scores)) & !is.null(colnames(z_scores)) & is(clusters,
       "factor") & !is.null(names(clusters)))))
   stopifnot(`rec_min_thresh must be a number between 0 and 1` = (is(rec_min_thresh, "numeric") &
     rec_min_thresh <= 1 & rec_min_thresh >= 0))
@@ -563,6 +567,7 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
 #' @param to Format of gene output (MGI, or HGNC)
 #' @param host Host to connect to. Defaults to https://www.ensembl.org following the useMart default, but can be changed to archived hosts if useMart fails to connect.
 #' @return A data frame with input genes as col 1 and output as col 2
+#' @keywords internal
 #' @export
 #' 
 convert_genes <- function(genes, from = c("ENSMUSG", "ENSG", "MGI", "HGNC"), to = c("MGI", "HGNC"),
@@ -649,17 +654,20 @@ add_rl_column <- function(map, map_ref, conv, new_name) {
 #' @param destination Name of the receptor with which each ligand interacts
 #' @return A data frame of ligand expression targeting the specified receptor
 #' @export
-#'
-mean_ligand_expression <- function(x, ligands, cell_ident, cell_barcodes, destination) {
+#' 
+mean_ligand_expression <- function(x, ligands, cell_ident, cell_barcodes, destination){
   # initiate data frame to store results
   df <- NULL
-  for (feat in ligands) {
+  for(feat in ligands){
     # index of ligand row
     lig_index <- grep(paste0("^", feat, "$"), rownames(x))
-    # column indices of cells belonging to cell_ident
+    # column indecies of cells belonging to cell_ident
     cell_index <- colnames(x) %in% cell_barcodes
-    cell_df <- data.frame(origin = paste0(cell_ident, "_", feat), destination = destination, mean.expression = mean(x[lig_index,
-      cell_index]))
+    cell_df <- data.frame(
+      origin = paste0(cell_ident, "_", feat),
+      destination = destination,
+      mean.expression = mean(x[lig_index, cell_index])
+    )
     df <- rbind(df, cell_df)
   }
   return(df)

R/plot_fxns.R

@@ -4,7 +4,6 @@
 #' @importFrom igraph graph V E layout_in_circle layout_on_sphere layout_randomly layout_with_fr layout_with_kk simplify
 #' @importFrom ggpubr ggscatter
 #' @import grDevices
-#' @importFrom NMF aheatmap
 #' 
 NULL
 
@@ -19,8 +18,8 @@ NULL
 #' @param max_thresh Maximum signaling threshold for plotting. Defaults to Inf for no threshold.
 #' @param scale How to scale the values (after thresholding). Options are 'none', 'sqrt' for square root, or 'log' for log10.
 #' @param normalize Options to normalize the matrix. Normalization is done after thresholding and scaling. Accepted inputs are 'none' for no normalization, 'rec_norm' to normalize to the maximum value with each receptor cluster, or 'lig_norm' to normalize to the maximum value within each ligand cluster 
-#' @param ... Other parameters to pass to  [NMF::aheatmap()]
-#' @return an aheatmap rendered to the active graphics device
+#' @param ... Other parameters to pass to  [ComplexHeatmap::Heatmap()]
+#' @return a Heatmap rendered to the active graphics device
 #' @export signaling_heatmap
 #' 
 signaling_heatmap <- function(dom, clusts = NULL, min_thresh = -Inf, max_thresh = Inf, scale = "none",
@@ -51,7 +50,11 @@ signaling_heatmap <- function(dom, clusts = NULL, min_thresh = -Inf, max_thresh
   } else if (normalize != "none") {
     stop("Do not recognize normalize input")
   }
-  aheatmap(mat, ...)
+  Heatmap(
+    mat,
+    name = "collective\nsignaling",
+    ...
+  )
 }
 #' Create a cluster incoming signaling heatmap
 #' 
@@ -69,9 +72,9 @@ signaling_heatmap <- function(dom, clusts = NULL, min_thresh = -Inf, max_thresh
 #' @param max_thresh Maximum signaling threshold for plotting. Defaults to Inf for no threshold.
 #' @param scale How to scale the values (after thresholding). Options are 'none', 'sqrt' for square root, or 'log' for log10.
 #' @param normalize Options to normalize the matrix. Accepted inputs are 'none' for no normalization, 'rec_norm' to normalize to the maximum value with each receptor cluster, or 'lig_norm' to normalize to the maximum value within each ligand cluster 
-#' @param title Either a string to use as the title or a boolean describing whether to include a title. In order to pass the 'main' parameter to  [NMF::aheatmap()]  you must set title to FALSE.
-#' @param ... Other parameters to pass to  [NMF::aheatmap()]. Note that to use the 'main' parameter of  [NMF::aheatmap()]  you must set title = FALSE
-#' @return an aheatmap rendered to the active graphics device
+#' @param title Either a string to use as the title or a boolean describing whether to include a title. In order to pass the 'main' parameter to  [ComplexHeatmap::Heatmap()]  you must set title to FALSE.
+#' @param ... Other parameters to pass to  [ComplexHeatmap::Heatmap()]. Note that to use the 'column_title' parameter of  [ComplexHeatmap::Heatmap()]  you must set title = FALSE
+#' @return a Heatmap rendered to the active graphics device
 #' @export incoming_signaling_heatmap
 #' 
 incoming_signaling_heatmap <- function(dom, rec_clust, clusts = NULL, min_thresh = -Inf, max_thresh = Inf,
@@ -111,11 +114,31 @@ incoming_signaling_heatmap <- function(dom, rec_clust, clusts = NULL, min_thresh
     stop("Do not recognize normalize input")
   }
   if (title == TRUE) {
-    aheatmap(mat, main = paste0("Expression of ligands targeting cluster ", rec_clust), ...)
+    return(
+      Heatmap(
+        mat, 
+        name = "expression",
+        column_title = paste0("Expression of ligands targeting cluster ", rec_clust), 
+        ...
+      )
+    )
   } else if (title == FALSE) {
-    aheatmap(mat, ...)
+    return(
+      Heatmap(
+        mat, 
+        name = "expression",
+        ...
+      )
+    )
   } else {
-    aheatmap(mat, main = title, ...)
+    return(
+      Heatmap(
+        mat,
+        name = "expression",
+        column_title = title, 
+        ...
+      )
+    )
   }
 }
 #' Create a cluster to cluster signaling network diagram
@@ -418,15 +441,15 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, class
 #' @param dom Domino object with network built ([build_domino()])
 #' @param bool Boolean indicating whether the heatmap should be continuous or boolean. If boolean then bool_thresh will be used to determine how to define activity as positive or negative.
 #' @param bool_thresh Numeric indicating the threshold separating 'on' or 'off' for feature activity if making a boolean heatmap.
-#' @param title Either a string to use as the title or a boolean describing whether to include a title. In order to pass the 'main' parameter to  [NMF::aheatmap()]  you must set title to FALSE.
+#' @param title Either a string to use as the title or a boolean describing whether to include a title. In order to pass the 'main' parameter to  [ComplexHeatmap::Heatmap()]  you must set title to FALSE.
 #' @param norm Boolean indicating whether or not to normalize the transcrption factors to their max value.
 #' @param feats Either a vector of features to include in the heatmap or 'all' for all features. If left NULL then the features selected for the signaling network will be shown.
 #' @param ann_cols Boolean indicating whether to include cell cluster as a column annotation. Colors can be defined with cols. If FALSE then custom annotations can be passed to NMF.
 #' @param cols Named vector of colors to annotate cells by cluster color. Values are taken as colors and names as cluster. If left as NULL then default ggplot colors will be generated.
 #' @param min_thresh Minimum threshold for color scaling if not a boolean heatmap
 #' @param max_thresh Maximum threshold for color scaling if not a boolean heatmap
-#' @param ... Other parameters to pass to  [NMF::aheatmap()] . Note that to use the 'main' parameter of  [NMF::aheatmap()]  you must set title = FALSE and to use 'annCol' or 'annColors' ann_cols must be FALSE.
-#' @return an aheatmap rendered to the active graphics device
+#' @param ... Other parameters to pass to  [ComplexHeatmap::Heatmap()] . Note that to use the 'main' parameter of  [ComplexHeatmap::Heatmap()]  you must set title = FALSE and to use 'annCol' or 'annColors' ann_cols must be FALSE.
+#' @return a Heatmap rendered to the active graphics device
 #' @export feat_heatmap
 #' 
 feat_heatmap <- function(dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, title = TRUE, norm = FALSE,
@@ -487,16 +510,45 @@ feat_heatmap <- function(dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, tit
       cols <- ggplot_col_gen(length(levels(cl)))
       names(cols) <- levels(cl)
     }
-    cols <- list(Cluster = cols)
+    # cols <- list(Cluster = cols)
+    feat_anno <- columnAnnotation(
+      Cluster = cl,
+      col = list(Cluster = cols)
+    )
   }
   if (title != FALSE & ann_cols != FALSE) {
-    aheatmap(mat, Colv = NA, annCol = ac, annColors = cols, main = title, ...)
+    Heatmap(
+      mat,
+      name = "feature\nactivity",
+      top_annotation = feat_anno,
+      cluster_columns = FALSE, show_column_names = FALSE,
+      column_title = title,
+      ...
+    )
   } else if (title == FALSE & ann_cols != FALSE) {
-    aheatmap(mat, Colv = NA, annCol = ac, annColors = cols, ...)
+    Heatmap(
+      mat,
+      name = "feature\nactivity",
+      top_annotation = feat_anno,
+      cluster_columns = FALSE, show_column_names = FALSE,
+      ...
+    )
   } else if (title != FALSE & ann_cols == FALSE) {
-    aheatmap(mat, Colv = NA, main = title, ...)
+    Heatmap(
+      mat,
+      name = "feature\nactivity",
+      top_annotation = feat_anno,
+      cluster_columns = FALSE, show_column_names = FALSE,
+      column_title = title,
+      ...
+    )
   } else if (title == FALSE & ann_cols == FALSE) {
-    aheatmap(mat, Colv = NA, ...)
+    Heatmap(
+      mat,
+      name = "feature\nactivity",
+      cluster_columns = FALSE, show_column_names = FALSE,
+      ...
+    )
   }
 }
 #' Create a heatmap of correlation between receptors and transcription factors
@@ -507,12 +559,12 @@ feat_heatmap <- function(dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, tit
 #' @param dom Domino object with network built ([build_domino()])
 #' @param bool Boolean indicating whether the heatmap should be continuous or boolean. If boolean then bool_thresh will be used to determine how to define activity as positive or negative.
 #' @param bool_thresh Numeric indicating the threshold separating 'on' or 'off' for feature activity if making a boolean heatmap.
-#' @param title Either a string to use as the title or a boolean describing whether to include a title. In order to pass the 'main' parameter to  [NMF::aheatmap()]  you must set title to FALSE.
+#' @param title Either a string to use as the title or a boolean describing whether to include a title. In order to pass the 'main' parameter to  [ComplexHeatmap::Heatmap()]  you must set title to FALSE.
 #' @param feats Either a vector of features to include in the heatmap or 'all' for all features. If left NULL then the features selected for the signaling network will be shown.
 #' @param recs Either a vector of receptors to include in the heatmap or 'all' for all receptors. If left NULL then the receptors selected in the signaling network connected to the features plotted will be shown.
 #' @param mark_connections Boolean indicating whether to add an 'x' in cells where there is a connected receptor or TF. Default FALSE.
-#' @param ... Other parameters to pass to  [NMF::aheatmap()] . Note that to use the 'main' parameter of  [NMF::aheatmap()]  you must set title = FALSE and to use 'annCol' or 'annColors' ann_cols must be FALSE.
-#' @return an aheatmap rendered to the active graphics device
+#' @param ... Other parameters to pass to  [ComplexHeatmap::Heatmap()] . Note that to use the 'main' parameter of  [ComplexHeatmap::Heatmap()]  you must set title = FALSE and to use 'annCol' or 'annColors' ann_cols must be FALSE.
+#' @return a Heatmap rendered to the active graphics device
 #' @export cor_heatmap
 #' 
 cor_heatmap <- function(dom, bool = FALSE, bool_thresh = 0.15, title = TRUE, feats = NULL, recs = NULL,
@@ -570,9 +622,24 @@ cor_heatmap <- function(dom, bool = FALSE, bool_thresh = 0.15, title = TRUE, fea
     }
   }
   if (title != FALSE & mark_connections) {
-    aheatmap(mat, main = title, txt = cons, ...)
+    Heatmap(
+      mat,
+      name = "rho",
+      column_title = title,
+      cell_fun = function(j, i, x, y, w, h, col){
+        grid.text(
+          cons[i,j], x, y,
+          gp = gpar(col = "#000000")
+        )
+      },
+      ...
+    )
   } else {
-    aheatmap(mat, ...)
+    Heatmap(
+      mat, 
+      name = "rho",
+      ...
+    )
   }
 }
 #' Create a correlation plot between transcription factor activation score and receptor

R/processing_fxns.R

@@ -126,7 +126,7 @@ build_domino <- function(dom, max_tf_per_clust = 5, min_tf_pval = 0.01, max_rec_
         inc_ligs <- unlist(inc_ligs_list)
       }
       lig_genes <- intersect(inc_ligs, rownames(dom@z_scores))
-      if (length(lig_genes) == 1) {
+      if (length(lig_genes) %in% c(0, 1)) {
         lig_genes <- numeric(0)
       }
       cl_sig_mat <- matrix(0, ncol = length(levels(dom@clusters)), nrow = length(lig_genes))