In this tutorial, we will walk through the workflow for performing differential signaling analysis using the dominoSignal package. While the other tutorials focus on signaling networks inferred from a single data set, many studies involve multiple subjects or conditions, requiring the comparison of signaling across these groups. The dominoSignal package implements the Differential Cell Signaling Test (DCST), a framework for identifying receptor, ligand, and transcription factor (TF) linkages that are significantly differentially present across groups of subjects. This method is described in detail in Differential cell signaling testing for cell-cell communication inference from single-cell data by dominoSignal (2026) published in Bioinformatics. The publication also includes simulation benchmarking and application of the method to pancreatic cancer and immunotherapy data sets. Refer to the publication for full details on the method and biological validation of the workflow.
The differential signaling workflow consists of the following steps:
domino() object per subject to be compared (see
the Getting
Started tutorial for details on building a domino()
object).summarize_linkages() function.linkage_summary object to include
only the subjects of interest.count_linkage() function.test_differential_linkages() function, which implements a
Fisher’s exact test.plot_differential_linkages() function.The DCST workflow can be performed with 2 or more groups of subjects. In this tutorial, we will demonstrate the workflow for both a two-group comparison (condition A vs condition B) and a comparison across more than two groups.
In this tutorial, we use the domino() object built in
the Getting
Started tutorial. If you have not yet built a domino()
object, you can do so by following the instructions on that page or use
the code below to load the pre-built domino() object.
set.seed(42)
library(dominoSignal)
library(dplyr)
library(kableExtra)
# BiocFileCache helps with file management across sessions
bfc <- BiocFileCache::BiocFileCache(ask = FALSE)
data_url <- "https://zenodo.org/records/10951634/files/pbmc_domino_built.rds"
tmp_path <- BiocFileCache::bfcrpath(bfc, data_url)
dom <- readRDS(tmp_path)Because the differential signaling workflow compares linkages across
multiple domino() objects (one per subject), we will create
six synthetic domino() objects from the same PBMC dataset.
For each synthetic subject, we randomly retain a subset of the cells
without replacement while keeping their true cluster labels and then
perturb the cluster labels of a subset of those cells by reshuffling
within that subset. This approach preserves some amount of per-cluster
signaling while making sure each domino() object is
distinct from the others, giving the example workflow something to
compare.1
Note that reshuffling is only for demonstration purposes and does not
represent any biological differences. In practice, you would build one
domino() object from each biological subject or sample to
be used in the comparison. To ensure results are comparable across
subjects, use the same parameters in create_domino() and
build_domino() for all subjects. Refer to the DCST
publication for an example of this workflow applied to real
data.
To start, we extract the inputs from the original
domino() object so we can supply them to the
create_domino() function for the synthetic subjects.
rl_map <- dom_database(dom, name_only = FALSE)
features <- dom_tf_activation(dom)
# Access counts slot directly to maintain sparse matrix format
counts <- slot(dom, "counts")
zscores <- dom_zscores(dom)
tf_targets <- dom_linkages(dom, link_type = "tf-target")
clusters <- dom_clusters(dom, labels = TRUE)We can use a helper function to resample the cells and shuffle the
cluster labels for a subset of them, and then rebuild the
domino() object with create_domino() and
build_domino() using the same parameters across all
subjects so that they can be compared.2
build_reps <- function(clusters, seed, keep_frac, shuffle_frac) {
# Local seed for reproducibility without affecting global seed
withr::with_seed(seed, {
n_cells <- length(clusters)
keep_cells <- sort(sample.int(n_cells, size = round(n_cells * keep_frac),
replace = FALSE))
resampled_clusters <- clusters[keep_cells]
features <- features[, keep_cells, drop = FALSE]
counts <- counts[, keep_cells, drop = FALSE]
zscores <- zscores[, keep_cells, drop = FALSE]
shuffled_idx <- sample(seq_along(resampled_clusters), size = round(shuffle_frac *
length(resampled_clusters)), replace = FALSE)
resampled_clusters[shuffled_idx] <- sample(resampled_clusters[shuffled_idx])
dom_shuffle <- create_domino(rl_map = rl_map, features = features, counts = counts,
z_scores = zscores, clusters = resampled_clusters, tf_targets = tf_targets,
use_complexes = TRUE, remove_rec_dropout = FALSE, verbose = FALSE)
dom_shuffle <- build_domino(dom = dom_shuffle, max_tf_pval = 0.05, max_tf_per_clust = 25,
max_rec_per_tf = 25, rec_tf_cor_threshold = 0.1, min_rec_percentage = 0.05)
})
return(dom_shuffle)
}
params <- data.frame(set = c(161803, 314159, 6631034, 6021023, 831446, 98067), keep = c(0.95,
0.95, 0.95, 0.95, 0.9, 0.9), shuffle = c(0.15, 0.1, 0.15, 0.1, 0.15, 0.1))
# This does take a few minutes to run
dom_list <- purrr::pmap(params, function(set, keep, shuffle) {
build_reps(clusters, seed = set, keep_frac = keep, shuffle_frac = shuffle)
})
names(dom_list) <- paste0("P", seq_along(dom_list))We now have six domino() objects that we can treat as
different subjects (P1-6).
linkage_summary classThe linkage_summary() class is the structure used for
the differential signaling workflow. It contains linkages present in
each cluster across the domino() objects it contains, along
with subject metadata. The object contains three slots:
subject_names: a factor with a unique name for each
subject (domino() object) represented in the summary.subject_meta: a data frame with one row per subject
including provided metadata (such as condition, treatment, or other
covariates) that can be used to group or subset subjects for
comparison.subject_linkages: a nested list organized as subject to
cluster to linkage type (tfs for active transcription
factors, rec for active receptors,
incoming_lig for ligands that can activate the cluster,
tfs_rec for transcription factor to receptor pairs, and
rec_lig for receptor to ligand pairs).To create a linkage_summary() object, we use the
summarize_linkages() function, which takes a named list of
domino() objects and a data frame of subject metadata. The
first column of the data frame must contain subject names that match the
names of the domino() objects in the list.
subject_metas <- data.frame(patient = c("P1", "P2", "P3", "P4", "P5", "P6"), time = c("pre",
"post", "pre", "post", "pre", "post"), condition = c("control", "treatment",
"treatment", "control", "control", "control"), batch = c("A", "A", "B", "B",
"C", "C"), stringsAsFactors = FALSE)
link_summary <- summarize_linkages(dom_list, subject_metas, subject_names = subject_metas$patient)
print(link_summary)
#> A linkage summary object of 6 subjects with 4 metadata annotations and linkages between 9 clusters.The print and show methods for the
linkage_summary() class report the number of subjects, the
number of metadata annotations, and the number of clusters present in
the summary.
Since we built this linkage summary from six domino()
objects with six distinct groups and two distinct conditions, we can
subset the summary to include only the subjects of interest.
To create a comparison for only the treatment condition, we can
subset the linkage_summary() object by the relevant
metadata column.
link_treat <- subset(link_summary, subset = condition == "treatment")
print(link_treat)
#> A linkage summary object of 2 subjects with 4 metadata annotations and linkages between 9 clusters.Subjects can also be selected by name or with compound expressions.
subset(link_summary, subset = subject_names %in% c("P1", "P3"))
#> A linkage summary object of 2 subjects with 4 metadata annotations and linkages between 9 clusters.
subset(link_summary, subset = condition == "treatment" | time == "pre")
#> A linkage summary object of 4 subjects with 4 metadata annotations and linkages between 9 clusters.As an alternative to creating separate subset objects, the other
functions in the differential signaling workflow accept a
subject_names argument that can be used to restrict the
subjects considered in the analysis without needing to create a separate
linkage_summary object.
The count_linkage() function tabulates the number of
subjects with a given linkage in a specified cluster, which can also be
broken down by a metadata column. The linkage argument
selects which of the linkage types to count.
If no grouping variable is specified and no subject names are provided, the function will return a data frame with the total number of subjects with each linkage in the selected cluster.
count_all <- count_linkage(link_summary, cluster = "CD8_T_cell", linkage = "rec_lig") %>%
slice_sample(n = 5)| feature | total_count |
|---|---|
| TNFRSF1B <- LTA | 1 |
| CD2 <- CD58 | 6 |
| LAMP1 <- FAM3C | 6 |
| CXCR3 <- CXCL11 | 6 |
| TNFRSF1B <- GRN | 1 |
If specific subjects are of interest, the counts can be restricted to those subjects by providing a vector of the subject names.
count_subjects <- count_linkage(link_summary, cluster = "CD8_T_cell", linkage = "tfs",
subject_names = c("P1", "P2"))| feature | total_count |
|---|---|
| ATF6 | 2 |
| CHURC1 | 1 |
| CREM | 2 |
| DBP | 2 |
| ELF3 | 2 |
The count_linkage() function can also return counts that
are grouped by a metadata column. For example, we can count the linkages
for subjects in the treatment condition grouped by the “time” column to
see how many subjects in each time point have each linkage.
counts_time_trt <- count_linkage(link_treat, cluster = "CD14_monocyte", linkage = "tfs_rec",
group.by = "time")| feature | total_count | post | pre |
|---|---|---|---|
| CEBPE <- TNFRSF1A | 2 | 1 | 1 |
| CEBPA <- CSF1R | 2 | 1 | 1 |
| CEBPE <- CSF1R | 2 | 1 | 1 |
| IRF8 <- TNFRSF1B | 2 | 1 | 1 |
| ATF3 <- TNFRSF1A | 2 | 1 | 1 |
This is equivalent to grouping by the “time” column and providing the subject names of only the subjects in the treatment condition.
counts_time_subj <- count_linkage(link_summary, cluster = "CD14_monocyte", linkage = "tfs_rec",
group.by = "time", subject_names = c("P2", "P3"))
identical(counts_time_trt, counts_time_subj)
#> [1] TRUEOf course, more than one subject can be in each group, and the counts will reflect the number of subjects in each group with each linkage.
counts_condition <- count_linkage(link_summary, cluster = "naive_CD4_T_cell", linkage = "incoming_lig",
group.by = "condition") %>%
slice_sample(n = 5)| feature | total_count | control | treatment |
|---|---|---|---|
| DHCR24 | 6 | 4 | 2 |
| CD70 | 6 | 4 | 2 |
| CDH1 | 6 | 4 | 2 |
| IL7 | 6 | 4 | 2 |
| PTGES2 | 1 | 1 | 0 |
And more than two groups can be compared at once.
counts_batch <- count_linkage(link_summary, cluster = "CD14_monocyte", linkage = "rec",
group.by = "batch") %>%
slice_sample(n = 5)| feature | total_count | A | B | C |
|---|---|---|---|---|
| MARCO | 6 | 2 | 2 | 2 |
| C5AR1 | 6 | 2 | 2 | 2 |
| LTBR | 4 | 2 | 0 | 2 |
| NOTCH2 | 6 | 2 | 2 | 2 |
| CD36 | 6 | 2 | 2 | 2 |
The count_linkage() function tabulates the number of
subjects with a given linkage in a specified cluster, but it does not
test whether the frequency of that linkage is significantly different
across groups. The test_differential_linkages() function
uses a contingency table of presence/absence of each linkage across
groups to perform a Fisher’s exact test. The p-values across all tested
linkages are adjusted for multiple comparisons using the false discovery
rate method.
The test_differential_linkages() function returns a data
frame that includes the cluster, linkage, grouping variable, test name
(always fisher’s exact test), the linkage being compared, the test
statistics (odds ratio, if applicable, and p-value, both raw and
adjusted), total number of subjects where linkage is active, number of
subjects in each group where linkage is active, number of total subjects
compared, and the total number of subjects in each group.
When exactly two groups are compared, the fisher.test()
function returns an odds ratio in addition to the p-value, reflecting
the relative odds of the linkage being present in one group compared to
the other.
dcst_time_trt <- test_differential_linkages(link_treat, cluster = "CD14_monocyte",
linkage = "tfs_rec", group.by = "time")| cluster | linkage | group.by | test_name | feature | odds.ratio | p.value | p.adj | total_count | post_count | pre_count | total_n | post_n | pre_n | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ZNF77 <- CD36 | CD14_monocyte | tfs_rec | time | fishers.exact | ZNF77 <- CD36 | 0 | 1 | 1 | 2 | 1 | 1 | 2 | 1 | 1 |
| FOSB <- CD86 | CD14_monocyte | tfs_rec | time | fishers.exact | FOSB <- CD86 | 0 | 1 | 1 | 2 | 1 | 1 | 2 | 1 | 1 |
| FOSL1 <- TNFRSF1B | CD14_monocyte | tfs_rec | time | fishers.exact | FOSL1 <- TNFRSF1B | 0 | 1 | 1 | 2 | 1 | 1 | 2 | 1 | 1 |
| KLF4 <- FPR1 | CD14_monocyte | tfs_rec | time | fishers.exact | KLF4 <- FPR1 | 0 | 1 | 1 | 2 | 1 | 1 | 2 | 1 | 1 |
| ZNF516 <- FCGR2A | CD14_monocyte | tfs_rec | time | fishers.exact | ZNF516 <- FCGR2A | 0 | 1 | 1 | 2 | 1 | 1 | 2 | 1 | 1 |
As with the count_linkage() function, the same results
can be obtained by providing the subject names of only the subjects in
the treatment condition.
dcst_time_subj <- test_differential_linkages(link_summary, cluster = "CD14_monocyte",
linkage = "tfs_rec", group.by = "time", subject_names = c("P2", "P3"))
identical(dcst_time_trt, dcst_time_subj)
#> [1] TRUEThe test_differential_linkages() function does require
that any grouping variable has at least two groups with at least one
subject in each group.
test_differential_linkages(link_treat, cluster = "CD14_monocyte", linkage = "tfs_rec",
group.by = "condition")
#> Error in `test_differential_linkages()`:
#> ! At least 2 groups of condition must be represented among the included subjects; found: treatmentIt is obviously preferable to have more than one subject in each group (though the function does return results for groups with only one subject). Note that there are no significant differences between the groups in this example, which is expected since the subjects were generated from the same data set.
dcst_time <- test_differential_linkages(link_summary, cluster = "CD8_T_cell", linkage = "rec",
group.by = "time") %>%
dplyr::arrange(desc(odds.ratio))| cluster | linkage | group.by | test_name | feature | odds.ratio | p.value | p.adj | total_count | post_count | pre_count | total_n | post_n | pre_n | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CD27 | CD8_T_cell | rec | time | fishers.exact | CD27 | 0 | 1 | 1 | 6 | 3 | 3 | 6 | 3 | 3 |
| SORL1 | CD8_T_cell | rec | time | fishers.exact | SORL1 | 0 | 1 | 1 | 6 | 3 | 3 | 6 | 3 | 3 |
| CD8A | CD8_T_cell | rec | time | fishers.exact | CD8A | 0 | 1 | 1 | 6 | 3 | 3 | 6 | 3 | 3 |
| CXCR4 | CD8_T_cell | rec | time | fishers.exact | CXCR4 | 0 | 1 | 1 | 6 | 3 | 3 | 6 | 3 | 3 |
| PTGDR | CD8_T_cell | rec | time | fishers.exact | PTGDR | 0 | 1 | 1 | 6 | 3 | 3 | 6 | 3 | 3 |
When more than two groups are compared, the
fisher.test() function does not return an odds ratio (so it
is set to NA in the DCST output), but the p-value still
reflects whether the frequency of the linkage is significantly different
across groups.
dcst_batch <- test_differential_linkages(link_summary, cluster = "naive_CD4_T_cell",
linkage = "rec", group.by = "batch") %>%
arrange(p.value)| cluster | linkage | group.by | test_name | feature | odds.ratio | p.value | p.adj | total_count | A_count | B_count | C_count | total_n | A_n | B_n | C_n | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CD44 | naive_CD4_T_cell | rec | batch | fishers.exact | CD44 | NA | 1 | 1 | 6 | 2 | 2 | 2 | 6 | 2 | 2 | 2 |
| PTGER2 | naive_CD4_T_cell | rec | batch | fishers.exact | PTGER2 | NA | 1 | 1 | 1 | 0 | 0 | 1 | 6 | 2 | 2 | 2 |
| TNFRSF25 | naive_CD4_T_cell | rec | batch | fishers.exact | TNFRSF25 | NA | 1 | 1 | 6 | 2 | 2 | 2 | 6 | 2 | 2 | 2 |
| RORA | naive_CD4_T_cell | rec | batch | fishers.exact | RORA | NA | 1 | 1 | 6 | 2 | 2 | 2 | 6 | 2 | 2 | 2 |
| HLA-C | naive_CD4_T_cell | rec | batch | fishers.exact | HLA-C | NA | 1 | 1 | 6 | 2 | 2 | 2 | 6 | 2 | 2 | 2 |
There is no limit to the number of groups that can be compared at once (as long as there is at least one subject in each group).
dcst_patient <- test_differential_linkages(link_summary, cluster = "naive_CD4_T_cell",
linkage = "rec_lig", group.by = "patient") %>%
arrange(p.value)| cluster | linkage | group.by | test_name | feature | odds.ratio | p.value | p.adj | total_count | P1_count | P2_count | P3_count | P4_count | P5_count | P6_count | total_n | P1_n | P2_n | P3_n | P4_n | P5_n | P6_n | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| IL27_receptor <- IL27 | naive_CD4_T_cell | rec_lig | patient | fishers.exact | IL27_receptor <- IL27 | NA | 1 | 1 | 4 | 1 | 1 | 1 | 1 | 0 | 0 | 6 | 1 | 1 | 1 | 1 | 1 | 1 |
| NCR3 <- NCR3LG1 | naive_CD4_T_cell | rec_lig | patient | fishers.exact | NCR3 <- NCR3LG1 | NA | 1 | 1 | 5 | 1 | 1 | 1 | 0 | 1 | 1 | 6 | 1 | 1 | 1 | 1 | 1 | 1 |
| ICOS <- TNF | naive_CD4_T_cell | rec_lig | patient | fishers.exact | ICOS <- TNF | NA | 1 | 1 | 5 | 0 | 1 | 1 | 1 | 1 | 1 | 6 | 1 | 1 | 1 | 1 | 1 | 1 |
| LAMP1 <- FAM3C | naive_CD4_T_cell | rec_lig | patient | fishers.exact | LAMP1 <- FAM3C | NA | 1 | 1 | 6 | 1 | 1 | 1 | 1 | 1 | 1 | 6 | 1 | 1 | 1 | 1 | 1 | 1 |
| NCR3 <- BAG6 | naive_CD4_T_cell | rec_lig | patient | fishers.exact | NCR3 <- BAG6 | NA | 1 | 1 | 5 | 1 | 1 | 1 | 0 | 1 | 1 | 6 | 1 | 1 | 1 | 1 | 1 | 1 |
The DCST paper includes additional information on the recommended number of subjects for robust results (ideally at least 15 per group), as well as an approach to simulate biological replicates using bootstrapping of cells for experiments with small numbers of replicates (scripts for the manuscript can be found at the corresponding GitHub repository).
In order to visualize the results of the DCST, the
plot_differential_linkages() function can be used to create
a heatmap of the test statistic along with row-level bar plots showing
the proportions of subjects in which linkages are active.
The function has arguments for the test statistic (the column of the DCST result that contains the statistic to be plotted), the range of values to be plotted, the heatmap gradient, and the color palette to use for the groups being compared.
An example of a two-group comparison is shown below, where the
p-value is used as the test statistic and rows are ordered by ascending
p-value. All p-values are 1 in this example, so every cell shows the
maximum gradient color.
plot_differential_linkages(differential_linkages = slice_sample(dcst_time_trt, n = 15),
test_statistic = "p.value", stat_range = c(0, 1), stat_ranking = "ascending")plot_differential_linkages(differential_linkages = slice_sample(dcst_time, n = 15),
test_statistic = "p.value", gradient = c(minimum = "#c44f4f", maximum = "#ffd771"),
stat_ranking = "ascending")If comparing between two groups, it is also possible to use the odds ratio as the test statistic, which will show the relative odds of a linkage being present in one group compared to the other.
plot_differential_linkages(differential_linkages = slice_sample(dcst_time, n = 15),
test_statistic = "odds.ratio", stat_range = c(0, 4), stat_ranking = "descending",
gradient = c(minimum = "gray90", maximum = "orange"))When comparing between more than two groups, the odds ratio is not available, so the p-value (or adjusted p-value) is the only option for the test statistic.
plot_differential_linkages(differential_linkages = slice_sample(dcst_batch, n = 15),
test_statistic = "p.value", stat_range = c(0, 1), stat_ranking = "ascending",
gradient = c(minimum = "#ff4c43", maximum = "#ffe6c0"))A custom color can be assigned to each group’s barplot annotation
using the group_palette argument, which takes a named
vector of colors where names correspond to the group names in the DCST
result.
plot_differential_linkages(differential_linkages = slice_sample(dcst_batch, n = 15),
test_statistic = "p.value", stat_range = c(0, 1), stat_ranking = "ascending",
gradient = c(minimum = "#7b43ff", maximum = "#c0c5ff"), group_palette = c(A = "#64a5e6",
B = "#71fcad", C = "#c968c3"))This can be useful for consistency across multiple plots or to match publication color schemes.
Here’s an example of what a workflow might look like from the list of
domino() objects to the final visualization of differential
linkages. We begin with the list of synthetic domino()
objects created above, create the metadata, create a
link_summary() object, examine the counts for the
linkages,and then run the DCST for a comparison by batch. In this
example, we focus on incoming ligands in the CD14 monocyte cluster.
subject_metas <- data.frame(patient = c("P1", "P2", "P3", "P4", "P5", "P6"), time = c("pre",
"post", "pre", "post", "pre", "post"), condition = c("control", "treatment",
"treatment", "control", "control", "control"), batch = c("A", "A", "B", "B",
"C", "C"), stringsAsFactors = FALSE)
link_summary <- summarize_linkages(dom_list, subject_metas, subject_names = subject_metas$patient)
counts_batch <- count_linkage(link_summary, cluster = "CD14_monocyte", linkage = "incoming_lig",
group.by = "batch")
dcst_batch <- test_differential_linkages(link_summary, cluster = "CD14_monocyte",
linkage = "incoming_lig", group.by = "batch")
plot_differential_linkages(dcst_batch, test_statistic = "p.value", stat_ranking = "ascending",
gradient = c(minimum = "#e02c2c", maximum = "#feffcd"), group_palette = c(A = "#2a6aaa",
B = "#07bd56", C = "#a244a3"))Other methods, such as shuffling all cluster labels, results in clusters that are all statistically similar which results in few linkages reaching significance. While this approach reduces the differences between domino objects, it preserves enough signaling to infer a reasonable number of (insignificant) linkages.↩︎
These parameters are fairly permissive to ensure enough linkages are inferred for demonstration purposes.↩︎