Merge pull request #62 from FertigLab/krishnan_bioccheck_notes

Add running examples

DESCRIPTION

@@ -13,7 +13,7 @@ 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),
+    R(>= 4.2.0),
 Imports:
     biomaRt,
     ComplexHeatmap,

NAMESPACE

@@ -1,9 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
+export(add_rl_column)
 export(build_domino)
 export(circos_ligand_receptor)
-export(collate_network_items)
-export(convert_genes)
 export(cor_heatmap)
 export(cor_scatter)
 export(count_linkage)
@@ -17,6 +16,7 @@ export(dom_database)
 export(dom_de)
 export(dom_info)
 export(dom_linkages)
+export(dom_network_items)
 export(dom_signaling)
 export(dom_tf_activation)
 export(dom_zscores)
@@ -29,7 +29,6 @@ export(rename_clusters)
 export(signaling_heatmap)
 export(signaling_network)
 export(summarize_linkages)
-export(table_convert_genes)
 export(test_differential_linkages)
 exportClasses(domino)
 exportClasses(linkage_summary)

R/class_definitions.R

@@ -3,12 +3,12 @@
 #'
 NULL
 #' 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,
 #' and a referenced receptor-ligand database formatted as a receptor-ligand map.
 #' Calculated intermediate values are also stored.
-#' 
+#'
 #' @slot db_info List of data sets from lr database.
 #' @slot counts Raw count gene expression data
 #' @slot z_scores Matrix of z-scored expression data with cells as columns
@@ -23,42 +23,42 @@ NULL
 #' @name domino-class
 #' @rdname domino-class
 #' @exportClass domino
-#' 
+#'
 domino <- methods::setClass(
   Class = "domino",
   slots = c(
-    db_info="list",
-    z_scores="matrix",
-    counts="dgCMatrix",
-    clusters="factor",
-    features="matrix",
-    cor="matrix",
-    linkages="list",
-    clust_de="matrix",
-    misc="list",
-    cl_signaling_matrices="list",
-    signaling="matrix"
+    db_info = "list",
+    z_scores = "matrix",
+    counts = "dgCMatrix",
+    clusters = "factor",
+    features = "matrix",
+    cor = "matrix",
+    linkages = "list",
+    clust_de = "matrix",
+    misc = "list",
+    cl_signaling_matrices = "list",
+    signaling = "matrix"
   ),
   prototype = list(
-    misc = list("build"=FALSE)
+    misc = list("build" = FALSE)
   )
 )
 #' 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-
 #' ligand data bases. A data frame summarizing meta features that describe the
 #' domino objects compared in the linkage summary facilitates comparisons of
 #' established linkages and differential signaling interactions across categorical
 #' sample covariates.
-#' 
+#'
 #' @slot subject_names unique names for each domino result included in the summary
 #' @slot subject_meta data.frame with each row describing one subject and columns describing features of the subjects by which to draw comparisons of signaling networks
 #' @slot subject_linkages nested list of linkages inferred for each subject. Lists are stored in a heirarchical structure of subject-cluster-linkage where linkages include transcription factors (tfs), linkages between transcription factors and receptors (tfs_rec), active receptors (rec), possible receptor-ligand interactions (rec_lig), and incoming ligands (incoming_lig)
 #' @name linkage_summary-class
 #' @rdname linkage_summary-class
 #' @exportClass linkage_summary
