Differential Expression Analysis for CCC & Downstream Signalling Networks#
Background#
Cell-cell communication (CCC) events play a critical role in diseases, often experiencing deregulation. To identify differential expression of CCC events between conditions, we can build upon standard differential expression analysis (DEA) approaches, such as DESeq2. While dimensionality reduction methods like extracting intercellular programmes with MOFA+ and Tensor-cell2cell reduce CCC into sets of loadings, hypothesis-driven DEA tests focus on individual gene changes, making them easier to understand and interpret.
In this tutorial, we perform DEA at the pseudobulk level to assess differential expression of genes between conditions. We then translate the results into deregulated complex-informed ligand-receptor interactions and analyze their connections to downstream signaling events.
For further information on pseudobulk DEA, please refer to the Differential Gene Expression chapter in the Single-cell Best Practices book, as well as Decoupler’s pseudobulk vignette. These resources provide more comprehensive details on the subject.
Load Packages#
import numpy as np
import pandas as pd
import scanpy as sc
import plotnine as p9
import liana as li
import decoupler as dc
import omnipath as op
# Import DESeq2
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats
# embed inline figures as JPEG (not the 2x 'retina' PNG default) to keep the notebook small
sc.set_figure_params(ipython_format='jpeg')
# Obtain TF regulons
net = dc.op.collectri(organism='human', remove_complexes=False, license='academic', verbose=False)
Load & Prep Data#
As a simple example, we will look at ~25k PBMCs from 8 pooled patient lupus samples, each before and after IFN-beta stimulation (Kang et al., 2018; GSE96583). Note that by focusing on PBMCs, for the purpose of this tutorial, we assume that coordinated events occur among them.
This dataset is downloaded from a link on Figshare; preprocessed for pertpy.
adata = li.ds.kang_2018()
adata
AnnData object with n_obs × n_vars = 24673 × 15706
obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'condition', 'cluster', 'cell_type', 'patient', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'sample', 'cell_abbr'
var: 'name'
obsm: 'X_pca', 'X_umap'
layers: 'counts'
Define columns of interest from .obs
Note that we use cell abbreviations because MOFA will use them as labels for the views.
sample_key = 'sample'
groupby = 'cell_abbr'
condition_key = 'condition'
Basic QC#
Note that this data has been largely pre-processed & annotated, we refer the user to the Quality Control and other relevant chapters from the best-practices book for information about pre-processing and annotation steps.
# filter cells and genes
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
Showcase the data#
# Show pre-computed UMAP
sc.pl.umap(adata, color=[condition_key, sample_key, 'cell_type', groupby], frameon=False, ncols=2)
Differential Testing#
First, we need to generate pseudobulk profiles for each cell type, and we do so using the decoupler package.
pdata = dc.pp.pseudobulk(
adata,
sample_col=sample_key,
groups_col=groupby,
layer='counts',
mode='sum'
)
pdata
AnnData object with n_obs × n_vars = 128 × 15701
obs: 'sample', 'cell_abbr', 'condition', 'cell_type', 'patient', 'psbulk_cells', 'psbulk_counts'
var: 'name', 'n_cells'
layers: 'psbulk_props'
# filter samples based on number of cells and counts
dc.pp.filter_samples(pdata, min_cells = 10, min_counts=1000)
We can plot the quality control metrics for each pseudobulk sample:
dc.pl.filter_samples(pdata, groupby=[sample_key, groupby], figsize=(11, 4))
Differential Expression Analysis#
Next, now that we have generated the pseudobulk profiles, we can perform some edgeR-like filtering using decoupler-py, and then differential expression analysis using the pydeseq2 package - a re-implementation of the original DESeq2 method (Love et al., 2014).
Here, we perform DEA on the pseudobulk profiles for each cell type, for more info check this tutorial: https://decoupler.readthedocs.io/en/latest/notebooks/scell/rna_psbk.html
%%capture
dea_results = {}
quiet = True
for cell_group in pdata.obs[groupby].unique():
# Select cell profiles
ctdata = pdata[pdata.obs[groupby] == cell_group].copy()
# Obtain genes that pass the edgeR-like thresholds
# NOTE: QC thresholds might differ between cell types, consider applying them by cell type
genes = dc.pp.filter_by_expr(ctdata,
group=condition_key,
min_count=5, # a minimum number of counts in a number of samples
min_total_count=10, # a minimum total number of reads across samples
inplace=False,
)
# Filter by these genes
ctdata = ctdata[:, genes].copy()
# Build DESeq2 object
# NOTE: this data is actually paired, so one could consider fitting the patient label as a confounder
dds = DeseqDataSet(
adata=ctdata,
design_factors=condition_key,
ref_level=[condition_key, 'ctrl'], # set control as reference
refit_cooks=True,
quiet=quiet
)
# Compute LFCs
dds.deseq2()
# Contrast between stim and ctrl
stat_res = DeseqStats(dds, contrast=[condition_key, 'stim', 'ctrl'], quiet=quiet)
stat_res.quiet = quiet
# Compute Wald test
stat_res.summary()
# Shrink LFCs
stat_res.lfc_shrink(coeff='condition[T.stim]') # {condition_key}_cond_vs_ref
dea_results[cell_group] = stat_res.results_df
This results in a wall of currently unavoidable verbose text and prints, as such I use %%capture to hide it.
One can use quiet to some of the functions but not logfc_shrinkage
# concat results across cell types
dea_df = pd.concat(dea_results)
dea_df = dea_df.reset_index().rename(columns={'level_0': groupby,'level_1':'index'}).set_index('index')
dea_df.head()
| cell_abbr | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|---|
| index | |||||||
| NOC2L | B | 13.944414 | -0.269461 | 0.215233 | -1.393996 | 1.633188e-01 | 3.865297e-01 |
| ISG15 | B | 687.496723 | 5.529444 | 0.154587 | 35.670704 | 1.125446e-278 | 3.457370e-275 |
| TNFRSF18 | B | 11.700916 | -1.144321 | 0.278038 | -4.437036 | 9.120602e-06 | 9.216608e-05 |
| SDF4 | B | 8.173816 | -0.476321 | 0.282784 | -1.948320 | 5.137671e-02 | 1.709959e-01 |
| UBE2J2 | B | 8.212013 | -0.365051 | 0.296765 | -1.450806 | 1.468339e-01 | 3.632300e-01 |
# PyDeseq Seems to intrdoce NAs for some p-values
# NOTE: there sometimes some NaN being introduced, best to double check that, in this case it's only for a single gene, but it might be a problem.
len(dea_df[dea_df.isna().any(axis=1)])
1
DEA to Ligand-Receptor Interactions#
Now that we have DEA results per gene, we can combine them into statistics of potentially deregulated ligand-receptor interactions.
To do so, liana provides a simple function li.mt.df_to_lr that calculates average expression as well as proportions based on the passed adata object, and combines those with the DEA results and a ligand-receptor resource. Since in this case we want to focus on gene statics relevant to the condition (stim), let’s subset the adata to those and normalize the counts.
adata = adata[adata.obs[condition_key]=='stim'].copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
Let’s combine the DEA results with the ligand-receptor interactions. We need to pass the names of the statistics from the DEA table in which we are interest to li.mt.df_to_lr, here we will use the adjusted p-values and Wald test statistic.
lr_res = li.mt.df_to_lr(adata,
dea_df=dea_df,
resource_name='consensus', # NOTE: uses HUMAN gene symbols!
expr_prop=0.1, # calculated for adata as passed - used to filter interactions
groupby=groupby,
stat_keys=['stat', 'pvalue', 'padj'],
use_raw=False,
complex_col='stat', # NOTE: we use the Wald Stat to deal with complexes
verbose=True,
return_all_lrs=False,
)
lr_res = lr_res.sort_values("interaction_stat", ascending=False, key=abs)
lr_res.head()
| ligand | receptor | ligand_complex | receptor_complex | source | ligand_stat | ligand_pvalue | ligand_padj | ligand_expr | ligand_props | ... | receptor_pvalue | receptor_padj | receptor_expr | receptor_props | interaction_stat | interaction_pvalue | interaction_padj | interaction_expr | interaction_props | interaction | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 218 | CXCL11 | CCR5 | CXCL11 | CCR5 | CD14 | 27.347028 | 1.171470e-164 | 1.223483e-161 | 3.106561 | 0.899096 | ... | 2.014382e-43 | 6.896889e-42 | 0.601853 | 0.332369 | 20.581987 | 1.007191e-43 | 3.448445e-42 | 1.854207 | 0.615732 | CXCL11^CCR5 |
| 239 | CCL8 | CCR5 | CCL8 | CCR5 | CD14 | 26.495511 | 1.091828e-154 | 8.145038e-152 | 4.313416 | 0.894756 | ... | 2.014382e-43 | 6.896889e-42 | 0.601853 | 0.332369 | 20.156228 | 1.007191e-43 | 3.448445e-42 | 2.457634 | 0.613562 | CCL8^CCR5 |
| 238 | CCL8 | CCR1 | CCL8 | CCR1 | CD14 | 26.495511 | 1.091828e-154 | 8.145038e-152 | 4.313416 | 0.894756 | ... | 4.931677e-35 | 1.300668e-33 | 1.079367 | 0.552260 | 19.422259 | 2.465839e-35 | 6.503338e-34 | 2.696391 | 0.723508 | CCL8^CCR1 |
| 480 | CCL8 | CCR1 | CCL8 | CCR1 | CD14 | 26.495511 | 1.091828e-154 | 8.145038e-152 | 4.313416 | 0.894756 | ... | 3.790010e-07 | 3.235442e-06 | 0.384328 | 0.247232 | 15.787362 | 1.895005e-07 | 1.617721e-06 | 2.348872 | 0.570994 | CCL8^CCR1 |
| 459 | CXCL11 | CCR5 | CXCL11 | CCR5 | CD14 | 27.347028 | 1.171470e-164 | 1.223483e-161 | 3.106561 | 0.899096 | ... | 2.868568e-05 | 1.829894e-04 | 0.309305 | 0.217712 | 15.765342 | 1.434284e-05 | 9.149472e-05 | 1.707933 | 0.558404 | CXCL11^CCR5 |
5 rows × 22 columns
Dealing with heteromeric complexes#
LIANA will filter lowly-expressed interactions, i.e. those for which any of the genes are not expressed in at least 0.1 of the cells (by default) in the AnnData object. This can be adjusted with the expr_prop parameter.
Moreover, to deal with complexes for each cell type, as either source or target of the potential CCC events, LIANA will find and assign the subunit of a complex with the lowest gene expression (by default) as the subunit of interest, and will then use the stats for that subunit as the stats of the whole protein complex.
To this end, we also provide the option to provide a complex_col parameter, which will be used to assign the complex subunit of interest. This column should be a part of the stat_keys. Note that the absolute minimum value is used (i.e. the value closest to 0 is thought to be the ‘worst’ result), so this will not work for statistics with ascending values (e.g. p-values).
Visualize the Results#
interaction_* columns returned by li.mt.df_to_lr are just the mean of the ligand and receptor columns of the corresponding statistic!
Please use with caution as this is just a summary of the interaction that we can use to e.g. to sort the interactions as done above.
Instead, we recommend to use the ligand and receptor statistics separately to filter and visualize the interactions.
Moreover, by averaging the statistics across the ligand and receptor, we are focusing on the interactions for which both the ligand and receptor are deregulated in the same direction, i.e. both up or both down. However, this might ignore interactions in which e.g. the the ligand is deregulated while the receptor is not, or such where they are deregulated in opposite directions. These could represent potential inhibitory mechanisms, but we leave this to the user to explore.
# Let's visualize how this looks like for all interactions (across all cell types)
lr_res = lr_res.sort_values("interaction_stat", ascending=False)
lr_res['interaction_stat'].hist(bins=50)
<Axes: >
Now that we have covered the basics, we can visualize our interactions in a few ways.
Let’s start with the top interactions according to their Wald statistic, and then plot the statistics for the ligands & receptors involved in those interactions across cell types, to do so LIANA+ provide li.pl.tileplot:
li.pl.tileplot(liana_res=lr_res,
fill = 'expr',
label='padj',
label_fn = lambda x: '*' if x < 0.05 else np.nan,
top_n=15,
orderby = 'interaction_stat',
orderby_ascending = False,
orderby_absolute = False,
source_title='Ligand',
target_title='Receptor',
)
If you want to plot the expression values for ligand-receptor interactions without the DEA statistics, you can set the return_all_lrs parameter to True in the li.mt.df_to_lr function. This will return a dataframe with all the ligand-receptor interactions, where missing DEA stats will be set as nan, while mean expression and proportions per cluster will be obtained via the AnnData object.
Ligand-Receptor Plot#
We can also use visualize of the stats, summarized at the level of the interaction, to prioritize the interactions, or any subunit statistics using li.pl.dotplot. For example, we can visualize the mean Wald statistic between the ligand & receptor, together with the pvalues for the ligand.
plot = li.pl.dotplot(liana_res=lr_res,
colour='interaction_stat',
size='ligand_pvalue',
inverse_size=True,
orderby='interaction_stat',
orderby_ascending=False,
orderby_absolute=True,
top_n=10,
size_range=(0.5, 4)
)
# customize plot
(
plot
+ p9.theme_bw(base_size=14)
# fill cmap blue to red, with 0 the middle
+ p9.scale_color_cmap('RdBu_r', limits=(-10, 10))
# rotate x
+ p9.theme(axis_text_x=p9.element_text(angle=90), figure_size=(11, 6))
)
Now that we have identified a set of interactions that are potentially deregulated we can look into the downstream signalling events that they might be involved in.
Intracellular Signaling Networks#
Cellular signaling networks govern the behavior of cells, allowing them to respond to external signals, including various cell-cell communication events. Thus, CCC events can be thought of as upstream perturbants of intracellular signaling networks that lead to deregulations of downstream signaling events. Such deregulations are expected to be associated with various conditions and disease. Thus, understanding intracellular signaling networks is critical to model the cellular mechanisms.
Here, we will combine several tools to identify plausible signalling cascades driven by CCC events.
Our approach includes the following steps:
Select a number of potentially deregulated ligand-receptor interactions (input nodes), in terms of summarized PyDESeq2 statistics.
Select a number of potentially deregulated TFs (output nodes). This is done via the use of Transcription factor (TF) activity inference. Carried out on differential gene expression data using TF regulon knowledge with decoupler
Obtain a prior knowledge network (PKNs), with signed protein-protein interactions from OmniPath.
Generate weights for the nodes in the PKN
Use CORNETO to identify a solution in the form of a causal (smallest sign-consistent signaling) network that explains the measured inputs and outputs
Import OmniPath#
For this part OmniPath is required.
# utily function to select top n interactions
def select_top_n(d, n=None):
d = dict(sorted(d.items(), key=lambda item: abs(item[1]), reverse=True))
return {k: v for i, (k, v) in enumerate(d.items()) if i < n}
Select Cell types of Interest#
One limitation of using DEA to identify interactions of interest is that it tells us little about deregulation at the level of cell types. However, from dimensionality reductions on CCC, as done with Tensor-cell2cell & MOFA on the same dataset, we can see there is a potential deregulation of CCC that involve CD14 monocytes both as sources (senders) and targets (or receivers) of intecellular communication. Thus, we will focus on the interactions and downstream signalling within that cell type.
source_label = 'CD14'
target_label = 'CD14'
# NOTE: We sort by the absolute value of the interaction stat
lr_stats = lr_res[lr_res['source'].isin([source_label]) & lr_res['target'].isin([target_label])].copy()
lr_stats = lr_stats.sort_values('interaction_stat', ascending=False, key=abs)
Select Receptors based on interaction stats#
These will be used as the input or start nodes for the network. In this case, we will use interactions potentially involved in autocrine signalling in CD14 monocytes.
lr_dict = lr_stats.set_index('receptor')['interaction_stat'].to_dict()
input_scores = select_top_n(lr_dict, n=10)
input_scores
{'CD40': 15.288575765820969,
'CD80': 11.76003734068972,
'SIRPA': 9.483281090234765,
'CCR1': 8.833055306399888,
'HLA-DPB1': 8.829743317474247,
'LILRB2': 7.633845587613143,
'LILRB1': 5.81618655584982,
'CCR5': 5.734035152179016,
'LILRA3': 5.2047813052906795,
'CD47': 5.050045257819155}
Select Transcription Factors of interest#
Before we select the transcription factors, we need to infer their activity. We will do so using decoupler with CollecTri regulons. Specifically, we will estimate TF activities using the Wald statistics (from PyDESeq2) for the genes in the regulons.
# First, let's transform the DEA statistics into a DF
# we will use these to estimate deregulated TF activity
dea_wide = dea_df[[groupby, 'stat']].reset_index(names='genes').pivot(index=groupby, columns='genes', values='stat')
dea_wide = dea_wide.fillna(0)
dea_wide
| genes | A1BG | AAAS | AAED1 | AAGAB | AAK1 | AAMDC | AAMP | AAR2 | AARS | AARSD1 | ... | ZSCAN16 | ZSCAN16-AS1 | ZSCAN18 | ZSCAN32 | ZSWIM6 | ZSWIM7 | ZUFSP | ZW10 | ZYX | ZZZ3 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| cell_abbr | |||||||||||||||||||||
| B | 0.000000 | 0.000000 | -0.487284 | 0.000000 | 0.000000 | 0.000000 | -1.906085 | 0.000000 | 0.556648 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | -1.891232 | 0.000000 | 0.000000 | 0.000000 | 0.00000 |
| CD14 | -2.990605 | 0.000000 | 9.714917 | 0.000000 | 1.076383 | 6.880305 | -2.973870 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | -1.028499 | -4.063465 | 1.667606 | 0.000000 | -1.732015 | 0.00000 |
| CD4T | -2.235063 | -1.901853 | 0.187765 | 1.047908 | -4.141667 | -0.094972 | -0.505706 | 1.956408 | 1.763777 | -0.262918 | ... | -0.194163 | -0.947079 | -1.499626 | -0.375836 | 0.000000 | -1.637439 | 3.342277 | -0.157619 | -1.456379 | -0.24172 |
| CD8T | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 |
| DCs | 0.000000 | 0.000000 | 3.171024 | 0.000000 | 0.000000 | 0.000000 | -1.614708 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | -2.563751 | 0.00000 |
| FGR3 | 0.000000 | 0.000000 | 2.790252 | 0.000000 | 0.000000 | 0.000000 | -3.204464 | 0.000000 | 0.000000 | 0.000000 | ... | 4.118250 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | -2.077878 | 0.00000 |
| NK | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | -2.005801 | 0.000000 | 0.000000 | 0.000000 | ... | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 |
7 rows × 7802 columns
# Run Enrichment Analysis
estimates, pvals = dc.mt.ulm(dea_wide, net=net)
estimates.T.sort_values(target_label, key=abs, ascending=False).head()
| B | CD14 | CD4T | CD8T | DCs | FGR3 | NK | |
|---|---|---|---|---|---|---|---|
| STAT1 | 18.516348 | 14.027463 | 17.311885 | 18.165816 | 18.266293 | 15.942676 | 17.198443 |
| IRF1 | 15.936464 | 11.398815 | 15.216455 | 13.256008 | 15.011091 | 12.502367 | 13.513206 |
| STAT2 | 15.379482 | 10.497779 | 12.846587 | 14.634710 | 15.763186 | 13.113093 | 14.014088 |
| IRF9 | 19.410173 | 9.332927 | 15.766561 | 20.848557 | 17.578566 | 14.507459 | 18.989734 |
| NFKB1 | 6.726540 | 7.847975 | 8.384435 | 5.889624 | 9.236922 | 8.669230 | 7.108677 |
Select top TFs#
Now that we have the potentially deregulated TFs, we focus on the top 10 TFs, based on their enrichment scores. In this case, we will look specifically at the top TFs deregulated in CD14 monocytes.
tf_data = estimates.copy()
tf_dict = tf_data.loc[target_label].to_dict()
output_scores = select_top_n(tf_dict, n=5)
Generate a Prior Knowledge Network#
Now we will obtain protein-protein interactions from OmniPath, filter them according to curation effort to ensure we only keep those that are of high quality, and convert them into a knowledge graph.
# obtain ppi network
ppis = op.interactions.OmniPath().get(genesymbols = True)
ppis['mor'] = ppis['is_stimulation'].astype(int) - ppis['is_inhibition'].astype(int)
ppis = ppis[(ppis['mor'] != 0) & (ppis['curation_effort'] >= 5) & ppis['consensus_direction']]
input_pkn = ppis[['source_genesymbol', 'mor', 'target_genesymbol']]
input_pkn.columns = ['source', 'mor', 'target']
input_pkn.head()
| source | mor | target | |
|---|---|---|---|
| 3 | CAV1 | 1 | TRPC1 |
| 6 | ITPR2 | 1 | TRPC1 |
| 9 | STIM1 | 1 | TRPC1 |
| 10 | TRPC1 | 1 | TRPC3 |
| 11 | TRPC3 | 1 | TRPC1 |
# convert the PPI network into a knowledge graph
prior_graph = li.rs.build_prior_network(input_pkn, input_scores, output_scores, verbose=True)
In this section we use Prior Knowledge Networks (PKNs) from OmniPath to generate network hypotheses based on the deregulated interactions considering both sign and direction. Specifically, we focus on highly curated protein-protein interactions, which often represent hubs in the network. Since such network approaches are highly dependent on prior knowledge, for a review on prior knowledge bias and similar network inference methods, including thier limitations, see Garrido-Rodriguez et al., 2022.
Calculate Node weights#
Calculate gene expression proportions within the target cell type; we will use those as node weights in the network.
temp = adata[adata.obs[groupby] == target_label].copy()
node_weights = pd.DataFrame(temp.X.getnnz(axis=0) / temp.n_obs, index=temp.var_names)
node_weights = node_weights.rename(columns={0: 'props'})
node_weights = node_weights['props'].to_dict()
Find Causal Network#
CORNETO (Rodriguez-Mier et al., In prep) generalizes biological network inference problems using convex and combinatorial optimization. Here, we use it to find the smallest sign-consistent network that explains the measured inputs and outputs, a network inference problem formulated in CARNIVAL.
To run CORNETO, we need to first install it; it’s very lightweight and can be installed via pip:
pip install corneto==0.9.1-alpha.6 cvxpy cylp
import corneto as cn
cn.info()
Installed version: v0.9.1-alpha.6 (latest: v0.0.0)
Available backends: CVXPY v1.9.2
Default backend (corneto.K): CVXPY
Installed solvers: CLARABEL, SCS, SCIP, SCIPY, HIGHS, OSQP
Graphviz version: v0.21
Repository: https://github.com/saezlab/corneto
df_res, problem = li.mt.find_causalnet(
prior_graph,
input_scores,
output_scores,
node_weights,
# penalize (max_penalty) nodes with counts in less than 0.1 of the cells
node_cutoff=0.1,
max_penalty=1,
# the penaly of those in > 0.1 prop of cells set to:
min_penalty=0.01,
edge_penalty=0.1,
verbose=True,
max_runs=1, # NOTE that this repeats the solving either until the max runs are reached (increase to 10+ for stable results)
stable_runs=1, # or until X number of consequitive stable runs are reached (i.e. no new edges are added, increase to 50+ for stable results)
solver='HIGHS' # free MILP solver bundled via corneto/cvxpy
)
===============================================================================
CVXPY
v1.9.2
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
Running HiGHS 1.15.1 (git hash: 04024d7): Copyright (c) 2026 under MIT licence terms
Includes third-party software components, see THIRD_PARTY_NOTICES.md for full details
MIP has 43301 rows; 11130 cols; 78175 nonzeros; 8195 integer variables (8195 binary)
Coefficient ranges:
Matrix [1e-03, 6e+02]
Cost [5e-03, 2e+01]
Bound [1e+00, 1e+00]
RHS [1e-03, 6e+02]
Presolving model
17855 rows, 10083 cols, 49087 nonzeros 0s
16873 rows, 9766 cols, 51302 nonzeros 0s
Presolve reductions: rows 16873(-26428); columns 9766(-1364); nonzeros 51302(-26873)
Solving MIP model with:
16873 rows
9766 cols (7265 binary, 0 integer, 0 implied int., 2501 continuous, 0 domain fixed)
51302 nonzeros
Thread count 32 (of 64 threads). Using 1 max workers. Parallel search off
Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic;
I => Shifting; J => Feasibility jump; L => Sub-MIP; P => Empty MIP; R => Randomized rounding;
S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution; Y => HiGHS solution;
Z => ZI Round; l => Trivial lower; p => Trivial point; u => Trivial upper; z => Trivial zero
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Src Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
J 0 0 0 0.00% -inf -16.44316752 Large 0 0 0 0 0.3s
0 0 0 0.00% -81.08915704 -16.44316752 393.15% 0 0 14 1238 0.4s
L 0 0 0 0.00% -80.0103794 -78.16629512 2.36% 108 22 140 3487 2.4s
L 0 0 0 0.00% -80.0103694 -78.36629512 2.10% 130 24 140 9752 3.5s
6.3% inactive integer columns, restarting
Model after restart has 12811 rows, 8697 cols (6321 bin., 0 int., 0 impl., 2376 cont., 0 dom.fix.), and 38953 nonzeros
0 0 0 0.00% -80.0103694 -78.36629512 2.10% 19 0 0 14413 3.7s
0 0 0 0.00% -80.0103694 -78.36629512 2.10% 19 19 6 15856 3.8s
L 0 0 0 0.00% -79.71388778 -78.60804291 1.41% 280 46 6 18153 5.6s
Symmetry detection completed in 0.1s
Found 3 generator(s) and 104 full orbitope(s) acting on 219 columns
12 0 2 50.20% -79.71388778 -78.60804291 1.41% 298 36 149 119231 17.7s
31 0 11 50.83% -79.71388778 -78.60804291 1.41% 336 36 884 150703 22.9s
89 6 37 54.39% -79.68203692 -78.60804291 1.37% 458 60 1264 171320 28.1s
152 11 62 57.52% -79.68203692 -78.60804291 1.37% 576 69 1722 195050 33.5s
225 14 93 63.69% -79.60700696 -78.60804291 1.27% 634 63 2219 222333 38.9s
278 20 120 66.90% -79.60700696 -78.60804291 1.27% 725 78 2577 245391 44.0s
315 20 137 66.91% -79.60700696 -78.60804291 1.27% 745 78 2839 270837 49.0s
349 25 155 66.91% -79.60627978 -78.60804291 1.27% 764 78 4096 295787 54.1s
394 19 178 71.94% -79.60627978 -78.60804291 1.27% 886 84 4502 313956 59.1s
L 400 19 181 77.80% -79.6052466 -78.60804291 1.27% 1233 218 4629 333127 73.4s
428 19 194 79.26% -79.6052466 -78.60804291 1.27% 1267 113 4894 369355 80.6s
451 9 211 89.03% -79.55235237 -78.60804291 1.20% 1215 141 5158 389130 85.8s
470 0 226 100.00% -78.61513288 -78.60804291 0.01% 1212 151 5300 400551 88.6s
Solving report
Status Optimal
Primal bound -78.6080429078
Dual bound -78.615132883
Gap 0.00902% (tolerance: 0.01%)
P-D integral 9.42142851887
Solution status feasible
-78.6080429078 (objective)
0 (bound viol.)
3.31401572851e-13 (int. viol.)
0 (row viol.)
Timing 88.61
0.86 (Presolve)
MIP time [calls] = 0.17 [1]
subMIP time [calls] = 0.69 [82]
87.73 (Solve)
MIP time [calls] = 79.58 [1]
subMIP time [calls] = 8.16 [70]
0.00 (Postsolve)
MIP time [calls] = 0.00 [1]
subMIP time [calls] = 0.00 [82]
Max sub-MIP depth 16
Nodes 470
Repair LPs 0
LP iterations 400551
181533 (strong br.)
18793 (separation)
35049 (heuristics)
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
===============================================================================
CVXPY
v1.9.2
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
Running HiGHS 1.15.1 (git hash: 04024d7): Copyright (c) 2026 under MIT licence terms
Includes third-party software components, see THIRD_PARTY_NOTICES.md for full details
MIP has 43301 rows; 11130 cols; 78175 nonzeros; 8195 integer variables (8195 binary)
Coefficient ranges:
Matrix [1e-03, 6e+02]
Cost [5e-03, 2e+01]
Bound [1e+00, 1e+00]
RHS [1e-03, 6e+02]
Presolving model
17855 rows, 10083 cols, 49087 nonzeros 0s
16870 rows, 9763 cols, 51330 nonzeros 0s
Presolve reductions: rows 16870(-26431); columns 9763(-1367); nonzeros 51330(-26845)
Solving MIP model with:
16870 rows
9763 cols (7265 binary, 0 integer, 0 implied int., 2498 continuous, 0 domain fixed)
51330 nonzeros
Thread count 32 (of 64 threads). Using 1 max workers. Parallel search off
Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic;
I => Shifting; J => Feasibility jump; L => Sub-MIP; P => Empty MIP; R => Randomized rounding;
S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution; Y => HiGHS solution;
Z => ZI Round; l => Trivial lower; p => Trivial point; u => Trivial upper; z => Trivial zero
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Src Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
J 0 0 0 0.00% -inf -16.80132942 Large 0 0 0 0 0.5s
0 0 0 0.00% -81.08915704 -16.80132942 382.64% 0 0 14 1153 0.6s
L 0 0 0 0.00% -80.0103794 -77.35652091 3.43% 177 25 200 3220 3.0s
L 0 0 0 0.00% -80.0103694 -78.60057505 1.79% 190 27 200 8706 5.1s
6.3% inactive integer columns, restarting
Model after restart has 12808 rows, 8694 cols (6321 bin., 0 int., 0 impl., 2373 cont., 0 dom.fix.), and 38981 nonzeros
0 0 0 0.00% -80.01036218 -78.60057505 1.79% 23 0 0 16532 5.4s
0 0 0 0.00% -80.01036218 -78.60057505 1.79% 23 23 6 18065 5.5s
Symmetry detection completed in 0.1s
Found 3 generator(s) and 104 full orbitope(s) acting on 219 columns
B 0 0 0 0.00% -79.7037422 -78.60804291 1.39% 580 47 37 21569 10.8s
T 14 2 4 2.15% -79.7037422 -78.60804291 1.39% 605 47 208 52431 12.1s
22 3 8 45.90% -79.7037422 -78.60804291 1.39% 618 47 292 97650 19.1s
34 7 12 49.41% -79.7037422 -78.60804291 1.39% 658 47 535 131355 24.6s
103 24 37 50.31% -79.65373459 -78.60804291 1.33% 716 50 1141 160629 29.9s
167 18 70 56.36% -79.54435133 -78.60804291 1.19% 1125 79 1625 182288 34.9s
249 12 113 61.80% -79.54435133 -78.60804291 1.19% 1063 54 2150 200508 39.9s
356 18 163 68.25% -79.53854833 -78.60804291 1.18% 1374 102 2500 222313 45.0s
423 4 203 93.55% -79.49550468 -78.60804291 1.13% 1359 78 2970 247508 50.2s
431 0 208 100.00% -78.61428762 -78.60804291 0.01% 1413 74 3289 257262 52.4s
Solving report
Status Optimal
Primal bound -78.6080429078
Dual bound -78.6142876235
Gap 0.00794% (tolerance: 0.01%)
P-D integral 10.2833460324
Solution status feasible
-78.6080429078 (objective)
0 (bound viol.)
2.72243030299e-12 (int. viol.)
0 (row viol.)
Timing 52.45
0.53 (Presolve)
MIP time [calls] = 0.32 [1]
subMIP time [calls] = 0.21 [15]
51.92 (Solve)
MIP time [calls] = 48.66 [1]
subMIP time [calls] = 3.26 [15]
0.00 (Postsolve)
MIP time [calls] = 0.00 [1]
subMIP time [calls] = 0.00 [15]
Max sub-MIP depth 16
Nodes 431
Repair LPs 0
LP iterations 257262
144703 (strong br.)
16380 (separation)
13105 (heuristics)
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
Visualize the Inferred Network#
Now that the solution has been found, we can visualize it using the cn.methods.carnival.visualize_network function.
cn.methods.carnival.visualize_network(df_res)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/backend/execute.py:76, in run_check(cmd, input_lines, encoding, quiet, **kwargs)
75 kwargs['stdout'] = kwargs['stderr'] = subprocess.PIPE
---> 76 proc = _run_input_lines(cmd, input_lines, kwargs=kwargs)
77 else:
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/backend/execute.py:96, in _run_input_lines(cmd, input_lines, kwargs)
95 def _run_input_lines(cmd, input_lines, *, kwargs):
---> 96 popen = subprocess.Popen(cmd, stdin=subprocess.PIPE, **kwargs)
98 stdin_write = popen.stdin.write
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/subprocess.py:1039, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group)
1036 self.stderr = io.TextIOWrapper(self.stderr,
1037 encoding=encoding, errors=errors)
-> 1039 self._execute_child(args, executable, preexec_fn, close_fds,
1040 pass_fds, cwd, env,
1041 startupinfo, creationflags, shell,
1042 p2cread, p2cwrite,
1043 c2pread, c2pwrite,
1044 errread, errwrite,
1045 restore_signals,
1046 gid, gids, uid, umask,
1047 start_new_session, process_group)
1048 except:
1049 # Cleanup if the child failed starting.
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/subprocess.py:1991, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group)
1990 if err_filename is not None:
-> 1991 raise child_exception_type(errno_num, err_msg, err_filename)
1992 else:
FileNotFoundError: [Errno 2] No such file or directory: PosixPath('dot')
The above exception was the direct cause of the following exception:
ExecutableNotFound Traceback (most recent call last)
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/IPython/core/formatters.py:1036, in MimeBundleFormatter.__call__(self, obj, include, exclude)
1033 method = get_real_method(obj, self.print_method)
1035 if method is not None:
-> 1036 return method(include=include, exclude=exclude)
1037 return None
1038 else:
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/jupyter_integration.py:98, in JupyterIntegration._repr_mimebundle_(self, include, exclude, **_)
96 include = set(include) if include is not None else {self._jupyter_mimetype}
97 include -= set(exclude or [])
---> 98 return {mimetype: getattr(self, method_name)()
99 for mimetype, method_name in MIME_TYPES.items()
100 if mimetype in include}
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/jupyter_integration.py:112, in JupyterIntegration._repr_image_svg_xml(self)
110 def _repr_image_svg_xml(self) -> str:
111 """Return the rendered graph as SVG string."""
--> 112 return self.pipe(format='svg', encoding=SVG_ENCODING)
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/piping.py:104, in Pipe.pipe(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
55 def pipe(self,
56 format: typing.Optional[str] = None,
57 renderer: typing.Optional[str] = None,
(...) 61 engine: typing.Optional[str] = None,
62 encoding: typing.Optional[str] = None) -> typing.Union[bytes, str]:
63 """Return the source piped through the Graphviz layout command.
64
65 Args:
(...) 102 '<?xml version='
103 """
--> 104 return self._pipe_legacy(format,
105 renderer=renderer,
106 formatter=formatter,
107 neato_no_op=neato_no_op,
108 quiet=quiet,
109 engine=engine,
110 encoding=encoding)
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/_tools.py:185, in deprecate_positional_args.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
177 wanted = ', '.join(f'{name}={value!r}'
178 for name, value in deprecated.items())
179 warnings.warn(f'The signature of {func_name} will be reduced'
180 f' to {supported_number} positional arg{s_}{qualification}'
181 f' {list(supported)}: pass {wanted} as keyword arg{s_}',
182 stacklevel=stacklevel,
183 category=category)
--> 185 return func(*args, **kwargs)
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/piping.py:121, in Pipe._pipe_legacy(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
112 @_tools.deprecate_positional_args(supported_number=1, ignore_arg='self')
113 def _pipe_legacy(self,
114 format: typing.Optional[str] = None,
(...) 119 engine: typing.Optional[str] = None,
120 encoding: typing.Optional[str] = None) -> typing.Union[bytes, str]:
--> 121 return self._pipe_future(format,
122 renderer=renderer,
123 formatter=formatter,
124 neato_no_op=neato_no_op,
125 quiet=quiet,
126 engine=engine,
127 encoding=encoding)
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/piping.py:149, in Pipe._pipe_future(self, format, renderer, formatter, neato_no_op, quiet, engine, encoding)
146 if encoding is not None:
147 if codecs.lookup(encoding) is codecs.lookup(self.encoding):
148 # common case: both stdin and stdout need the same encoding
--> 149 return self._pipe_lines_string(*args, encoding=encoding, **kwargs)
150 try:
151 raw = self._pipe_lines(*args, input_encoding=self.encoding, **kwargs)
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/backend/piping.py:212, in pipe_lines_string(engine, format, input_lines, encoding, renderer, formatter, neato_no_op, quiet)
206 cmd = dot_command.command(engine, format,
207 renderer=renderer,
208 formatter=formatter,
209 neato_no_op=neato_no_op)
210 kwargs = {'input_lines': input_lines, 'encoding': encoding}
--> 212 proc = execute.run_check(cmd, capture_output=True, quiet=quiet, **kwargs)
213 return proc.stdout
File /data/ddimitrov/software/miniforge3/envs/liana313/lib/python3.13/site-packages/graphviz/backend/execute.py:81, in run_check(cmd, input_lines, encoding, quiet, **kwargs)
79 except OSError as e:
80 if e.errno == errno.ENOENT:
---> 81 raise ExecutableNotFound(cmd) from e
82 raise
84 if not quiet and proc.stderr:
ExecutableNotFound: failed to execute PosixPath('dot'), make sure the Graphviz executables are on your systems' PATH
<graphviz.graphs.Digraph at 0x7fb442172660>
We can see that the network above, largely captures a potential regulatory cascade with inhibitory (–|) and stimulatory (–>) interactions, related to JAK-STAT signalling. The network, in this case, starts from a receptor (triangle), coming from the top interactions, and ends with the deregulated TFs (square). The remainder of the nodes (circles) were inferred, taking their weights into account, and were not necessarily included in the input or output nodes.
In this example, we represent the directionality of signalling such that intracellular signalling is downstream of intercellular communication events. However, in biology cellular response is an admixture of both; thus such approaches are a simplification of biological reality.
Describe Results#
Let’s examine the result of the subnetwork search - it provides information about the predicted signs of nodes and edges.
df_res.head()
| source | source_type | source_weight | source_pred_val | target | target_type | target_weight | target_pred_val | edge_type | edge_pred_val | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | STAT1 | unmeasured | 0.0 | 1.0 | IRF1 | output | 11.398815 | 1.0 | 1 | 1.0 |
| 1 | NFKBIB | unmeasured | 0.0 | -1.0 | NFKB1 | output | 7.847975 | 1.0 | -1 | 1.0 |
| 2 | IKBKG | unmeasured | 0.0 | -1.0 | NFKBIB | unmeasured | 0.000000 | -1.0 | 1 | -1.0 |
| 3 | JAK1 | unmeasured | 0.0 | 1.0 | STAT2 | output | 10.497779 | 1.0 | 1 | 1.0 |
| 4 | GNAI2 | unmeasured | 0.0 | 1.0 | SRC | unmeasured | 0.000000 | 1.0 | 1 | 1.0 |
Nodes#
source: source nodes in the Protein-Protein Interaction (PPI) network. Suffixes such as “_s,” “_pert_c0,” and “_meas_c0” indicate specific experimental conditions or measurement types (they are there simply because of how the ILP problem is formulated, and can be ignored).
target: target nodes in the PPI network, with suffixes similar to the source node.
source_type (unmeasured, input):
input: start nodes (provided by the users, here receptors).
output: end nodes (provides by the user, here transcription factors).
unmeasured: Nodes that are neither input nor output - i.e. those that predicted by the algorithm.
source_weight and target_weight: Inputs to the causal net method, indicating the influence of “measured” nodes within the network. Only the sign is taken into account.
source_pred_val (1, 0, -1): Regulatory state of the node:
1: Upregulated
0: No differential expression
-1: Downregulated
target_pred_val (1, -1): Regulatory state of the target node:
1: Upregulated
-1: Downregulated
Edges (interaction)#
edge_type (1, -1, 0): Type of interaction from prior knowledge:
1: Activating interaction (e.g., A -> B)
-1: Inhibitory interaction
edge_pred_val (1, -1): Predicted effect of the interaction on the target node:
1: Upregulation
-1: Downregulation
Installing the Gurobi Solver: A Step-by-Step Guide#
While in this small example HIGHS (an internal scipy solver works), for larger networks we recommend using a solver such as Gurobi.
Gurobi is a powerful optimization solver used in various mathematical programming problems. Here’s how you can install it:
1. Download Gurobi for Your Operating System:#
Visit the Gurobi download page and select the version compatible with your OS.
2. Unzip and Update Path:#
After downloading, unzip the file. Locate the /bin folder inside the unzipped directory and add it to your system’s $PATH variable. This step is crucial as it allows your system to recognize and run Gurobi from anywhere.
3. Register for an Academic License:#
If you’re an academic user, you can obtain a free license. Register and request an academic license through the Gurobi portal. Follow the prompts to complete your registration.
4. Install Gurobi Python Interface:#
Open your command prompt or terminal and run:
pip install gurobipy
This command installs the necessary Python interface to interact with Gurobi.
5. Configure the Solver:#
In your code, ensure the solver parameter is set to gurobi to direct your program to use the Gurobi solver.
By following these steps, you should have Gurobi installed and ready to tackle complex optimization problems.