diff --git a/DESCRIPTION b/DESCRIPTION
index 9aaa7e4..db8b0c0 100644
--- a/DESCRIPTION
+++ b/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,
diff --git a/NAMESPACE b/NAMESPACE
index ac9ddec..5b41425 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/class_definitions.R b/R/class_definitions.R
index 3190181..41a19b7 100644
--- a/R/class_definitions.R
+++ b/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 = ""
+ )
}
})
diff --git a/R/convenience_fxns.R b/R/convenience_fxns.R
index d826ffd..985917d 100644
--- a/R/convenience_fxns.R
+++ b/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,58 +47,12 @@ 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)
@@ -102,17 +60,22 @@ collate_network_items <- function(dom, clusters = NULL, return = NULL) {
#' 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')
diff --git a/R/import_fxns.R b/R/import_fxns.R
index c3621a9..ae48d57 100644
--- a/R/import_fxns.R
+++ b/R/import_fxns.R
@@ -7,9 +7,9 @@
NULL
#' Create a receptor-ligand map from a cellphonedb signaling database
-#'
+#'
#' Generates a data frame of ligand-receptor interactions from a CellPhoneDB database annotating the genes encoding the interacting ligands and receptors to be queried in transcriptomic data.
-#'
+#'
#' @param genes dataframe or file path to table of gene names in uniprot, hgnc_symbol, or ensembl format in cellphonedb database format
#' @param proteins dataframe or file path to table of protein features in cellphonedb format
#' @param interactions dataframe or file path to table of protein-protein interactions in cellphonedb format
@@ -21,28 +21,41 @@ NULL
#' @param alternate_convert_table supplied table for non-ensembl method of conversion
#' @return Data frame where each row describes a possible receptor-ligand interaction
#' @export create_rl_map_cellphonedb
+#' @examples
+#' rl_map_tiny <- create_rl_map_cellphonedb(genes = domino2:::genes_tiny,
+#' proteins = domino2:::proteins_tiny, interactions = domino2:::interactions_tiny,
+#' complexes = domino2:::complexes_tiny)
#'
-create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes = NULL, database_name = "CellPhoneDB",
- gene_conv = NULL, gene_conv_host = "https://www.ensembl.org", alternate_convert = FALSE, alternate_convert_table = NULL) {
+create_rl_map_cellphonedb <- function(
+ genes, proteins, interactions, complexes = NULL, database_name = "CellPhoneDB",
+ gene_conv = NULL, gene_conv_host = "https://www.ensembl.org", alternate_convert = FALSE, alternate_convert_table = NULL) {
# Check input structures:
- stopifnot(`genes argument must be file path or dataframe` = (is(genes, "data.frame") | is(genes,
- "character")))
- stopifnot(`proteins argument must be file path or dataframe` = (is(proteins, "data.frame") | is(proteins,
- "character")))
+ stopifnot(`genes argument must be file path or dataframe` = (is(genes, "data.frame") | is(
+ genes,
+ "character"
+ )))
+ stopifnot(`proteins argument must be file path or dataframe` = (is(proteins, "data.frame") | is(
+ proteins,
+ "character"
+ )))
stopifnot(`interactions argument must be file path or dataframe` = (is(interactions, "data.frame") |
is(interactions, "character")))
- stopifnot(`complexes argument must be NULL, file path or dataframe` = (is.null(complexes) | is(complexes,
- "data.frame") | is(complexes, "character")))
+ stopifnot(`complexes argument must be NULL, file path or dataframe` = (is.null(complexes) | is(
+ complexes,
+ "data.frame"
+ ) | is(complexes, "character")))
stopifnot(`Database name must be a string` = is(database_name, "character") & length(database_name) ==
1)
stopifnot(`Gene conversion must be NULL or a character vector with 2 items` = (is.null(gene_conv) |
(is(gene_conv, "character") & length(gene_conv) == 2)))
stopifnot(`Gene conversion host must be a string` = is(gene_conv_host, "character") & length(gene_conv_host) ==
1)
- stopifnot(`Alternate conversion argument (not recommended) must be TRUE or FALSE` = is(alternate_convert,
- "logical"))
- if(alternate_convert & is.null(alternate_convert_table)) {
- stop("If using alternate conversion table (not recommended), a table must be provided")
+ stopifnot(`Alternate conversion argument (not recommended) must be TRUE or FALSE` = is(
+ alternate_convert,
+ "logical"
+ ))
+ 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:
@@ -60,32 +73,42 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
}
# replace empty cells in columns annotating gene properties with 'False' There are some
# unannotated genes in database v2.0 that seem to have been fixed in v4.0
- gene_features <- c("transmembrane", "peripheral", "secreted", "secreted_highlight", "receptor",
- "integrin", "other")
+ gene_features <- c(
+ "transmembrane", "peripheral", "secreted", "secreted_highlight", "receptor",
+ "integrin", "other"
+ )
proteins[proteins$receptor == "", colnames(proteins) %in% gene_features] <- "False"
# change cases of True/False syntax from Python to TRUE/FALSE R syntax
for (x in colnames(genes)) {
- if (identical(unique(genes[[x]]), c("True", "False")) | identical(unique(genes[[x]]), c("False",
- "True"))) {
+ if (identical(unique(genes[[x]]), c("True", "False")) | identical(unique(genes[[x]]), c(
+ "False",
+ "True"
+ ))) {
genes[[x]] <- ifelse(genes[[x]] == "True", TRUE, FALSE)
}
}
for (x in colnames(proteins)) {
- if (identical(unique(proteins[[x]]), c("True", "False")) | identical(unique(proteins[[x]]),
- c("False", "True"))) {
+ if (identical(unique(proteins[[x]]), c("True", "False")) | identical(
+ unique(proteins[[x]]),
+ c("False", "True")
+ )) {
proteins[[x]] <- ifelse(proteins[[x]] == "True", TRUE, FALSE)
}
}
for (x in colnames(interactions)) {
- if (identical(unique(interactions[[x]]), c("True", "False")) | identical(unique(interactions[[x]]),
- c("False", "True"))) {
+ if (identical(unique(interactions[[x]]), c("True", "False")) | identical(
+ unique(interactions[[x]]),
+ c("False", "True")
+ )) {
interactions[[x]] <- ifelse(interactions[[x]] == "True", TRUE, FALSE)
}
}
if (!is.null(complexes)) {
for (x in colnames(complexes)) {
- if (identical(unique(complexes[[x]]), c("True", "False")) | identical(unique(complexes[[x]]),
- c("False", "True"))) {
+ if (identical(unique(complexes[[x]]), c("True", "False")) | identical(
+ unique(complexes[[x]]),
+ c("False", "True")
+ )) {
complexes[[x]] <- ifelse(complexes[[x]] == "True", TRUE, FALSE)
}
}
@@ -94,8 +117,10 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
if (!is.null(gene_conv) & !identical(gene_conv[1], gene_conv[2])) {
# obtain conversion dictionary
if (alternate_convert) {
- conv_dict <- table_convert_genes(genes$gene_name, from = gene_conv[1], to = gene_conv[2],
- alternate_convert_table)
+ conv_dict <- table_convert_genes(genes$gene_name,
+ from = gene_conv[1], to = gene_conv[2],
+ alternate_convert_table
+ )
} else {
conv_dict <- convert_genes(genes$gene_name, from = gene_conv[1], to = gene_conv[2], host = gene_conv_host)
}
@@ -122,11 +147,11 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
# if the original gene trying to be converted is not in the gene dictionary the
# interaction is not included in the final rl_map
if (sum(g %in% conv_dict[, 1]) < length(g)) {
- for (gn in g) {
- conversion_flag[[gn]] <- TRUE
- }
+ for (gn in g) {
+ conversion_flag[[gn]] <- TRUE
+ }
} else {
- g <- paste(unique(conv_dict[conv_dict[, 1] %in% g, 2]), collapse = ";")
+ g <- paste(unique(conv_dict[conv_dict[, 1] %in% g, 2]), collapse = ";")
}
}
return(g)
@@ -146,7 +171,7 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
# interaction is not included in the final rl_map
if (sum(gene_a %in% conv_dict[, 1]) < length(gene_a)) {
for (gn in gene_a) {
- conversion_flag[[gn]] <- TRUE
+ conversion_flag[[gn]] <- TRUE
}
} else {
gene_a <- unique(conv_dict[conv_dict[, 1] %in% gene_a, 2])
@@ -160,8 +185,8 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
next
}
if (length(conversion_flag)) {
- print(paste("No gene orthologs found for:", names(conversion_flag), collapse = " "))
- print(paste("Skipping interaction:", partner_a, partner_b, collapse = " "))
+ message(paste("No gene orthologs found for:", names(conversion_flag), collapse = " "))
+ message(paste("Skipping interaction:", partner_a, partner_b, collapse = " "))
next
}
a_df <- as.data.frame(a_features)
@@ -178,11 +203,11 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
# if the original gene trying to be converted is not in the gene dictionary the
# interaction is not included in the final rl_map
if (sum(g %in% conv_dict[, 1]) < length(g)) {
- for (gn in g) {
- conversion_flag[[gn]] <- TRUE
- }
+ for (gn in g) {
+ conversion_flag[[gn]] <- TRUE
+ }
} else {
- g <- paste(unique(conv_dict[conv_dict[, 1] %in% g, 2]), collapse = ";")
+ g <- paste(unique(conv_dict[conv_dict[, 1] %in% g, 2]), collapse = ";")
}
}
return(g)
@@ -202,7 +227,7 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
# interaction is not included in the final rl_map
if (sum(gene_b %in% conv_dict[, 1]) < length(gene_b)) {
for (gn in gene_b) {
- conversion_flag[[gn]] <- TRUE
+ conversion_flag[[gn]] <- TRUE
}
} else {
gene_b <- unique(conv_dict[conv_dict[, 1] %in% gene_b, 2])
@@ -216,8 +241,8 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
next
}
if (length(conversion_flag)) {
- print(paste("No gene orthologs found for:", names(conversion_flag), collapse = " "))
- print(paste("Skipping interaction:", partner_a, partner_b, collapse = " "))
+ message(paste("No gene orthologs found for:", names(conversion_flag), collapse = " "))
+ message(paste("Skipping interaction:", partner_a, partner_b, collapse = " "))
next
}
b_df <- as.data.frame(b_features)
@@ -232,21 +257,25 @@ create_rl_map_cellphonedb <- function(genes, proteins, interactions, complexes =
rl_map <- rl_map[!(rl_map$type_A == "R" & rl_map$type_B == "R") & !(rl_map$type_A == "L" & rl_map$type_B ==
"L"), ]
# specify column order
- rl_map <- rl_map[, c("int_pair", "name_A", "uniprot_A", "gene_A", "type_A", "name_B", "uniprot_B",
- "gene_B", "type_B", "annotation_strategy", "source", "database_name")]
+ rl_map <- rl_map[, c(
+ "int_pair", "name_A", "uniprot_A", "gene_A", "type_A", "name_B", "uniprot_B",
+ "gene_B", "type_B", "annotation_strategy", "source", "database_name"
+ )]
return(rl_map)
}
#' Create a list of genes in regulons inferred by SCENIC
-#'
+#'
#' Generates a list of transcription factors and the genes targeted by the transcription factor as part of their regulon inferred by pySCENIC
-#'
+#'
#' @param regulons Dataframe or file path to the table of the output of the grn (gene regulatory network) function from pySCENIC
#' @return A list where names are transcription factors and the stored values are character vectors of genes in the inferred regulons
#' @export create_regulon_list_scenic
-#'
+#' @examples
+#' regulon_list_tiny <- create_regulon_list_scenic(regulons = domino2:::regulons_tiny)
+#'
create_regulon_list_scenic <- function(regulons) {
- if (is(regulons, "character")){
+ if (is(regulons, "character")) {
regulons <- read.csv(regulons)
}
TFS <- unique(regulons[["TF"]])
@@ -270,15 +299,15 @@ create_regulon_list_scenic <- function(regulons) {
}
#' Create a domino object and prepare it for network construction
-#'
-#' 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 rl_map Data frame where each row describes a receptor-ligand interaction with required columns gene_A & gene_B including the gene names for the receptor and ligand and type_A & type_B annotating if genes A and B are a ligand (L) or receptor (R)
#' @param features Either a path to a csv containing cell level features of interest (ie. the auc matrix from pySCENIC) or named matrix with cells as columns and features as rows.
#' @param ser Seurat object containing scaled RNA expression data in the RNA assay slot and cluster identity. Either a ser object OR z_scores and clusters must be provided. If ser is present z_scores and clusters will be ignored.
@@ -295,21 +324,39 @@ create_regulon_list_scenic <- function(regulons) {
#' @param tf_variance_quantile What proportion of variable features to take if using variance to threshold features. Default is 0.5. Higher numbers will keep more features. Ignored if tf_selection_method is not 'variable'
#' @return A domino object
#' @export create_domino
-#'
-create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores = NULL, clusters = NULL,
- use_clusters = TRUE, tf_targets = NULL, verbose = TRUE, use_complexes = TRUE, rec_min_thresh = 0.025,
- remove_rec_dropout = TRUE, tf_selection_method = "clusters", tf_variance_quantile = 0.5) {
+#' @examples
+#' pbmc_dom_tiny_all <- pbmc_dom_tiny <- create_domino(rl_map = domino2:::rl_map_tiny,
+#' features = domino2:::auc_tiny, counts = domino2:::RNA_count_tiny, z_scores = domino2:::RNA_zscore_tiny,
+#' clusters = domino2:::clusters_tiny, tf_targets = domino2:::regulon_list_tiny, use_clusters = FALSE,
+#' use_complexes = FALSE, rec_min_thresh = 0.1, remove_rec_dropout = TRUE,
+#' tf_selection_method = "all")
+#'
+#' pbmc_dom_tiny_clustered <- create_domino(rl_map = domino2:::rl_map_tiny,
+#' features = domino2:::auc_tiny, counts = domino2:::RNA_count_tiny, z_scores = domino2:::RNA_zscore_tiny,
+#' clusters = domino2:::clusters_tiny, tf_targets = domino2:::regulon_list_tiny,
+#' use_clusters = TRUE, use_complexes = TRUE, remove_rec_dropout = FALSE)
+#'
+create_domino <- function(
+ rl_map, features, ser = NULL, counts = NULL, z_scores = NULL, clusters = NULL,
+ use_clusters = TRUE, tf_targets = NULL, verbose = TRUE, use_complexes = TRUE, rec_min_thresh = 0.025,
+ remove_rec_dropout = TRUE, tf_selection_method = "clusters", tf_variance_quantile = 0.5) {
# Check inputs:
- stopifnot(`rl_map must be a data.frame with column names gene_A, gene_B, type_A, and type_B` = (is(rl_map,
- "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)) &
+ stopifnot(`rl_map must be a data.frame with column names gene_A, gene_B, type_A, and type_B` = (is(
+ rl_map,
+ "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 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)))))
+ 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))
# Create object
@@ -424,19 +471,20 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
rownames(p_vals) <- rownames(features)
colnames(p_vals) <- levels(dom@clusters)
if (verbose) {
- print("Calculating feature enrichment by cluster")
+ message("Calculating feature enrichment by cluster")
clust_n <- length(levels(dom@clusters))
}
for (clust in levels(dom@clusters)) {
if (verbose) {
cur <- which(levels(dom@clusters) == clust)
- print(paste0(cur, " of ", clust_n))
+ message(paste0(cur, " of ", clust_n))
}
cells <- which(dom@clusters == clust)
for (feat in rownames(dom@features)) {
p_vals[feat, clust] <- stats::wilcox.test(
- dom@features[feat, cells], dom@features[feat, -cells], alternative = "g"
- )$p.value
+ dom@features[feat, cells], dom@features[feat, -cells],
+ alternative = "g"
+ )$p.value
}
}
dom@clust_de <- p_vals
@@ -447,7 +495,7 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
if (tf_selection_method == "variable") {
dom@clusters <- factor()
variances <- apply(dom@features, 1, function(x) {
- sd(x)/mean(x)
+ sd(x) / mean(x)
})
keep_n <- length(variances) * tf_variance_quantile
keep_id <- which(rank(variances) > keep_n)
@@ -459,7 +507,7 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
warning("tf_targets is not a list. No regulons stored")
} else {
dom@linkages[["tf_targets"]] <- tf_targets
- }
+ }
# Calculate correlation matrix between features and receptors.
dom@counts <- counts
zero_sum <- rowSums(counts == 0)
@@ -469,7 +517,7 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
rownames(rho) <- ser_receptors
colnames(rho) <- rownames(dom@features)
if (verbose) {
- print("Calculating correlations")
+ message("Calculating correlations")
n_tf <- nrow(dom@features)
}
for (module in rownames(dom@features)) {
@@ -477,10 +525,10 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
# correlation equal to 0.
if (verbose) {
cur <- which(rownames(dom@features) == module)
- print(paste0(cur, " of ", n_tf))
+ message(paste0(cur, " of ", n_tf))
}
if (!is.null(dom@linkages$tf_targets)) {
- tf <- gsub(pattern = "\\.\\.\\.", replacement = "", module) # correction for AUC values from pySCENIC that append an elipses to TF names due to (+) characters in the orignial python output
+ tf <- gsub(pattern = "\\.\\.\\.", replacement = "", module) # correction for AUC values from pySCENIC that append an elipses to TF names due to (+) characters in the orignial python output
module_targets <- tf_targets[[tf]]
module_rec_targets <- intersect(module_targets, ser_receptors)
} else {
@@ -524,7 +572,7 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
r_genes <- r
}
if (sum(rownames(rho) %in% r_genes) != length(r_genes)) {
- print(paste0(r, " has component genes that did not pass testing parameters"))
+ message(paste0(r, " has component genes that did not pass testing parameters"))
cor_list[[r]] <- rep(0, ncol(rho))
next
}
@@ -547,7 +595,7 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
for (rec in ser_receptors) {
rec_percent <- sapply(X = levels(dom@clusters), FUN = function(x) {
# percentage of cells in cluster with non-zero expression of receptor gene
- sum(dom@counts[rec, dom@clusters == x] > 0)/length(dom@counts[rec, dom@clusters ==
+ sum(dom@counts[rec, dom@clusters == x] > 0) / length(dom@counts[rec, dom@clusters ==
x])
})
cl_rec_percent <- rbind(cl_rec_percent, rec_percent)
@@ -559,23 +607,25 @@ create_domino <- function(rl_map, features, ser = NULL, counts = NULL, z_scores
}
#' Use biomaRt to convert genes
-#'
+#'
#' This function reads in a vector of genes and converts the genes to specified symbol type
-#'
+#'
#' @param genes Vector of genes to convert.
#' @param from Format of gene input (ENSMUSG, ENSG, MGI, or HGNC)
#' @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"),
- host = "https://www.ensembl.org") {
+#'
+convert_genes <- function(
+ genes, from = c("ENSMUSG", "ENSG", "MGI", "HGNC"), to = c("MGI", "HGNC"),
+ host = "https://www.ensembl.org") {
# 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(`Host must be web host to connect to` = (is(host, "character") & length(host) == 1))
if (from == "ENSMUSG") {
@@ -602,21 +652,29 @@ convert_genes <- function(genes, from = c("ENSMUSG", "ENSG", "MGI", "HGNC"), to
tarMart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl", host = host)
tarAtts <- "hgnc_symbol"
}
- genesV2 <- getLDS(attributes = sourceAtts, filters = sourceAtts, values = genes, mart = srcMart,
- attributesL = tarAtts, martL = tarMart, uniqueRows = FALSE)
+ genesV2 <- getLDS(
+ attributes = sourceAtts, filters = sourceAtts, values = genes, mart = srcMart,
+ attributesL = tarAtts, martL = tarMart, uniqueRows = FALSE
+ )
return(genesV2)
}
+
#' Adds a column to the RL signaling data frame.
-#'
-#' This function adds a column to the internal rl 'map' used to map all
+#'
+#' This function adds a column to the internal rl 'map' used to map all
#' receptor and receptor complexes to all ligand and ligand complexes.
-#'
+#'
#' @param map RL signaling data frame.
#' @param map_ref Name of column to match new data to
#' @param conv Data frame matching current data in map to new data.
#' @param new_name Name of new column to be created in RL map
#' @return An updated RL signaling data frame
-#'
+#' @export
+#' @examples
+#' lr_name <- data.frame("abbrev" = c("L", "R"), "full" = c("Ligand", "Receptor"))
+#' rl_map_expanded <- add_rl_column(map = domino2:::rl_map_tiny, map_ref = "type_A",
+#' conv = lr_name, new_name = "type_A_full")
+#'
add_rl_column <- function(map, map_ref, conv, new_name) {
map_in_ref <- match(map[[map_ref]], conv[, 1])
not_in_ref <- which(is.na(map_in_ref))
@@ -643,10 +701,10 @@ add_rl_column <- function(map, map_ref, conv, new_name) {
}
#' Calculate mean ligand expression as a data.frame for plotting in circos plot
-#'
+#'
#' Creates a data frame of mean ligand expression for use in plotting a circos
#' plot of ligand expression and saving tables of mean expression.
-#'
+#'
#' @param x Gene by cell expression matrix
#' @param ligands Character vector of ligand genes to be quantified
#' @param cell_ident Vector of cell type (identity) names for which to calculate mean ligand gene expression
@@ -654,6 +712,11 @@ 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
+#' @examples
+#' counts <- dom_counts(domino2:::pbmc_dom_built_tiny)
+#' mean_exp <- mean_ligand_expression(counts,
+#' ligands = c("PTPRC", "FASLG"), cell_ident = "CD14_monocyte",
+#' cell_barcodes = colnames(counts), destination = "FAS")
#'
mean_ligand_expression <- function(x, ligands, cell_ident, cell_barcodes, destination){
# initiate data frame to store results
diff --git a/R/plot_fxns.R b/R/plot_fxns.R
index 8e17c42..ad2d731 100644
--- a/R/plot_fxns.R
+++ b/R/plot_fxns.R
@@ -8,10 +8,10 @@
NULL
#' Create a network heatmap
-#'
-#' Creates a heatmap of the signaling network. Alternatively, the network
+#'
+#' Creates a heatmap of the signaling network. Alternatively, the network
#' matrix can be accessed directly in the signaling slot of a domino object.
-#'
+#'
#' @param dom Domino object with network built ([build_domino()])
#' @param clusts Vector of clusters to be included. If NULL then all clusters are used.
#' @param min_thresh Minimum signaling threshold for plotting. Defaults to -Inf for no threshold.
@@ -29,8 +29,9 @@ NULL
#' #normalize
#' signaling_heatmap(domino2:::pbmc_dom_built_tiny, normalize = "rec_norm")
#'
-signaling_heatmap <- function(dom, clusts = NULL, min_thresh = -Inf, max_thresh = Inf, scale = "none",
- normalize = "none", ...) {
+signaling_heatmap <- function(
+ dom, clusts = NULL, min_thresh = -Inf, max_thresh = Inf, scale = "none",
+ normalize = "none", ...) {
if (!dom@misc[["build"]]) {
stop("Please run domino_build prior to generate signaling network.")
}
@@ -63,15 +64,16 @@ signaling_heatmap <- function(dom, clusts = NULL, min_thresh = -Inf, max_thresh
...
)
}
+
#' Create a cluster incoming signaling heatmap
-#'
+#'
#' Creates a heatmap of a cluster incoming signaling matrix. Each cluster has a
#' list of ligands capable of activating its enriched transcription factors. The
-#' function creates a heatmap of cluster average expression for all of those
+#' function creates a heatmap of cluster average expression for all of those
#' ligands. A list of all cluster incoming signaling matrices can be found in
#' the cl_signaling_matrices slot of a domino option as an alternative to this
#' plotting function.
-#'
+#'
#' @param dom Domino object with network built ([build_domino()])
#' @param rec_clust Which cluster to select as the receptor. Must match naming of clusters in the domino object.
#' @param clusts Vector of clusters to be included. If NULL then all clusters are used.
@@ -86,9 +88,10 @@ signaling_heatmap <- function(dom, clusts = NULL, min_thresh = -Inf, max_thresh
#' @examples
#' #incoming signaling of the CD8 T cells
#' incoming_signaling_heatmap(domino2:::pbmc_dom_built_tiny, "CD8_T_cell")
-#'
-incoming_signaling_heatmap <- function(dom, rec_clust, clusts = NULL, min_thresh = -Inf, max_thresh = Inf,
- scale = "none", normalize = "none", title = TRUE, ...) {
+#'
+incoming_signaling_heatmap <- function(
+ dom, rec_clust, clusts = NULL, min_thresh = -Inf, max_thresh = Inf,
+ scale = "none", normalize = "none", title = TRUE, ...) {
if (!dom@misc[["build"]]) {
stop("Please run domino_build prior to generate signaling network.")
}
@@ -97,7 +100,7 @@ incoming_signaling_heatmap <- function(dom, rec_clust, clusts = NULL, min_thresh
}
mat <- dom@cl_signaling_matrices[[rec_clust]]
if (dim(mat)[1] == 0) {
- print("No signaling found for this cluster under build parameters.")
+ message("No signaling found for this cluster under build parameters.")
return()
}
if (!is.null(clusts)) {
@@ -151,12 +154,13 @@ incoming_signaling_heatmap <- function(dom, rec_clust, clusts = NULL, min_thresh
)
}
}
+
#' Create a cluster to cluster signaling network diagram
-#'
+#'
#' Creates a network diagram of signaling between clusters. Nodes are clusters
#' and directed edges indicate signaling from one cluster to another. Edges are
#' colored based on the color scheme of the ligand expressing cluster.
-#'
+#'
#' @param dom Domino object with network built ([build_domino()])
#' @param cols Named vector indicating the colors for clusters. Values are colors and names must match clusters in the domino object. If left as NULL then ggplot colors are generated for the clusters.
#' @param edge_weight Weight for determining thickness of edges on plot. Signaling values are multiplied by this value.
@@ -165,9 +169,9 @@ incoming_signaling_heatmap <- function(dom, rec_clust, clusts = NULL, min_thresh
#' @param showIncomingSignalingClusts Vector of clusters to plot the incoming signaling on
#' @param min_thresh Minimum signaling threshold. Values lower than the threshold will be set to the threshold. Defaults to -Inf for no threshold.
#' @param max_thresh Maximum signaling threshold for plotting. Values higher than the threshold will be set to the threshold. Defaults to Inf for no threshold.
-#' @param normalize Options to normalize the signaling 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 normalize Options to normalize the signaling 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 scale How to scale the values (after thresholding). Options are 'none', 'sqrt' for square root, 'log' for log10, or 'sq' for square.
-#' @param layout Type of layout to use. Options are 'random', 'sphere', 'circle', 'fr' for Fruchterman-Reingold force directed layout, and 'kk' for Kamada Kawai for directed layout.
+#' @param layout Type of layout to use. Options are 'random', 'sphere', 'circle', 'fr' for Fruchterman-Reingold force directed layout, and 'kk' for Kamada Kawai for directed layout.
#' @param scale_by How to size vertices. Options are 'lig_sig' for summed outgoing signaling, 'rec_sig' for summed incoming signaling, and 'none'. In the former two cases the values are scaled with asinh after summing all incoming or outgoing signaling.
#' @param vert_scale Integer used to scale size of vertices with our without variable scaling from size_verts_by.
#' @param plot_title Text for the plot's title.
@@ -181,9 +185,10 @@ incoming_signaling_heatmap <- function(dom, rec_clust, clusts = NULL, min_thresh
#' signaling_network(domino2:::pbmc_dom_built_tiny, showOutgoingSignalingClusts = "CD14_monocyte", scale = "none",
#' norm = "none", layout = "fr", scale_by = "none", vert_scale = 5)
#'
-signaling_network <- function(dom, cols = NULL, edge_weight = 0.3, clusts = NULL, showOutgoingSignalingClusts = NULL,
- showIncomingSignalingClusts = NULL, min_thresh = -Inf, max_thresh = Inf, normalize = "none", scale = "sq",
- layout = "circle", scale_by = "rec_sig", vert_scale = 3, plot_title = NULL, ...) {
+signaling_network <- function(
+ dom, cols = NULL, edge_weight = 0.3, clusts = NULL, showOutgoingSignalingClusts = NULL,
+ showIncomingSignalingClusts = NULL, min_thresh = -Inf, max_thresh = Inf, normalize = "none", scale = "sq",
+ layout = "circle", scale_by = "rec_sig", vert_scale = 3, plot_title = NULL, ...) {
if (!length(dom@clusters)) {
stop("This domino object was not built with clusters so there is no intercluster signaling.")
}
@@ -263,7 +268,7 @@ signaling_network <- function(dom, cols = NULL, edge_weight = 0.3, clusts = NULL
# Get vert angle for labeling circos plot
if (layout == "circle") {
v_angles <- 1:length(igraph::V(graph))
- v_angles <- -2 * pi * (v_angles - 1)/length(v_angles)
+ v_angles <- -2 * pi * (v_angles - 1) / length(v_angles)
igraph::V(graph)$label.degree <- v_angles
}
names(v_cols) <- c()
@@ -295,17 +300,18 @@ signaling_network <- function(dom, cols = NULL, edge_weight = 0.3, clusts = NULL
}
plot(graph, layout = l, main = plot_title, ...)
}
+
#' Create a gene association network
-#'
-#' Create a gene association network for genes from a given cluster. The
+#'
+#' Create a gene association network for genes from a given cluster. The
#' selected cluster acts as the receptor for the gene association network, so
#' only ligands, receptors, and features associated with the receptor cluster
#' will be included in the plot.
-#'
+#'
#' @param dom Domino object with network built ([build_domino()])
#' @param clust Receptor cluster to create the gene association network for. A vector of clusters may be provided.
#' @param OutgoingSignalingClust Vector of clusters to plot the outgoing signaling from
-#' @param class_cols Named vector of colors used to color classes of vertices. Values must be colors and names must be classes ('rec', 'lig', and 'feat' for receptors, ligands, and features.).
+#' @param class_cols Named vector of colors used to color classes of vertices. Values must be colors and names must be classes ('rec', 'lig', and 'feat' for receptors, ligands, and features.).
#' @param cols Named vector of colors for individual genes. Genes not included in this vector will be colored according to class_cols.
#' @param lig_scale FALSE or a numeric value to scale the size of ligand vertices based on z-scored expression in the data set.
#' @param layout Type of layout to use. Options are 'grid', 'random', 'sphere', 'circle', 'fr' for Fruchterman-Reingold force directed layout, and 'kk' for Kamada Kawai for directed layout.
@@ -316,8 +322,9 @@ signaling_network <- function(dom, cols = NULL, edge_weight = 0.3, clusts = NULL
#' #basic usage
#' gene_network(domino2:::pbmc_dom_built_tiny, clust = "CD8_T_cell", OutgoingSignalingClust = "CD14_monocyte")
#'
-gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, class_cols = c(lig = "#FF685F",
- rec = "#47a7ff", feat = "#39C740"), cols = NULL, lig_scale = 1, layout = "grid", ...) {
+gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL,
+ class_cols = c(lig = "#FF685F",rec = "#47a7ff", feat = "#39C740"),
+ cols = NULL, lig_scale = 1, layout = "grid", ...) {
if (!dom@misc[["build"]]) {
warning("Please build a signaling network with domino_build prior to plotting.")
}
@@ -334,7 +341,7 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, class
# Check if signaling exists for target cluster
mat <- dom@cl_signaling_matrices[[cl]]
if (dim(mat)[1] == 0) {
- print(paste("No signaling found for", cl, "under build parameters."))
+ message(paste("No signaling found for", cl, "under build parameters."))
(next)()
}
all_sums <- c(all_sums, rowSums(mat))
@@ -344,7 +351,7 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, class
all_sums <- all_sums[!duplicated(names(all_sums))]
# If no signaling for target clusters then don't do anything
if (length(tfs) == 0) {
- print("No signaling found for provided clusters")
+ message("No signaling found for provided clusters")
return()
}
} else {
@@ -378,7 +385,7 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, class
allowed_ligs <- names(mat[mat > 0])
all_sums <- mat[mat > 0]
} else {
- allowed_ligs <- rownames(mat[rowSums(mat) > 0, ]) #I remove any ligands with zeroes for all clusters
+ allowed_ligs <- rownames(mat[rowSums(mat) > 0, ]) # I remove any ligands with zeroes for all clusters
all_sums <- rowSums(mat[rowSums(mat) > 0, ])
}
} else {
@@ -434,9 +441,9 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, class
l[all_ligs, 1] <- -0.75
l[all_recs, 1] <- 0
l[all_tfs, 1] <- 0.75
- l[all_ligs, 2] <- (1:length(all_ligs)/mean(1:length(all_ligs)) - 1) * 2
- l[all_recs, 2] <- (1:length(all_recs)/mean(1:length(all_recs)) - 1) * 2
- l[all_tfs, 2] <- (1:length(all_tfs)/mean(1:length(all_tfs)) - 1) * 2
+ l[all_ligs, 2] <- (1:length(all_ligs) / mean(1:length(all_ligs)) - 1) * 2
+ l[all_recs, 2] <- (1:length(all_recs) / mean(1:length(all_recs)) - 1) * 2
+ l[all_tfs, 2] <- (1:length(all_tfs) / mean(1:length(all_tfs)) - 1) * 2
rownames(l) <- c()
} else if (layout == "random") {
l <- igraph::layout_randomly(graph)
@@ -452,11 +459,12 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, class
plot(graph, layout = l, main = paste0("Signaling ", OutgoingSignalingClust, " to ", clust), ...)
return(invisible(list(graph = graph, layout = l)))
}
+
#' Create a heatmap of features organized by cluster
-#'
+#'
#' Creates a heatmap of feature expression (typically transcription factor
#' activation scores) by cells organized by cluster.
-#'
+#'
#' @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.
@@ -476,8 +484,9 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, class
#' #using thresholds
#' feat_heatmap(domino2:::pbmc_dom_built_tiny, min_thresh = 0.1, max_thresh = 0.6, norm = TRUE, bool = FALSE)
#'
-feat_heatmap <- function(dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, title = TRUE, norm = FALSE,
- cols = NULL, ann_cols = TRUE, min_thresh = NULL, max_thresh = NULL, ...) {
+feat_heatmap <- function(
+ dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, title = TRUE, norm = FALSE,
+ cols = NULL, ann_cols = TRUE, min_thresh = NULL, max_thresh = NULL, ...) {
if (!length(dom@clusters)) {
warning("This domino object wasn't built with clusters. Cells will not be ordered.")
ann_cols <- FALSE
@@ -518,7 +527,7 @@ feat_heatmap <- function(dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, tit
na <- which(is.na(mid))
na_feats <- paste(feats[na], collapse = " ")
if (length(na) != 0) {
- print(paste("Unable to find", na_feats))
+ message(paste("Unable to find", na_feats))
feats <- feats[-na]
}
} else if (feats == "all") {
@@ -575,11 +584,12 @@ feat_heatmap <- function(dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, tit
)
}
}
+
#' Create a heatmap of correlation between receptors and transcription factors
-#'
-#' Creates a heatmap of correlation values between receptors and transcription
+#'
+#' Creates a heatmap of correlation values between receptors and transcription
#' factors either with boolean threshold or with continuous values displayed
-#'
+#'
#' @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.
@@ -597,9 +607,10 @@ feat_heatmap <- function(dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, tit
#' cor_heatmap(domino2:::pbmc_dom_built_tiny, bool = TRUE, bool_thresh = 0.25)
#' #identify combinations that are connected
#' cor_heatmap(domino2:::pbmc_dom_built_tiny, bool = FALSE, mark_connections = TRUE)
-#'
-cor_heatmap <- function(dom, bool = FALSE, bool_thresh = 0.15, title = TRUE, feats = NULL, recs = NULL,
- mark_connections = FALSE, ...) {
+#'
+cor_heatmap <- function(
+ dom, bool = FALSE, bool_thresh = 0.15, title = TRUE, feats = NULL, recs = NULL,
+ mark_connections = FALSE, ...) {
mat <- dom@cor
if (bool) {
cp <- mat
@@ -622,7 +633,7 @@ cor_heatmap <- function(dom, bool = FALSE, bool_thresh = 0.15, title = TRUE, fea
na <- which(is.na(mid))
na_feats <- paste(feats[na], collapse = " ")
if (length(na) != 0) {
- print(paste("Unable to find", na_feats))
+ message(paste("Unable to find", na_feats))
feats <- feats[-na]
}
} else if (identical(feats, "all")) {
@@ -673,10 +684,11 @@ cor_heatmap <- function(dom, bool = FALSE, bool_thresh = 0.15, title = TRUE, fea
)
}
}
+
#' Create a correlation plot between transcription factor activation score and receptor
-#'
+#'
#' Create a correlation plot between transcription factor activation score and receptor
-#'
+#'
#' @param dom Domino object with network built ([build_domino()])
#' @param tf Target TF module for plotting with receptor
#' @param rec Target receptor for plotting with TF
@@ -684,10 +696,9 @@ cor_heatmap <- function(dom, bool = FALSE, bool_thresh = 0.15, title = TRUE, fea
#' @param ... Other parameters to pass to ggscatter.
#' @return a ggplot object
#' @export cor_scatter
-#' @examples
-#' #basic usage
+#' @examples
#' cor_scatter(domino2:::pbmc_dom_built_tiny, "ATF4","CD22")
-#'
+#'
cor_scatter <- function(dom, tf, rec, remove_rec_dropout = TRUE, ...) {
if (remove_rec_dropout) {
keep_id <- which(dom@counts[rec, ] > 0)
@@ -701,11 +712,12 @@ cor_scatter <- function(dom, tf, rec, remove_rec_dropout = TRUE, ...) {
ggscatter(dat, x = "rec", y = "tf", add = "reg.line", conf.int = FALSE, cor.coef = FALSE,
cor.method = "pearson", xlab = rec, ylab = tf, size = 0.25, ...)
}
+
#' Plot expression of a receptor's ligands by other cell types as a chord plot
-#'
+#'
#' Creates a chord plot of expression of ligands that can activate a specified
#' receptor where chord widths correspond to mean ligand expression by the cluster.
-#'
+#'
#' @param dom Domino object that has undergone network building with [build_domino()]
#' @param receptor Name of a receptor active in at least one cell type in the domino object
#' @param ligand_expression_threshold Minimum mean expression value of a ligand by a cell type for a chord to be rendered between the cell type and the receptor
@@ -721,8 +733,9 @@ cor_scatter <- function(dom, tf, rec, remove_rec_dropout = TRUE, ...) {
#' names(cols) = levels(domino2:::pbmc_dom_built_tiny@clusters)
#' circos_ligand_receptor(domino2:::pbmc_dom_built_tiny, receptor = "FAS", cell_colors = cols)
#'
-circos_ligand_receptor <- function(dom, receptor, ligand_expression_threshold = 0.01, cell_idents = NULL,
- cell_colors = NULL) {
+circos_ligand_receptor <- function(
+ dom, receptor, ligand_expression_threshold = 0.01, cell_idents = NULL,
+ cell_colors = NULL) {
ligands <- dom@linkages$rec_lig[[receptor]]
signaling_df <- NULL
if (is.null(cell_idents)) {
@@ -749,7 +762,7 @@ circos_ligand_receptor <- function(dom, receptor, ligand_expression_threshold =
signaling_df$mean.expression[is.na(signaling_df$mean.expression)] <- 0
# create a scaled mean expression plot for coord widths greater than 1 by dividing by the max
# expression [range (0-1)] scaled.mean will only be used when the max expression is > 1
- signaling_df$scaled.mean.expression <- signaling_df$mean.expression/max(signaling_df$mean.expression)
+ signaling_df$scaled.mean.expression <- signaling_df$mean.expression / max(signaling_df$mean.expression)
# exit function if no ligands are expressed above ligand expression threshold
if (sum(signaling_df[["mean.expression"]] > ligand_expression_threshold) == 0) {
stop(paste0("No ligands of ", receptor, " exceed ligand expression threshold."))
@@ -758,13 +771,14 @@ circos_ligand_receptor <- function(dom, receptor, ligand_expression_threshold =
arc_df <- signaling_df[, c("origin", "destination")]
arc_df["ligand.arc"] <- 1
# receptor arc will always sum to 4 no matter how many ligands and cell idents are plotted
- arc_df["receptor.arc"] <- 4/(nrow(signaling_df))
+ arc_df["receptor.arc"] <- 4 / (nrow(signaling_df))
# name grouping based on [cell_ident]
nm <- c(receptor, arc_df$origin)
group <- structure(c(nm[1], gsub("-.*", "", nm[-1])), names = nm)
# order group as a factor with the receptor coming first
- group <- factor(group, levels = c(receptor, sort(unique(gsub("-.*", "", nm))[-1]) # alphabetical order of the other cell idents
-))
+ group <- factor(group, levels = c(
+ receptor, sort(unique(gsub("-.*", "", nm))[-1]) # alphabetical order of the other cell idents
+ ))
# colors for ligand chords
lig_colors <- ggplot_col_gen(length(ligands))
names(lig_colors) <- ligands
@@ -773,15 +787,17 @@ circos_ligand_receptor <- function(dom, receptor, ligand_expression_threshold =
cell_colors <- ggplot_col_gen(length(cell_idents))
names(cell_colors) <- cell_idents
}
- grid_col <- c("#FFFFFF") # hide the arc corresponding to the receptor by coloring white
+ grid_col <- c("#FFFFFF") # hide the arc corresponding to the receptor by coloring white
for (i in 1:length(ligands)) {
grid_col <- c(grid_col, rep(lig_colors[i], length(cell_idents)))
}
names(grid_col) <- c(receptor, signaling_df$origin)
circlize::circos.clear()
circlize::circos.par(start.degree = 0)
- circlize::chordDiagram(arc_df, group = group, grid.col = grid_col, link.visible = FALSE, annotationTrack = c("grid"),
- preAllocateTracks = list(track.height = circlize::mm_h(4), track.margin = c(circlize::mm_h(2), 0)), big.gap = 2)
+ circlize::chordDiagram(arc_df,
+ group = group, grid.col = grid_col, link.visible = FALSE, annotationTrack = c("grid"),
+ preAllocateTracks = list(track.height = circlize::mm_h(4), track.margin = c(circlize::mm_h(2), 0)), big.gap = 2
+ )
for (send in signaling_df$origin) {
if (signaling_df[signaling_df$origin == send, ][["mean.expression"]] > ligand_expression_threshold) {
if (max(signaling_df[["mean.expression"]]) > 1) {
@@ -791,8 +807,10 @@ circos_ligand_receptor <- function(dom, receptor, ligand_expression_threshold =
expr <- signaling_df[signaling_df$origin == send, ][["mean.expression"]]
max_width <- 1
}
- circlize::circos.link(send, c(0.5 - (expr/2), 0.5 + (expr/2)), receptor, 2, col = paste0(grid_col[[send]],
- "88"))
+ circlize::circos.link(send, c(0.5 - (expr / 2), 0.5 + (expr / 2)), receptor, 2, col = paste0(
+ grid_col[[send]],
+ "88"
+ ))
}
}
sector_names <- circlize::get.all.sector.index()
@@ -800,41 +818,55 @@ circos_ligand_receptor <- function(dom, receptor, ligand_expression_threshold =
for (cell in cell_sectors) {
row_pick <- sector_names[grepl(paste0("^", cell), sector_names)]
if (length(row_pick)) {
- circlize::highlight.sector(sector_names[grepl(paste0("^", cell, "-"), sector_names)], track.index = 1,
+ circlize::highlight.sector(sector_names[grepl(paste0("^", cell, "-"), sector_names)],
+ track.index = 1,
col = cell_colors[cell], text = cell, cex = 1, facing = "inside", text.col = "black",
- niceFacing = FALSE, text.vjust = -1.5)
+ niceFacing = FALSE, text.vjust = -1.5
+ )
}
}
# highlight receptor sector
- circlize::highlight.sector(sector_names[grepl(paste0("^", receptor, "$"), sector_names)], track.index = 1,
+ circlize::highlight.sector(sector_names[grepl(paste0("^", receptor, "$"), sector_names)],
+ track.index = 1,
col = "#FFFFFF", text = receptor, cex = 1.5, facing = "clockwise", text.col = "black", niceFacing = TRUE,
- pos = 4)
+ pos = 4
+ )
# create legends
- lgd_cells <- ComplexHeatmap::Legend(at = as.character(cell_idents), type = "grid", legend_gp = grid::gpar(fill = cell_colors),
- title_position = "topleft", title = "cell identity")
- lgd_ligands <- ComplexHeatmap::Legend(at = ligands, type = "grid", legend_gp = grid::gpar(fill = lig_colors), title_position = "topleft",
- title = "ligand")
- chord_width <- 10/(4 + length(cell_idents) * length(ligands))
- lgd_chord <- ComplexHeatmap::Legend(at = c(ligand_expression_threshold, max_width), col_fun = circlize::colorRamp2(c(ligand_expression_threshold,
- max_width), c("#DDDDDD", "#DDDDDD")), legend_height = grid::unit(chord_width, "in"), title_position = "topleft",
- title = "ligand expression")
+ lgd_cells <- ComplexHeatmap::Legend(
+ at = as.character(cell_idents), type = "grid", legend_gp = grid::gpar(fill = cell_colors),
+ title_position = "topleft", title = "cell identity"
+ )
+ lgd_ligands <- ComplexHeatmap::Legend(
+ at = ligands, type = "grid", legend_gp = grid::gpar(fill = lig_colors), title_position = "topleft",
+ title = "ligand"
+ )
+ chord_width <- 10 / (4 + length(cell_idents) * length(ligands))
+ lgd_chord <- ComplexHeatmap::Legend(
+ at = c(ligand_expression_threshold, max_width), col_fun = circlize::colorRamp2(c(
+ ligand_expression_threshold,
+ max_width
+ ), c("#DDDDDD", "#DDDDDD")), legend_height = grid::unit(chord_width, "in"), title_position = "topleft",
+ title = "ligand expression"
+ )
lgd_list_vertical <- ComplexHeatmap::packLegend(lgd_cells, lgd_ligands, lgd_chord)
ComplexHeatmap::draw(lgd_list_vertical, x = grid::unit(0.02, "npc"), y = grid::unit(0.98, "npc"), just = c("left", "top"))
}
+
#' Plot differential linkages among domino results ranked by a comparative statistic
-#'
+#'
#' Plot differential linkages among domino results ranked by a comparative statistic
-#'
+#'
#' @param differential_linkages a data.frame output from the test_differential_linkages function
#' @param test_statistic column name of differential_linkages where the test statistic used for ranking linkages is stored (ex. 'p.value')
#' @param stat_range a two value vector of the minimum and maximum values of test_statistic for plotting linkage features
-#' @param stat_ranking 'ascending' (lowest value of test statisic is colored red and plotted at the top) or 'descending' (highest value of test statistic is colored red and plotted at the top).
+#' @param stat_ranking 'ascending' (lowest value of test statisic is colored red and plotted at the top) or 'descending' (highest value of test statistic is colored red and plotted at the top).
#' @param group_palette a named vector of colors to use for each group being compared
#' @return a Heatmap-class object of features ranked by test_statistic annotated with the proportion of subjects that showed active linkage of the features.
#' @export
-#'
-plot_differential_linkages <- function(differential_linkages, test_statistic, stat_range = c(0, 1),
- stat_ranking = c("ascending", "descending"), group_palette = NULL) {
+#'
+plot_differential_linkages <- function(
+ differential_linkages, test_statistic, stat_range = c(0, 1),
+ stat_ranking = c("ascending", "descending"), group_palette = NULL) {
if (!test_statistic %in% colnames(differential_linkages)) {
stop(paste0("test statistic '", test_statistic, "' not present in colnames(differential_linkages)"))
}
@@ -865,20 +897,27 @@ plot_differential_linkages <- function(differential_linkages, test_statistic, st
g_names_full <- colnames(df)[grepl("_n$", colnames(df)) & !grepl("^total_", colnames(df))]
g_names <- gsub("_n", "", g_names_full)
# proportion bar for linkage feature in all subjects
- ha_subject <- ComplexHeatmap::HeatmapAnnotation(subjects = ComplexHeatmap::anno_barplot(matrix(ncol = 2, c(df[["total_count"]],
- df[["total_n"]] - df[["total_count"]])), gp = grid::gpar(fill = c("black", "white"))), which = "row",
- annotation_name_gp = grid::gpar(fontsize = 8))
+ ha_subject <- ComplexHeatmap::HeatmapAnnotation(
+ subjects = ComplexHeatmap::anno_barplot(matrix(ncol = 2, c(
+ df[["total_count"]],
+ df[["total_n"]] - df[["total_count"]]
+ )), gp = grid::gpar(fill = c("black", "white"))), which = "row",
+ annotation_name_gp = grid::gpar(fontsize = 8)
+ )
ha_subject@anno_list$subjects@label <- "All\nSubjects"
# row annotation of linkage feature names
ha_name <- ComplexHeatmap::rowAnnotation(feat = ComplexHeatmap::anno_text(df[["feature"]], location = 0, rot = 0))
# plotted statistic for ordering results
mat <- matrix(df[[test_statistic]], ncol = 1)
rownames(mat) <- df[["feature"]]
- plot <- ComplexHeatmap::Heatmap(matrix = mat, cluster_rows = FALSE, left_annotation = ha_name, cell_fun = function(j,
- i, x, y, width, height, fill) {
+ plot <- ComplexHeatmap::Heatmap(matrix = mat, cluster_rows = FALSE, left_annotation = ha_name, cell_fun = function(
+ j,
+ i, x, y, width, height, fill) {
grid::grid.text(sprintf("%.3f", mat[i, j]), x, y, gp = grid::gpar(fontsize = 6))
- }, column_title = paste0(cluster, ": ", test_statistic), name = test_statistic, col = circlize::colorRamp2(breaks = stat_range,
- colors = stat_gradient), height = nrow(mat) * grid::unit(0.25, "in"), width = grid::unit(1, "in")) + ha_subject
+ }, column_title = paste0(cluster, ": ", test_statistic), name = test_statistic, col = circlize::colorRamp2(
+ breaks = stat_range,
+ colors = stat_gradient
+ ), height = nrow(mat) * grid::unit(0.25, "in"), width = grid::unit(1, "in")) + ha_subject
# generate an heatmap annotation for each category
if (is.null(group_palette)) {
group_palette <- ggplot_col_gen(length(g_names))
@@ -888,43 +927,48 @@ plot_differential_linkages <- function(differential_linkages, test_statistic, st
g <- g_names[i]
g_count <- paste0(g, "_count")
g_n <- paste0(g, "_n")
- ha <- ComplexHeatmap::HeatmapAnnotation(group = ComplexHeatmap::anno_barplot(matrix(ncol = 2, c(df[[g_count]], df[[g_n]] -
- df[[g_count]])), gp = grid::gpar(fill = c(group_palette[g], "#FFFFFF"))), name = g, which = "row",
- annotation_name_gp = grid::gpar(fontsize = 8))
+ ha <- ComplexHeatmap::HeatmapAnnotation(
+ group = ComplexHeatmap::anno_barplot(matrix(ncol = 2, c(df[[g_count]], df[[g_n]] -
+ df[[g_count]])), gp = grid::gpar(fill = c(group_palette[g], "#FFFFFF"))), name = g, which = "row",
+ annotation_name_gp = grid::gpar(fontsize = 8)
+ )
ha@anno_list$group@label <- g
plot <- plot + ha
}
return(plot)
}
+
#' Normalize a matrix to its max value by row or column
-#'
+#'
#' Normalizes a matrix to its max value by row or column
-#'
+#'
#' @param mat Matrix to be normalized
-#' @param dir Direction to normalize the matrix c('row', 'col')
+#' @param dir Direction to normalize the matrix c('row', 'col')
#' @return A normalized matrix in the direction specified.
#' @keywords internal
-#'
+#'
do_norm <- function(mat, dir) {
if (dir == "row") {
mat <- t(apply(mat, 1, function(x) {
- x/max(x)
+ x / max(x)
}))
return(mat)
} else if (dir == "col") {
mat <- apply(mat, 2, function(x) {
- x/max(x)
+ x / max(x)
})
return(mat)
}
}
+
#' Generate ggplot colors
-#'
+#'
#' Accepts a number of colors to generate and generates a ggplot color spectrum.
-#'
+#'
#' @param n Number of colors to generate
#' @return A vector of colors according to ggplot color generation.
#' @keywords internal
+#'
ggplot_col_gen <- function(n) {
hues <- seq(15, 375, length = n + 1)
return(grDevices::hcl(h = hues, l = 65, c = 100)[1:n])
diff --git a/R/processing_fxns.R b/R/processing_fxns.R
index 2b147d3..66d4e08 100644
--- a/R/processing_fxns.R
+++ b/R/processing_fxns.R
@@ -1,9 +1,9 @@
#' Calculate a signaling network for a domino object
-#'
+#'
#' This function calculates a signaling network. It requires a domino object
-#' preprocessed from create_domino and returns a domino object prepared for
+#' preprocessed from create_domino and returns a domino object prepared for
#' plotting with the various plotting functions in this package.
-#'
+#'
#' @param dom Domino object from create_domino.
#' @param max_tf_per_clust Maximum number of transcription factors called active in a cluster.
#' @param min_tf_pval Minimum p-value from differential feature score test to call a transcription factor active in a cluster.
@@ -11,16 +11,22 @@
#' @param rec_tf_cor_threshold Minimum pearson correlation used to consider a receptor linked with a transcription factor. Increasing this will decrease the number of receptors linked to each transcription factor.
#' @param min_rec_percentage Minimum percentage of cells in cluster expressing a receptor for the receptor to be linked to trancription factors in that cluster.
#' @return A domino object with a signaling network built
-#' @export
+#' @export
+#' @examples
+#' pbmc_dom_tiny_built <- build_domino(dom = domino2:::pbmc_dom_tiny, min_tf_pval = .001, max_tf_per_clust = 25,
+#' max_rec_per_tf = 25, rec_tf_cor_threshold = .25, min_rec_percentage = 0.1)
#'
-build_domino <- function(dom, max_tf_per_clust = 5, min_tf_pval = 0.01, max_rec_per_tf = 5, rec_tf_cor_threshold = 0.15,
- min_rec_percentage = 0.1) {
+build_domino <- function(
+ dom, max_tf_per_clust = 5, min_tf_pval = 0.01, max_rec_per_tf = 5, rec_tf_cor_threshold = 0.15,
+ min_rec_percentage = 0.1) {
if (dom@misc[["create"]] == FALSE) {
stop("Please run domino_create to create the domino object.")
}
dom@misc[["build"]] <- TRUE
- dom@misc[["build_vars"]] <- c(max_tf_per_clust = max_tf_per_clust, min_tf_pval = min_tf_pval,
- max_rec_per_tf = max_rec_per_tf, rec_tf_cor_threshold = rec_tf_cor_threshold, min_rec_percentage = min_rec_percentage)
+ dom@misc[["build_vars"]] <- c(
+ max_tf_per_clust = max_tf_per_clust, min_tf_pval = min_tf_pval,
+ max_rec_per_tf = max_rec_per_tf, rec_tf_cor_threshold = rec_tf_cor_threshold, min_rec_percentage = min_rec_percentage
+ )
if (length(dom@clusters)) {
# Get transcription factors for each cluster
clust_tf <- list()
@@ -31,8 +37,10 @@ build_domino <- function(dom, max_tf_per_clust = 5, min_tf_pval = 0.01, max_rec_
zeros <- names(ordered)[which(ordered == 0)]
fcs <- c()
for (zero in zeros) {
- fc <- mean(dom@features[zero, which(dom@clusters == clust)]) - mean(dom@features[zero,
- which(dom@clusters != clust)])
+ fc <- mean(dom@features[zero, which(dom@clusters == clust)]) - mean(dom@features[
+ zero,
+ which(dom@clusters != clust)
+ ])
fcs <- c(fcs, fc)
}
names(fcs) <- zeros
@@ -117,9 +125,9 @@ build_domino <- function(dom, max_tf_per_clust = 5, min_tf_pval = 0.01, max_rec_
# if complexes were used
inc_ligs_list <- lapply(inc_ligs, function(l) {
if (l %in% names(dom@linkages$complexes)) {
- return(dom@linkages$complexes[[l]])
+ return(dom@linkages$complexes[[l]])
} else {
- return(l)
+ return(l)
}
})
names(inc_ligs_list) <- inc_ligs
@@ -151,12 +159,12 @@ build_domino <- function(dom, max_tf_per_clust = 5, min_tf_pval = 0.01, max_rec_
# if complexes were used
cl_sig_list <- lapply(seq_along(inc_ligs_list), function(x) {
if (all(inc_ligs_list[[x]] %in% lig_genes)) {
- # Some of the ligands in the list object may not be present in the data
- if (length(inc_ligs_list[[x]]) > 1) {
- return(colMeans(cl_sig_mat[inc_ligs_list[[x]], ]))
- } else {
- return(cl_sig_mat[inc_ligs_list[[x]], ])
- }
+ # Some of the ligands in the list object may not be present in the data
+ if (length(inc_ligs_list[[x]]) > 1) {
+ return(colMeans(cl_sig_mat[inc_ligs_list[[x]], ]))
+ } else {
+ return(cl_sig_mat[inc_ligs_list[[x]], ])
+ }
}
})
names(cl_sig_list) <- names(inc_ligs_list)
@@ -192,15 +200,16 @@ build_domino <- function(dom, max_tf_per_clust = 5, min_tf_pval = 0.01, max_rec_
}
return(dom)
}
+
#' Pulls all items from a list pooled into a single vector
-#'
+#'
#' Helper function to convert from a nested series of lists to a single vector.
-#'
+#'
#' @param list List to pull items from
#' @param list_names Names of items in list to pool
#' @return A vector contaning all items in the list by list_names
#' @keywords internal
-#'
+#'
lc <- function(list, list_names) {
vec <- c()
for (name in list_names) {
diff --git a/R/utils.R b/R/utils.R
index cf6a1f9..efafb76 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,16 +1,19 @@
#' Access database
-#'
+#'
#' A function to pull database information from a domino object
-#'
+#'
#' @param dom A domino object that has been created
-#' @param name_only A boolean for whether to return only the name of the database used
+#' @param name_only A boolean for whether to return only the name of the database used
#' or the entire database that is stored. Default TRUE.
-#' @return A vector of unique databases used in building the domino object OR
+#' @return A vector of unique databases used in building the domino object OR
#' a data frame that includes the database information used in the domino object creation
#' @export
-
+#' @examples
+#' database_name <- dom_database(domino2:::pbmc_dom_built_tiny)
+#' full_database <- dom_database(domino2:::pbmc_dom_built_tiny, name_only = FALSE)
+#'
dom_database <- function(dom, name_only = TRUE) {
- db = slot(dom, "db_info")
+ db <- slot(dom, "db_info")
if (name_only) {
return(unique(db$database_name))
} else {
@@ -19,39 +22,47 @@ dom_database <- function(dom, name_only = TRUE) {
}
#' Access z-scores
-#'
+#'
#' A function to pull z-scored expression from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
#' @return A matrix containing the z-scored gene expression values for each gene (row) by cell (column)
#' @export
-
+#' @examples
+#' zscores <- dom_zscores(domino2:::pbmc_dom_built_tiny)
+#'
dom_zscores <- function(dom) {
slot(dom, "z_scores")
}
#' Access counts
-#'
+#'
#' A function to pull gene expression from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
#' @return A matrix containing the gene expression values for each gene (row) by cell (column)
#' @export
+#' @examples
+#' counts <- dom_counts(domino2:::pbmc_dom_built_tiny)
+#'
dom_counts <- function(dom) {
as.matrix(slot(dom, "counts"))
}
#' Access clusters
-#'
+#'
#' A function to pull cluster information from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
#' @param labels A boolean for whether to return the cluster labels for each cell or the clusters used for inferring communication
#' @return A vector containing either the names of the clusters used OR factors of the cluster label for each individual cell
#' @export
-
+#' @examples
+#' cluster_names <- dom_clusters(domino2:::pbmc_dom_built_tiny)
+#' cell_cluster_label <- dom_clusters(domino2:::pbmc_dom_built_tiny, labels = TRUE)
+#'
dom_clusters <- function(dom, labels = FALSE) {
- clust = slot(dom, "clusters")
+ clust <- slot(dom, "clusters")
if (labels) {
return(clust)
} else {
@@ -60,41 +71,53 @@ dom_clusters <- function(dom, labels = FALSE) {
}
#' Access transcription factor activation
-#'
+#'
#' A function to pull transcription factor activation scores from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
#' @return A matrix containing the transcription factor activation scores for each feature (row) by cell (column)
#' @export
+#' @examples
+#' tf_activation <- dom_tf_activation(domino2:::pbmc_dom_built_tiny)
+#'
dom_tf_activation <- function(dom) {
slot(dom, "features")
}
#' Access correlations
-#'
+#'
#' A function to pull receptor-transcription factor correlations from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
#' @return A matrix containing the correlation values for each receptor (row) by transcription factor (column)
#' @export
+#' @examples
+#' cor_matrix <- dom_correlations(domino2:::pbmc_dom_built_tiny)
+#'
dom_correlations <- function(dom) {
slot(dom, "cor")
}
#' Access linkages
-#'
+#'
#' A function to pull linkages from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
-#' @param link_type One value (out of "complexes", "receptor-ligand",
+#' @param link_type One value (out of "complexes", "receptor-ligand",
#' "tf-target", "tf-receptor", "receptor", "incoming-ligand") used
#' to select the desired type of linkage
#' @param by_cluster A boolean to indicate whether the linkages should be returned overall or by cluster
#' @return A list containing linkages between some combination of receptors, ligands, transcription factors, and clusters
#' @export
-dom_linkages <- function(dom, link_type = c("complexes", "receptor-ligand",
- "tf-target", "tf-receptor", "receptor", "incoming-ligand"), by_cluster = FALSE) {
- links = slot(dom, "linkages")
+#' @examples
+#' complexes <- dom_linkages(domino2:::pbmc_dom_built_tiny, "complexes")
+#' tf_rec_by_cluster <- dom_linkages(domino2:::pbmc_dom_built_tiny, "tf-receptor", TRUE)
+#'
+dom_linkages <- function(dom, link_type = c(
+ "complexes", "receptor-ligand",
+ "tf-target", "tf-receptor", "receptor", "incoming-ligand"
+ ), by_cluster = FALSE) {
+ links <- slot(dom, "linkages")
if (by_cluster) {
if (link_type == "tf-receptor") {
return(links$clust_tf)
@@ -110,9 +133,9 @@ dom_linkages <- function(dom, link_type = c("complexes", "receptor-ligand",
return(links$complexes)
} else if (link_type == "receptor-ligand") {
return(links$rec_lig)
- } else if (link_type == "tf_targets") {
+ } else if (link_type == "tf-target") {
return(links$tf_targets)
- } else if (link_type == "tf_receptor") {
+ } else if (link_type == "tf-receptor") {
return(links$tf_rec)
} else {
stop("This linkage type is not available.")
@@ -121,16 +144,19 @@ dom_linkages <- function(dom, link_type = c("complexes", "receptor-ligand",
}
#' Access signaling
-#'
+#'
#' A function to pull signaling matrices from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
#' @param cluster Either NULL to indicate global signaling or a specific cluster for which a
#' @return A data.frame containing the signaling score through each ligand (row) by each cluster (column) OR
#' a data.frame containing the global summed signaling scores between receptors (rows) and ligands (columns) of each cluster
#' @export
+#' @examples
+#' monocyte_signaling <- dom_signaling(domino2:::pbmc_dom_built_tiny, cluster = "CD14_monocyte")
+#'
dom_signaling <- function(dom, cluster = NULL) {
- if(is.null(cluster)) {
+ if (is.null(cluster)) {
as.data.frame(slot(dom, "signaling"))
} else {
as.data.frame(slot(dom, "cl_signaling_matrices")[cluster])
@@ -138,26 +164,86 @@ dom_signaling <- function(dom, cluster = NULL) {
}
#' Access differential expression
-#'
+#'
#' A function to pull differential expression p values from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
#' @return A matrix containing the p values for differential expression of transcription factors (rows) in each cluster (columns)
#' @export
+#' @examples
+#' de_mat <- dom_de(domino2:::pbmc_dom_built_tiny)
+#'
dom_de <- function(dom) {
slot(dom, "clust_de")
}
#' Access build information
-#'
+#'
#' A function to pull the parameters used when running [build_domino()] from a domino object
-#'
+#'
#' @param dom A domino object that has been created with [create_domino()]
#' @return A list containing booleans for whether the object has been created and build, and a list of the
#' build parameters that were used in [build_domino()] to infer the signaling network
#' @export
+#' @examples
+#' build_details <- dom_info(domino2:::pbmc_dom_built_tiny)
+#'
dom_info <- function(dom) {
- info = slot(dom, "misc")
- return(list("create" = info$create, "build"= info$build,
- "build_variables"= info$build_vars))
-}
\ No newline at end of file
+ info <- slot(dom, "misc")
+ return(list(
+ "create" = info$create, "build" = info$build,
+ "build_variables" = info$build_vars
+ ))
+}
+
+#' Access 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()] was 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
+#' @examples
+#' monocyte_receptors <- dom_network_items(domino2:::pbmc_dom_built_tiny, "CD14_monocyte", "receptors")
+#' all_tfs <- dom_network_items(domino2:::pbmc_dom_built_tiny, return = "features")
+#'
+dom_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]])
+ }
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 60ea17b..8b3e355 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -28,7 +28,6 @@ reference:
contents:
- domino-class
- linkage_summary
- - pbmc_dom
- title: "Analysis Functions"
desc: >
Functions for performing cell-cell communication analysis
@@ -68,14 +67,13 @@ reference:
- dom_de
- dom_info
- dom_linkages
+ - dom_network_items
- dom_signaling
- dom_tf_activation
- dom_zscores
- subtitle: "Misc. utility functions"
contents:
- add_rl_column
- - collate_network_items
- - convert_genes
- count_linkage
- mean_ligand_expression
- rename_clusters
diff --git a/docs/articles/domino2.html b/docs/articles/domino2.html
index 93a0495..eb48e90 100644
--- a/docs/articles/domino2.html
+++ b/docs/articles/domino2.html
@@ -98,67 +98,36 @@
Version 0.2
-
Using domino2 v0.2.1 for cell-cell communication inference
+
Using domino2 v0.2.2 for cell-cell communication inference
domino2 is a tool for analysis of intra- and intercellular signaling in single cell RNA sequencing (scRNAseq) data based on transcription factor (TF) activation. Here, we show a basic example pipeline for generating communication networks as a domino object. This vignette will demonstrate how to use domino2 on the 10X Genomics Peripheral Blood Mononuclear Cells (PBMC) data set of 2,700 cells. The data can be downloaded here.
Preprocessing scRNAseq data
-
Analysis of cell-cell communication with domino2 often follows initial processing and annotation of scRNAseq data. We are following the Satija Lab’s Guided Clustering Tutorial to prepare the data set for analysis with domino2.
+
Analysis of cell-cell communication with domino2 often follows initial processing and annotation of scRNAseq data. We used the Satija Lab’s Guided Clustering Tutorial to prepare the data set for analysis with domino2. The complete processing script is available in the data-raw directory of the domino2 package. The processed data can be downloaded from Zenodo.
-# Preprocessing of PBMC 3K tutorial data set using Seurat functions
-pbmc.data<-Read10X(data.dir ="../inst/extdata/pbmc3k_filtered_gene_bc_matrices")
-pbmc<-CreateSeuratObject(counts =pbmc.data, project ="pbmc3k", min.cells =3, min.features =200)
-pbmc[["percent.mt"]]<-PercentageFeatureSet(pbmc, pattern ="^MT-")
-pbmc<-NormalizeData(pbmc, normalization.method ="LogNormalize", scale.factor =10000)
-pbmc<-FindVariableFeatures(pbmc, selection.method ="vst", nfeatures =2000)
-pbmc<-ScaleData(pbmc, features =rownames(pbmc))
-pbmc<-RunPCA(pbmc, features =VariableFeatures(object =pbmc))
-pbmc<-FindNeighbors(pbmc, dims =1:10)
-pbmc<-FindClusters(pbmc, resolution =0.5)
-#> Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
-#>
-#> Number of nodes: 2700
-#> Number of edges: 97892
-#>
-#> Running Louvain algorithm...
-#> Maximum modularity in 10 random starts: 0.8719
-#> Number of communities: 9
-#> Elapsed time: 0 seconds
-pbmc<-RunUMAP(pbmc, dims =1:10)
-# Annotate clusters with cell phenotypes
-cell_dict<-data.frame(
- cluster =c(0:8),
- cell_type =c(
-"naive_CD4_T_cell", "CD14_monocyte", "memory_CD4_T_cell", "B_cell", "CD8_T_cell",
-"CD16_monocyte", "NK_cell", "dendritic_cell", "platelet"
-)
-)
-
-pbmc$cell_type<-
-plyr::mapvalues(
-pbmc$seurat_clusters,
- from =cell_dict$cluster,
- to =cell_dict$cell_type
-)
domino2 was designed to be compatible with Seurat objects in addition to accepting matrices and vectors of the data set as function inputs. Therefore, it’s a good idea to save the processed Seurat object for further analysis. However, pySCENIC is used for TF activation inference, and it requires as input a cell by gene matrix. This is the opposite orientation of Seurat objects but is the default for Python based tools such as scanpy. Saving our data as a loom file using loomR will be helpful for that step and is therefore recommended as well.
+
domino2 was designed to be compatible with Seurat objects in addition to accepting matrices and vectors of the data set as function inputs. Therefore, it’s a good idea to save the processed Seurat object for further analysis. However, pySCENIC is used for TF activation inference, and it requires as input a cell by gene matrix. This is the opposite orientation of Seurat objects but is the default for Python based tools such as scanpy. The RNA counts matrix can be saved as a tab-seperated value (.tsv) or comma-sperated value (.csv) file after transposing the matrix to cell by gene orientation.
tsv and csv files are very inefficient for storing large matrices. As an alternative, the counts matrix can be saved as a loom file and directly passed to pySCENIC in this format. Generating a loom file in R requires use of the loomR package. The package automatically handles the transposition of the counts matrix to cell by gene orientation as well. We recommend this approach to passing a counts matrix to pySCENIC. However, users should be aware that loomR is not hosted on CRAN or BioConductor at the time of this vignette’s creation.
+
+library(loomR)# save loom counts matrixpbmc_counts<-GetAssayData(object =pbmc, slot ="counts")
-pbmc_loom<-loomR::create(filename ="../inst/extdata/pbmc3k_counts.loom", data =pbmc_counts)
-pbmc_loom$close_all()# Remember to manually close connection to loom files!
+pbmc_loom<-loomR::create(filename =paste0(input_dir, "/pbmc3k_counts.loom"), data =pbmc_counts)
+# connection to the loom file must be closed to complete file creation
+pbmc_loom$close_all()
The cell x gene counts matrix may also be saved as a tab-separated value (.tsv) file if loomR cannot be used on your computing resource, but be aware that dense formatting of a .tsv can lead to very large file sizes and slow processing for large data sets. We have provided the necessary code below, but due to the large file size, this .tsv file is not provided in our data directory.
If you would like to skip ahead to the rest of the analysis, the processed objects from these steps are available in the data directory of this package.
@@ -168,7 +137,7 @@
Installationremotes package. The current stable version is v0.2.2, so we will install that branch.
A singularity image of SCENIC v0.12.1 can be installed from the DockerHub image as a source (example scripts for use with Docker are provided in the scenic_bash directory as well). SCENIC requires a list of TFs, motif annotations, and cisTarget motifs which are all available from the authors of SCENIC for human (HGNC), mouse (MGI), and fly. The following will download everything necessary for an analysis of a data set with HGNC gene labels for the hg38 genome.
-
SCENIC_DIR="../inst/extdata/scenic"
+
SCENIC_DIR="'temp_dir'/scenic"mkdir${SCENIC_DIR}# Build singularity imagesingularity build "${SCENIC_DIR}/aertslab-pyscenic-0.12.1.sif" docker://aertslab/pyscenic:0.12.1
-# Matrix containing motifs as rows and genes as columns and ranking position for each gene and motif (based on CRM scores) as values
-# below is a feather v1 format link; feather v2 required for this version of pySCENIC
-# "https://resources.aertslab.org/cistarget/databases/old/homo_sapiens/hg38/refseq_r80/mc9nr/gene_based/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather"
-# Feather v2 link:
-curl"https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather" \
- -o "${SCENIC_DIR}/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather"
-
-# below is a feather v1 format link; feather v2 required for this version of pySCENIC
-# "https://resources.aertslab.org/cistarget/databases/old/homo_sapiens/hg38/refseq_r80/mc9nr/gene_based/hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather"
-# Feather v2 link:
-curl"https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather" \
- -o "${SCENIC_DIR}/hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather"
-
-# List of genes encoding TFs in the HG38 reference genome
-curl"https://resources.aertslab.org/cistarget/tf_lists/allTFs_hg38.txt" \
- -o "${SCENIC_DIR}/allTFs_hg38.txt"
-
-# Motif annotations based on the 2017 cisTarget motif collection. Use these files if you are using the mc9nr databases.
-curl"https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl" \
- -o "${SCENIC_DIR}/motifs-v9-nr.hgnc-m0.001-o0.0.tbl"
+# Matrix containing motifs as rows and genes as columns and ranking position for each gene and motif (based on CRM scores) as values. URLs provided link to v2 feather files required for the 0.12.1 version of pySENIC.
+
+curl"https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather" \
+ -o "${SCENIC_DIR}/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather"
+
+curl"https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather" \
+ -o "${SCENIC_DIR}/hg38_500bp_up_100bp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather"
+
+# List of genes encoding TFs in the HG38 reference genome
+curl"https://resources.aertslab.org/cistarget/tf_lists/allTFs_hg38.txt" \
+ -o "${SCENIC_DIR}/allTFs_hg38.txt"
+
+# Motif annotations based on the 2017 cisTarget motif collection. Use these files if you are using the mc9nr databases.
+curl"https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl" \
+ -o "${SCENIC_DIR}/motifs-v9-nr.hgnc-m0.001-o0.0.tbl"
Downloads for CellPhoneDB:
Database files from CellPhoneDB v4.0.0 for human scRNAseq data can be installed from a public Github repository from the Tiechmann Group that developed CellPhoneDB.
-
CELLPHONE_DIR="../inst/extdata/cellphoneDB_v4"
-mkdir${CELLPHONE_DIR}
-
-# Annotation of protein complexes
-curl"https://raw.githubusercontent.com/Teichlab/cellphonedb-data/master/data/complex_input.csv" \
- -o "${CELLPHONE_DIR}/complex_input.csv"
-# Annotation of genes encoding database proteins
-curl"https://raw.githubusercontent.com/Teichlab/cellphonedb-data/master/data/gene_input.csv" \
- -o "${CELLPHONE_DIR}/gene_input.csv"
-# Annotated protein interactions
-curl"https://raw.githubusercontent.com/Teichlab/cellphonedb-data/master/data/interaction_input.csv" \
- -o "${CELLPHONE_DIR}/interaction_input.csv"
-# Annotated protein features
-curl"https://raw.githubusercontent.com/Teichlab/cellphonedb-data/master/data/protein_input.csv" \
- -o "${CELLPHONE_DIR}/protein_input.csv"
+
+#
+cellphone_url<-"https://github.com/ventolab/cellphonedb-data/archive/refs/tags/v4.0.0.tar.gz"
+
+# download compressed database
+cellphone_tar<-paste0(temp_dir, "/cellphoneDB_v4.tar.gz")
+download.file(url =cellphone_url, destfile =cellphone_tar)
+
+# move contents of the compressed file to a new directory
+untar(tarfile =cellphone_tar, exdir =cellphone_dir)
+cellphone_data<-paste0(cellphone_dir, "/cellphonedb-data-4.0.0/data")
@@ -299,24 +259,19 @@
Loading SCENIC Results
The TF activities are a required input for domino2, and the regulons learned by SCENIC are an important but optional input needed to annotate TF-target interactions and to prune TF-receptor linkages where the receptor is a target of the TF. This prevents the distinction of receptor expression driving TF activity or the TF inducing the receptor’s expression.
The initial regulons data frame read into R from the ctx function has 2 rows of column names that need to be replaced with one succinct description. domino2 has changed the input format for TF regulons to be a list storing vectors of target genes in each regulon, where the names of the list are the TF genes. This facilitates the use of alternative methods for TF activity quantification methods. We provide a helper function, create_regulon_list_scenic() for easy retrieval of TF regulons from the output of the pySCENIC ctx function.
The initial regulons data frame read into R from the ctx function has 2 rows of column names that need to be replaced with one succinct description. domino2 has changed the input format for TF regulons to be a list storing vectors of target genes in each regulon, where the names of the list are the TF genes. This facilitates the use of alternative methods for TF activity quantification. We provide a helper function, create_regulon_list_scenic() for easy retrieval of TF regulons from the output of the pySCENIC ctx function.
Users should be aware that the AUC matrix from SCENIC is loaded in a cell x TF orientation and should be transposed to TF x cell orientation. pySCENIC also appends “(+)” to TF names that are converted to “…” upon loading into R. These characters can be included without affecting the results of domino2 analysis, but can be confusing when querying TF features in the data. We recommend comprehensive removal of the “…” characters using the gsub() function.
auc_in<-as.data.frame(t(auc))
-# Remove pattern "..." from the end of all rownames:
+# Remove pattern '...' from the end of all rownames:rownames(auc_in)<-gsub("\\.\\.\\.$", "", rownames(auc_in))
The change to use rl_map formatting also enables users to manually append interactions of interest that are not included in the interaction database if need be. This can be attained by formatting the desired interactions as a data.frame with the same column headers as the rl_map and using the rbind() function.
# Integrin complexes are not annotated as receptors in CellPhoneDB_v4.0
-# collagen-integrin interactions between cells may be missed unless tables from the CellPhoneDB reference are edited or the interactions are manually added
+# collagen-integrin interactions between cells may be missed unless tables from
+# the CellPhoneDB reference are edited or the interactions are manually added
-col_int_df<-data.frame(
-"int_pair"="a11b1 complex & COLA1_HUMAN",
-"name_A"="a11b1 complex", "uniprot_A"="P05556,Q9UKX5", "gene_A"="ITB1,ITA11", "type_A"="R",
-"name_B"="COLA1_HUMAN", "uniprot_B"="P02452,P08123", "gene_B"="COL1A1,COL1A2", "type_B"="L",
-"annotation_strategy"="manual", "source"="manual", "database_name"="manual"
-)
+col_int_df<-data.frame(int_pair ="a11b1 complex & COLA1_HUMAN", name_A ="a11b1 complex",
+ uniprot_A ="P05556,Q9UKX5", gene_A ="ITB1,ITA11", type_A ="R", name_B ="COLA1_HUMAN",
+ uniprot_B ="P02452,P08123", gene_B ="COL1A1,COL1A2", type_B ="L", annotation_strategy ="manual",
+ source ="manual", database_name ="manual")rl_map_append<-rbind(col_int_df, rl_map)knitr::kable(head(rl_map_append))
-
-
-
+
+
+
+
+
-
-
-
-
+
+
@@ -528,78 +489,78 @@
Optional: Adding interactions man
manual
-
111
-
FGF4 & FGFR1
-
FGF4
-
P08620
-
FGF4
-
L
-
FGFR1
-
P11362
-
FGFR1
+
4
+
RAreceptor_RXRG & atRetinoicAcid_byALDH1A3
+
RAreceptor_RXRG
+
P48443,P29373
+
RXRG,CRABP2
R
-
I2D
-
+
atRetinoicAcid_byALDH1A3
+
P47895
+
ALDH1A3
+
L
+
curated
+
uniprot;reactome
CellPhoneDB_v4.0
-
2
-
ADORA3 & ENTPD1
-
ADORA3
-
P0DMS8
-
ADORA3
+
5
+
RAreceptor_RXRG & atRetinoicAcid_byALDH1A2
+
RAreceptor_RXRG
+
P48443,P29373
+
RXRG,CRABP2
R
-
ENTPD1
-
P49961
-
ENTPD1
+
atRetinoicAcid_byALDH1A2
+
O94788
+
ALDH1A2
L
curated
-
uniprot
+
uniprot;reactome
CellPhoneDB_v4.0
-
3
-
EPHB3 & EFNB1
-
EPHB3
-
P54753
-
EPHB3
+
6
+
RAreceptor_RXRG & atRetinoicAcid_byALDH1A1
+
RAreceptor_RXRG
+
P48443,P29373
+
RXRG,CRABP2
R
-
EFNB1
-
P98172
-
EFNB1
+
atRetinoicAcid_byALDH1A1
+
P00352
+
ALDH1A1
L
curated
-
PMID: 15114347
+
uniprot;reactome
CellPhoneDB_v4.0
-
4
-
TNFSF18 & TNFRSF18
-
TNFSF18
-
Q9UNG2
-
TNFSF18
-
L
-
TNFRSF18
-
Q9Y5U5
-
TNFRSF18
+
7
+
RAreceptor_RXRB & atRetinoicAcid_byALDH1A3
+
RAreceptor_RXRB
+
P28702,P29373
+
RXRB,CRABP2
R
+
atRetinoicAcid_byALDH1A3
+
P47895
+
ALDH1A3
+
L
curated
-
uniprot
+
uniprot;reactome
CellPhoneDB_v4.0
-
5
-
EFNA1 & EPHA1
-
EFNA1
-
P20827
-
EFNA1
-
L
-
EPHA1
-
P21709
-
EPHA1
+
8
+
RAreceptor_RXRB & atRetinoicAcid_byALDH1A2
+
RAreceptor_RXRB
+
P28702,P29373
+
RXRB,CRABP2
R
+
atRetinoicAcid_byALDH1A2
+
O94788
+
ALDH1A2
+
L
curated
-
PMID: 15114347
+
uniprot;reactome
CellPhoneDB_v4.0
@@ -639,37 +600,20 @@
Create Domino objectcreate_domino() function can be used to make the object. Parameters of note include “use_clusters” which is required to assess signaling between cell types rather than linkage of TFs and receptors broadly across the data set. “use_complexes” decides if receptors that function in heteromeric complexes will be considered in testing linkages between TFs and receptors. If TRUE, a receptor complex is only linked to a TF if a majority of the component genes meet the Spearman correlation threshold. “remove_rec_dropout” decides whether receptor values of 0 will be considered in correlation calculations; this measure is intended to reduce the effect of dropout on receptors with low expression.
-pbmc_dom<-create_domino(
- rl_map =rl_map,
- features =auc_in,
- counts =counts,
- z_scores =z_scores,
- clusters =clusters,
- tf_targets =regulon_list,
- use_clusters =TRUE,
- use_complexes =TRUE,
- remove_rec_dropout =FALSE
-)
-# rl_map: receptor-ligand map data frame
-# features: TF scores (AUC matrix)
-# counts: counts matrix
-# z_scores: scaled expression data
-# clusters: named vector of cell cluster assignments
-# tf_targets: list of TFs and their regulons
-# use_clusters: assess receptor activation and ligand expression on a per-cluster basis
-# use_complexes: include receptors and genes that function as a complex in results
-# remove_rec_dropout: whether to remove zeroes from correlation calculations
+pbmc_dom<-create_domino(rl_map =rl_map, features =auc_in, counts =counts, z_scores =z_scores,
+ clusters =clusters, tf_targets =regulon_list, use_clusters =TRUE, use_complexes =TRUE,
+ remove_rec_dropout =FALSE)
+# rl_map: receptor-ligand map data frame features: TF scores (AUC matrix)
+# counts: counts matrix z_scores: scaled expression data clusters: named vector
+# of cell cluster assignments tf_targets: list of TFs and their regulons
+# use_clusters: assess receptor activation and ligand expression on a
+# per-cluster basis use_complexes: include receptors and genes that function as
+# a complex in results remove_rec_dropout: whether to remove zeroes from
+# correlation calculations
The above results will be equivalent to using the Seurat object directly. Note that inclusion of both the matrix inputs and the Seurat object will prioritize the Seurat object.
Build Dominobuild_domino() finalizes the construction of the domino object by setting parameters for identifying TFs with differential activation between clusters, receptor linkage with TFs based on magnitude of positive correlation, and the minimum percentage of cells within a cluster that have expression of a receptor for the receptor to be called as active.
There are also options for thresholds of the number of TFs that may be called active in a cluster and the number of receptors that may be linked to any one TF. For thresholds of n TFs and m receptors, the bottom n TFs by lowest p-values from the wilcoxon rank sum test and the top m receptors by Spearman correlation coefficient are chosen.
-pbmc_dom<-build_domino(
- dom =pbmc_dom,
- min_tf_pval =.001,
- max_tf_per_clust =25,
- max_rec_per_tf =25,
- rec_tf_cor_threshold =.25,
- min_rec_percentage =0.1
-)
-# min_tf_pval: Threshold for p-value of DE for TFs
-# rec_tf_cor_threshold: Minimum correlation between receptor and TF
-# min_rec_percentage: Minimum percent of cells that must express receptor
+pbmc_dom<-build_domino(dom =pbmc_dom, min_tf_pval =0.001, max_tf_per_clust =25,
+ max_rec_per_tf =25, rec_tf_cor_threshold =0.25, min_rec_percentage =0.1)
+# min_tf_pval: Threshold for p-value of DE for TFs rec_tf_cor_threshold:
+# Minimum correlation between receptor and TF min_rec_percentage: Minimum
+# percent of cells that must express receptor
Both of the thresholds for the number of receptors and TFs can be sent to infinity (Inf) to collect all receptors and TFs that meet statistical significance thresholds.
Cumulative signaling between cell types
@@ -719,7 +651,7 @@
Cumulative signaling between ce
The cumulative degree of signaling between clusters is assessed as the sum of the scaled expression of ligands targeting active receptors on another cluster. This can be visualized in graph format using the signaling_network() function. Nodes represent each cell cluster and edges scale with the magnitude of signaling between the clusters. The color of the edge corresponds to the sender cluster for that signal.
Signaling networks can also be drawn with the edges only rendering the signals directed towards a given cell type or signals from one cell type directed to others. To see those options in use, as well as other plotting functions and their options, please see our plotting vignette.
@@ -728,22 +660,21 @@
Specific Signaling Int
Beyond the aggregated degree of signaling between cell types, the degrees of signaling through specific ligand-receptor interactions can be assessed. gene_network() provides a graph of linkages between active TFs in a cluster, the linked receptors in that cluster, and the possible ligands of these active receptors.
New to domino2, gene_network() can be used between two clusters to determine if any of the possible ligands for a given receptor are expressed by a putative outgoing signaling cluster.
Another form of comprehensive ligand expression assessment is available for individual active receptors in the form of circos plots new to domino2. The outer arcs correspond to clusters in the domino object with inner arcs representing each possible ligand of the plotted receptor. Arcs are drawn between ligands on a cell type and the receptor if the ligand is expressed above the specified threshold. Arc widths correspond to the mean express of the ligand by the cluster with the widest arc width scaling to the maximum expression of the ligand within the data.
In order to better implement some of the changes that have been made in domino2, the structure of the domino object has been modified. Here, we describe the data that is stored in the domino object, as well as how best to access it for any other use you may have in mind.
+
In order to better implement some of the changes that have been made in domino2, the structure of the domino object has been modified. Here, we describe the data that is stored in the domino object, as well as how best to access this data.
Object contents
-
The contents of the domino object class include a few broad groups: - input data, such as the L-R database used, counts, and z-scores - calculations, such as correlations or differential expression results - linkages, which show connections between TFs, receptors, and ligands across the data set and within clusters - signaling matrices, which show signaling from ligands to receptors in clusters (or a summary of overall signaling) - build information, which includes the parameters used to build the object in the build_domino() functions
+
The contents of the domino object class include a few broad groups: - input data, such as the L-R database used, counts, and z-scores - calculations, such as correlations or differential expression results - linkages, which show connections between TFs, receptors, and ligands across the data set and within clusters - signaling matrices, which show expression of ligands signaling to active receptors in clusters (or a summary of overall signaling) - build information, which includes the parameters used to build the object in the build_domino() functions
For commonly accessed information (the number of cells, clusters, and some build information), the show and print methods for domino objects can be used.
-
+
dom#> A domino object of 2700 cells#> Built with signaling between 9 clusters
-
+
print(dom)
-#> A domino object of 2700 cells
-#> Contains signaling between 9 clusters
-#> Built with a maximum of 25 TFs per cluster
-#> and a maximum of 25 receptors per TF
+#> A domino object of 2700 cells
+#> Built with signaling between 9 clusters
Access functions
In order to facilitate access to the information stored in the domino object, we have provided a collection of functions to retrieve specific items. These functions begin with “dom_” and can be listed using ls().
When creating a domino object with the create_domino() function, several inputs are required which are then stored in the domino object itself. These include cluster labels, the counts matrix, z-scored counts matrix, transcription factor activation scores, and the R-L database used in create_rl_map_cellphonedb().
For example, to access the cluster names in the domino object:
Input datadom_database() function. By default, the function returns the name(s) of the database(s) used:
-
+
Information about the database referenced for ligand-receptor pairs and composition of protien complexes can be extracted from the dom_database() function. By default, the function returns the name(s) of the database(s) used:
Active transcription factors in each cluster are determined by conducting wilcoxon rank sum tests for each transcription factor where the trascription factor activity scores amongst all cells in the cluster are tested against the activity scores of all cells outside of the cluster. The p-value for the one-sided test for greater activity within the cluster compared to outside can be accessed with the dom_de() function.
-
+
Active transcription factors in each cluster are determined by conducting wilcoxon rank sum tests for each transcription factor where the transcription factor activity scores amongst all cells in a cluster are tested against the activity scores of all cells outside of the cluster. The p-values for the one-sided test for greater activity within the cluster compared to outside can be accessed with the dom_de() function.
Linkages between various ligands, receptors, and transcription factors can be accessed in several different ways, depending on the specific link and the scope desired. The dom_linkages() function has three arguments - the first, like all of our access functions, is for the domino object. The second, link_type, is used to specify which linkages are desired (options are complexes, receptor-ligand, tf-target, or tf-receptor). The third argument, by_cluster, determines whether the linkages returned are arranged by cluster (though this does change the available linkage types to tf-receptor, receptor, or incoming-ligand). For example, to access the complexes used across the dataset:
-
+
Linkages between ligands, receptors, and transcription factors can be accessed in several different ways, depending on the specific link and the scope desired. The dom_linkages() function has three arguments - the first, like all of our access functions, is for the domino object. The second, link_type, is used to specify which linkages are desired (options are complexes, receptor-ligand, tf-target, or tf-receptor). The third argument, by_cluster, determines whether the linkages returned are arranged by cluster (though this does change the available linkage types to tf-receptor, receptor, or incoming-ligand). For example, to access the complexes used across the dataset:
+
complex_links<-dom_linkages(dom, link_type ="complexes")# Look for components of NODAL receptor complexcomplex_links$NODAL_receptor#> [1] "ACVR1B" "ACVR2B" "TDGF1" "character(0)"
If, for some reason, you find yourself in need of the entire linkage structure (not recommended), it can be accessed through its slot name; domino objects are S4 objects.
-
+
all_linkages<-slot(dom, "linkages")# Names of all sub-structures:names(all_linkages)#> [1] "complexes" "rec_lig" "tf_targets" #> [4] "clust_tf" "tf_rec" "clust_tf_rec" #> [7] "clust_rec" "clust_incoming_lig"
-
Alternately, to obtain a simplified list of receptors, ligands, and/or features in the domino object, use the collate_network_items() function. To pull all receptors associated with the dendritic cell cluster:
Alternately, to obtain a simplified list of receptors, ligands, and/or features in the domino object, use the dom_network_items() function. To pull all transcription factors associated with the dendritic cell cluster:
To keep track of the options set when running build_domino(), they are stored within the domino object itself. To view these options, use the dom_info() function.
In addition to the numerous plotting options available in the original Domino, domino2 has added more functionality and new methods to improve visualization and interpretation of analysis results. Here, we will go over the different plotting functions available, as well as different options that can be utilized, and options for customization. We will be using the domino object (loaded as dom) from our example analysis in Getting Started.
In addition to displaying the scores for the correlations, this function can also be used to identify correlations above a certain value (using the bool and bool_thresh arguments) or to identify the combinations of receptors and transcription factors (TFs) that are connected (with argument mark_connections).
The heatmap functions in domino2 are based on ComplexHeatmap::Heatmap() and will also accept additional arguments meant for that function. For example, while an argument for clustering rows or columns of the heatmap is not explicitly stated, they can still be passed to ComplexHeatmap::Heatmap() through cor_heatmap().
It functions similarly to cor_heatmap(), with arguments to select a specific vector of features, to use a boolean view with a threshold, and to pass other arguments to NMF::aheatmap(). Specific to this function though are arguments for setting the range of values to be visualized and one to choose whether or not to normalize the scores to the max value.
incoming_signaling_heatmap() can be used to visualize the cluster average expression of the ligands capable of activating the TFs enriched in the cluster. For example, to view the incoming signaling of the CD8 T cells:
We can also select for specific clusters of interest that are signaling to the CD8 T cells. If we are only interested in viewing the monocyte signaling:
As with our other heatmap functions, options are available for a minimum threshold, maximum threshold, whether to scale the values after thresholding, whether to normalize the matrix, and the ability to pass further arguments to NMF::aheatmap().
@@ -164,11 +169,11 @@
Heatmap of Incoming Signali
Heatmap of Signaling Between Clusters
signaling_heatmap() makes a heatmap showing the signaling strength of ligands from each cluster to receptors of each cluster based on averaged expression.
gene_network() makes a network plot to display signaling in selected clusters such that ligands, receptors and features associated with those clusters are displayed as nodes with edges as linkages. To look at signaling to the CD16 Monocytes from the CD14 Monocytes:
In addition to changes such as scaling, thresholds, or layouts, this plot can be modified by selecting incoming or outgoing clusters (or both!). An example to view signaling from the CD14 Monocytes to other cluster:
A new addition to domino2, circos_ligand_receptor() creates a chord plot showing ligands that can activate a specific receptor, displaying mean cluster expression of the ligand with the width of the chord.
This function can be given cluster colors to match other plots you may make with your data. In addition, the plot can be adjusted by changing the threshold of ligand expression required for a linkage to be visualized or selecting clusters of interest.
cor_scatter() can be used to plot each cell based on expression of the selected TF and receptor. This produces a scatter plot as well as a line of best fit to look at receptor-TF correlation.
Do keep in mind that there is an argument for remove_rec_dropout that should match the parameter that was used when the domino object was built. In this case, we did not use that build parameter, so we will leave the value in this argument as its default value of FALSE.
diff --git a/docs/reference/circos_ligand_receptor-1.png b/docs/reference/circos_ligand_receptor-1.png
new file mode 100644
index 0000000..27b513e
Binary files /dev/null and b/docs/reference/circos_ligand_receptor-1.png differ
diff --git a/docs/reference/circos_ligand_receptor.html b/docs/reference/circos_ligand_receptor.html
index 0e3462b..fadf9a7 100644
--- a/docs/reference/circos_ligand_receptor.html
+++ b/docs/reference/circos_ligand_receptor.html
@@ -108,6 +108,22 @@
Value
renders a circos plot to the active graphics device
+
+
Examples
+
cols<-c(
+"red", "orange", "green", "blue", "pink", "purple",
+"slategrey", "firebrick", "hotpink"
+)
+names(cols)<-dom_clusters(domino2:::pbmc_dom_built_tiny, labels =FALSE)
+circos_ligand_receptor(domino2:::pbmc_dom_built_tiny, receptor ="FCER2", cell_colors =cols)
+#> There are more than one numeric columns in the data frame. Take the
+#> first two numeric columns and draw the link ends with unequal width.
+#>
+#> Type `circos.par$message = FALSE` to suppress the message.
+
+
+
+
diff --git a/docs/reference/cor_heatmap-1.png b/docs/reference/cor_heatmap-1.png
new file mode 100644
index 0000000..b547e3d
Binary files /dev/null and b/docs/reference/cor_heatmap-1.png differ
diff --git a/docs/reference/cor_heatmap.html b/docs/reference/cor_heatmap.html
index bd9b58b..1b65e63 100644
--- a/docs/reference/cor_heatmap.html
+++ b/docs/reference/cor_heatmap.html
@@ -123,6 +123,15 @@
pbmc_dom_tiny_all<-pbmc_dom_tiny<-create_domino(rl_map =domino2:::rl_map_tiny,
+ features =domino2:::auc_tiny, counts =domino2:::RNA_count_tiny, z_scores =domino2:::RNA_zscore_tiny,
+ clusters =domino2:::clusters_tiny, tf_targets =domino2:::regulon_list_tiny, use_clusters =FALSE,
+ use_complexes =FALSE, rec_min_thresh =0.1, remove_rec_dropout =TRUE,
+ tf_selection_method ="all")
+#> Reading in and processing signaling database
+#> Database provided from source: CellPhoneDB
+#> Getting z_scores, clusters, and counts
+#> Calculating correlations
+#> 1 of 3
+#> 2 of 3
+#> 3 of 3
+#> FAS has component genes that did not pass testing parameters
+#> CD22 has component genes that did not pass testing parameters
+#> IGF1R has component genes that did not pass testing parameters
+
+pbmc_dom_tiny_clustered<-create_domino(rl_map =domino2:::rl_map_tiny,
+ features =domino2:::auc_tiny, counts =domino2:::RNA_count_tiny, z_scores =domino2:::RNA_zscore_tiny,
+ clusters =domino2:::clusters_tiny, tf_targets =domino2:::regulon_list_tiny,
+ use_clusters =TRUE, use_complexes =TRUE, remove_rec_dropout =FALSE)
+#> Reading in and processing signaling database
+#> Database provided from source: CellPhoneDB
+#> Getting z_scores, clusters, and counts
+#> Calculating feature enrichment by cluster
+#> 1 of 3
+#> 2 of 3
+#> 3 of 3
+#> Calculating correlations
+#> 1 of 3
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#> 2 of 3
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#> 3 of 3
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#>Warning: Cannot compute exact p-value with ties
+#> IGF1R has component genes that did not pass testing parameters
+
+
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.
diff --git a/docs/reference/feat_heatmap-1.png b/docs/reference/feat_heatmap-1.png
new file mode 100644
index 0000000..4dcd592
Binary files /dev/null and b/docs/reference/feat_heatmap-1.png differ
diff --git a/docs/reference/feat_heatmap.html b/docs/reference/feat_heatmap.html
index 8c3a971..515bff5 100644
--- a/docs/reference/feat_heatmap.html
+++ b/docs/reference/feat_heatmap.html
@@ -138,6 +138,14 @@
Value
a Heatmap rendered to the active graphics device
+
+
Examples
+
feat_heatmap(domino2:::pbmc_dom_built_tiny, min_thresh =0.1, max_thresh =0.6, norm =TRUE)
+#>Warning: You are using norm with min_thresh and max_thresh. Note that values will be thresholded AFTER normalization.
+
+
+
+
diff --git a/docs/reference/gene_network-1.png b/docs/reference/gene_network-1.png
new file mode 100644
index 0000000..c66720e
Binary files /dev/null and b/docs/reference/gene_network-1.png differ
diff --git a/docs/reference/gene_network.html b/docs/reference/gene_network.html
index 784d406..9b64962 100644
--- a/docs/reference/gene_network.html
+++ b/docs/reference/gene_network.html
@@ -129,6 +129,13 @@
Value
an igraph rendered to the active graphics device
+
a printed description of the number of cells in a domino object and its build status
+
+
Examples
+
domino2:::pbmc_dom_built_tiny
+#> A domino object of 360 cells
+#> Built with signaling between 3 clusters
+
+show(domino2:::pbmc_dom_built_tiny)
+#> A domino object of 360 cells
+#> Built with signaling between 3 clusters
+
+
+
diff --git a/docs/reference/signaling_heatmap-1.png b/docs/reference/signaling_heatmap-1.png
new file mode 100644
index 0000000..b7d7fbb
Binary files /dev/null and b/docs/reference/signaling_heatmap-1.png differ
diff --git a/docs/reference/signaling_heatmap.html b/docs/reference/signaling_heatmap.html
index 2d50581..4298e66 100644
--- a/docs/reference/signaling_heatmap.html
+++ b/docs/reference/signaling_heatmap.html
@@ -118,6 +118,13 @@
Value
a Heatmap rendered to the active graphics device
+