-#' 
+#'
 linkage_summary <- setClass(
   Class = "linkage_summary",
   slots = c(
@@ -69,39 +69,53 @@ linkage_summary <- setClass(
 )
 
 #' Print domino object
-#' 
+#'
 #' Prints a summary of a domino object
-#' 
+#'
 #' @param x Domino object
 #' @return a printed description of the number of cell clusters in the object
 #' @keywords internal
+#' @examples
+#' print(domino2:::pbmc_dom_built_tiny)
+#' 
 setMethod("print", "domino", function(x, ...) {
   if (x@misc$build) {
-    cat("A domino object of ", length(x@clusters), " cells
+    message(
+      "A domino object of ", length(x@clusters), " cells
                 Contains signaling between",
       length(levels(x@clusters)), "clusters
                 Built with a maximum of", as.integer(x@misc$build_vars["max_tf_per_clust"]),
       "TFs per cluster
                 and a maximum of", as.integer(x@misc$build_vars["max_rec_per_tf"]),
-      "receptors per TF\n")
+      "receptors per TF\n"
+    )
   } else {
-    cat(c("A domino object of ", length(x@clusters), " cells\n", "A signaling network has not been built\n"),
-      sep="")
+    message(c("A domino object of ", length(x@clusters), " cells\n", "A signaling network has not been built\n"),
+      sep = ""
+    )
   }
 })
 #' Show domino object information
-#' 
+#'
 #' Shows content overview of domino object
-#' 
+#'
 #' @param object Domino object
 #' @return a printed description of the number of cells in a domino object and its build status
 #' @keywords internal
+#' @examples
+#' domino2:::pbmc_dom_built_tiny
+#' 
+#' show(domino2:::pbmc_dom_built_tiny)
+#' 
 setMethod("show", "domino", function(object) {
   if (object@misc$build) {
-    cat(c("A domino object of ", length(object@clusters), " cells\n", "Built with signaling between ",
-      length(levels(object@clusters)), " clusters\n"), sep = "")
+    cat(c(
+      "A domino object of ", length(object@clusters), " cells\n", "Built with signaling between ",
+      length(levels(object@clusters)), " clusters\n"
+    ), sep = "")
   } else {
     cat(c("A domino object of ", length(object@clusters), " cells\n", "A signaling network has not been built\n"),
-      sep = "")
+      sep = ""
+    )
   }
 })

R/convenience_fxns.R

@@ -1,24 +1,28 @@
 #' @import plyr
 #' @import methods
-#' 
+#'
 NULL
 
 #' Renames clusters in a domino object
-#' 
-#' This function reads in a receptor ligand signaling database, cell level 
-#' features of some kind (ie. output from pySCENIC), z-scored single cell data, 
-#' and cluster id for single cell data, calculates a correlation matrix between 
-#' receptors and other features (this is transcription factor module scores if 
-#' using pySCENIC), and finds features enriched by cluster. It will return a 
-#' domino object prepared for [build_domino()], which will calculate a signaling 
+#'
+#' This function reads in a receptor ligand signaling database, cell level
+#' features of some kind (ie. output from pySCENIC), z-scored single cell data,
+#' and cluster id for single cell data, calculates a correlation matrix between
+#' receptors and other features (this is transcription factor module scores if
+#' using pySCENIC), and finds features enriched by cluster. It will return a
+#' domino object prepared for [build_domino()], which will calculate a signaling
 #' network.
-#' 
+#'
 #' @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 
-#' 
+#' @export
+#' @examples 
+#' new_clust <- c("CD8_T_cell" = "CD8+ T Cells",
+#'  "CD14_monocyte" = "CD14+ Monocytes", "B_cell" = "B Cells")
+#' pbmc_dom_built_tiny <- rename_clusters(domino2:::pbmc_dom_built_tiny, new_clust)
+#'
 rename_clusters <- function(dom, clust_conv) {
   if (is.null(dom@clusters)) {
     stop("There are no clusters in this domino object")
@@ -43,76 +47,35 @@ rename_clusters <- function(dom, clust_conv) {
   }
   return(dom)
 }
-#' Extracts all features, receptors, or ligands present in a signaling network.
-#' 
-#' This function collates all of the features, receptors, or ligands found in a
-#' signaling network anywhere in a list of clusters. This can be useful for
-#' comparing signaling networks across two separate conditions. In order to run
-#' this [build_domino()] must be run on the object previously.
-#' 
-#' @param dom Domino object containing a signaling network (i.e. [build_domino()] run)
-#' @param return String indicating where to collate 'features', 'receptors', or 'ligands'. If 'all' then a list of all three will be returned.
-#' @param clusters Vector indicating clusters to collate network items from. If left as NULL then all clusters will be included.
-#' @return A vector containing all features, receptors, or ligands in the data set or a list containing all three.
-#' @export 
-#' 
-collate_network_items <- function(dom, clusters = NULL, return = NULL) {
-  if (!dom@misc[["build"]]) {
-    stop("Please run domino_build prior to generate signaling network.")
-  }
-  if (is.null(clusters) & is.null(dom@clusters)) {
-    stop("There are no clusters in this domino object. Please provide clusters.")
-  }
-  if (is.null(clusters)) {
-    clusters <- levels(dom@clusters)
-  }
-  # Get all enriched TFs and correlated + expressed receptors for specified clusters
-  all_recs <- c()
-  all_tfs <- c()
-  all_ligs <- c()
-  for (cl in clusters) {
-    all_recs <- c(all_recs, unlist(dom@linkages$clust_tf_rec[[cl]]))
-    tfs <- names(dom@linkages$clust_tf_rec[[cl]])
-    tf_wo_rec <- which(sapply(dom@linkages$clust_tf_rec[[cl]], length) == 0)
-    if (length(tf_wo_rec > 0)) {
-      tfs <- tfs[-tf_wo_rec]
-    }
-    all_tfs <- c(all_tfs, tfs)
-    all_ligs <- c(all_ligs, rownames(dom@cl_signaling_matrices[[cl]]))
-  }
-  all_recs <- unique(all_recs)
-  all_tfs <- unique(all_tfs)
-  all_ligs <- unique(all_ligs)
-  # Make list and return whats asked for
-  list_out <- list(features = all_tfs, receptors = all_recs, ligands = all_ligs)
-  if (is.null(return)) {
-    return(list_out)
-  } else {
-    return(list_out[[return]])
-  }
-}
+
 #' Convert Genes Using Table
-#' 
-#' Takes a vector of gene inputs and returns converted gene table
-#' 
+#'
+#' Takes a vector of gene inputs and a conversion table ([an example]("http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt")
+#' and returns converted gene table
+#'
 #' @param genes The genes to convert.
 #' @param from  Gene symbol type of the input (ENSG, ENSMUSG, HGNC, MGI)
 #' @param to    Desired gene symbol type for the output (HGNC, MGI)
 #' @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")))
-  stopifnot(`From must be one of ENSMUSG, ENSG, MGI, or HGNC` = from %in% c("ENSMUSG", "ENSG", "MGI",
-    "HGNC"))
+  stopifnot(`From must be one of ENSMUSG, ENSG, MGI, or HGNC` = from %in% c(
+    "ENSMUSG", "ENSG", "MGI",
+    "HGNC"
+  ))
   stopifnot(`To must be one of MGI or HGNC` = to %in% c("MGI", "HGNC"))
-  stopifnot(`Conversion table must be provided with at least two of column names mm.ens, hs.ens, mgi and/or hgnc` = (is(conversion_table,
-    "data.frame") & length(which(colnames(conversion_table) %in% c("mm.ens", "hs.ens", "mgi",
-    "hgnc"))) > 1))
+  stopifnot(`Conversion table must be provided with at least two of column names mm.ens, hs.ens, mgi and/or hgnc` = (is(
+    conversion_table,
+    "data.frame"
+  ) & length(which(colnames(conversion_table) %in% c(
+    "mm.ens", "hs.ens", "mgi",
+    "hgnc"
+  ))) > 1))
   if (from == "ENSMUSG") {
     col1 <- conversion_table$mm.ens
   }
@@ -134,5 +97,3 @@ table_convert_genes <- function(genes, from, to, conversion_table) {
   genesV2 <- cbind(col1[which(col1 %in% genes)], col2[which(col1 %in% genes)])
   return(genesV2)
 }
-## Back up gene table for back up function: mouse_human_genes =
-## read.csv('http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt',sep='\t')