Multi-condition signalling network inference#
In this notebook we show how to infer signalling networks for a multicondition setting.
Here we will use mouse RNA-seq data of a multicondition study, where we will compare normal to spontaneous leukemia and Sleeping-Beauty (SB) mediated leukemia.
In the first part, we will show how to estimate Transcription Factor activities from gene expression data, following the Decoupler tutorial for functional analysis.
Then, we will infer 2 networks for each of the 2 conditions and analyse differences.
Author: Irene Rigato (Francesca Finotello’s group)
Reviewers: Sophia Müller-Dott, Pablo Rodriguez-Mier
# --- Saezlab tools ---
# https://decoupler-py.readthedocs.io/
import gzip
import os
import shutil
import tempfile
import urllib.request
import decoupler as dc
import numpy as np
# https://omnipathdb.org/
import omnipath as op
# Additional packages
import pandas as pd
# --- Additional libs ---
# Pydeseq for differential expression analysis
from pydeseq2.dds import DefaultInference, DeseqDataSet
from pydeseq2.ds import DeseqStats
# https://saezlab.github.io/
import corneto as cn
cn.info()
|
|
max_time = 300
seed = 0
# loading GEO GSE148679 dataset
url = "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE148679&format=file&file=GSE148679%5Fcounts%5Fgsea%5Fanalysis%2Etxt%2Egz"
adata = None
with tempfile.TemporaryDirectory() as tmpdirname:
# Path for the gzipped file in the temp folder
gz_file_path = os.path.join(tmpdirname, "counts.txt.gz")
# Download the file
with urllib.request.urlopen(url) as response:
with open(gz_file_path, "wb") as out_file:
shutil.copyfileobj(response, out_file)
# Decompress the file
decompressed_file_path = gz_file_path[:-3] # Removing '.gz' extension
with gzip.open(gz_file_path, "rb") as f_in:
with open(decompressed_file_path, "wb") as f_out:
shutil.copyfileobj(f_in, f_out)
adata = pd.read_csv(decompressed_file_path, index_col=0, sep="\t").T
adata.head()
| ID | Itm2a | Sergef | Fam109a | Dhx9 | Fam71e2 | Ssu72 | Olfr1018 | Eif2b2 | Mks1 | Hebp2 | ... | Olfr372 | Gosr1 | Ctsw | Ryk | Rhd | Pxmp4 | Gm25500 | 4930455C13Rik | Prss39 | Reg4 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Ebf1.1 | 265 | 332 | 306 | 27823 | 0 | 2969 | 0 | 2614 | 231 | 6 | ... | 0 | 2967 | 674 | 150 | 708 | 830 | 29 | 9 | 7 | 0 |
| Ebf1.2 | 175 | 471 | 335 | 35978 | 0 | 3108 | 0 | 2924 | 192 | 0 | ... | 0 | 2724 | 155 | 142 | 869 | 867 | 5 | 4 | 16 | 0 |
| Ebf1.3 | 266 | 382 | 325 | 30367 | 0 | 3207 | 0 | 2940 | 211 | 7 | ... | 0 | 2547 | 87 | 54 | 291 | 767 | 7 | 7 | 6 | 0 |
| Ebf1.4 | 212 | 388 | 331 | 30406 | 0 | 3168 | 0 | 2871 | 209 | 8 | ... | 0 | 2568 | 169 | 115 | 1244 | 829 | 5 | 1 | 4 | 0 |
| PE.1 | 321 | 312 | 290 | 22272 | 0 | 2631 | 0 | 2236 | 109 | 10 | ... | 0 | 2466 | 837 | 191 | 921 | 888 | 2 | 2 | 8 | 0 |
5 rows × 24420 columns
from anndata import AnnData
adata = AnnData(adata, dtype=np.float32)
adata.var_names_make_unique()
adata
AnnData object with n_obs × n_vars = 53 × 24420
# renaming conditions for clarity
adata.obs["condition"] = np.select(
[
adata.obs.index.str.contains("WT"),
adata.obs.index.str.contains("Pax5"),
adata.obs.index.str.contains("Ebf1"),
adata.obs.index.str.contains("Leuk"),
adata.obs.index.str.contains("SB"),
],
["Normal", "Normal", "Normal", "Leukemic_spontaneous", "Leukemic_SB"],
default="Normal",
)
# Visualize metadata
adata.obs["condition"].value_counts()
condition
Leukemic_SB 31
Normal 15
Leukemic_spontaneous 7
Name: count, dtype: int64
Splitting data for network inference and validation#
We will validate the network for the comparison Leukemic_SB vs Normal (enough samples)
adata1 = adata[adata.obs["condition"].isin(["Normal", "Leukemic_spontaneous"])].copy()
adata2 = adata[adata.obs["condition"].isin(["Normal", "Leukemic_SB"])].copy()
Differential expression analysis per condition#
Contrast 1: Leukemic_spontaneous vs Normal#
# Obtain genes that pass the thresholds
dc.pp.filter_by_expr(
adata1,
group="condition",
min_count=10,
min_total_count=15,
large_n=10,
min_prop=0.8,
)
adata1
AnnData object with n_obs × n_vars = 22 × 13585
obs: 'condition'
# Estimation of differential expression
inference1 = DefaultInference()
dds1 = DeseqDataSet(
adata=adata1,
design_factors="condition",
refit_cooks=True,
inference=inference1,
)
dds1.deseq2()
Using None as control genes, passed at DeseqDataSet initialization
stat_res_leu1 = DeseqStats(dds1, contrast=["condition", "Leukemic_spontaneous", "Normal"], inference=inference1)
stat_res_leu1.summary()
Log2 fold change & Wald test p-value: condition Leukemic_spontaneous vs Normal
baseMean log2FoldChange lfcSE stat pvalue \
ID
Itm2a 319.854773 1.040954 0.173014 6.016588 1.781313e-09
Sergef 518.491011 1.165972 0.102390 11.387579 4.822033e-30
Fam109a 353.074249 0.270614 0.099379 2.723045 6.468333e-03
Dhx9 28622.656109 -0.009588 0.058553 -0.163747 8.699304e-01
Ssu72 3007.691599 -0.047053 0.044517 -1.056967 2.905267e-01
... ... ... ... ... ...
Ctsw 204.692803 -1.015983 0.513348 -1.979131 4.780130e-02
Ryk 82.285327 -2.689443 0.325351 -8.266294 1.381879e-16
Rhd 519.107839 -4.933983 0.449260 -10.982454 4.641417e-28
Pxmp4 791.774976 -0.367077 0.068277 -5.376307 7.602927e-08
Gm25500 9.166514 0.879435 0.429192 2.049048 4.045745e-02
padj
ID
Itm2a 5.852272e-09
Sergef 9.577094e-29
Fam109a 1.008288e-02
Dhx9 8.900440e-01
Ssu72 3.376801e-01
... ...
Ctsw 6.506168e-02
Ryk 8.517614e-16
Rhd 7.596825e-27
Pxmp4 2.123474e-07
Gm25500 5.581542e-02
[13585 rows x 6 columns]
results_df1 = stat_res_leu1.results_df
results_df1.sort_values(by="padj", ascending=True, inplace=False).head()
| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ID | ||||||
| Polm | 2883.081884 | -7.255498 | 0.281539 | -25.770836 | 1.883068e-146 | 1.279074e-142 |
| Egfl6 | 958.998149 | -6.296339 | 0.244241 | -25.779215 | 1.516794e-146 | 1.279074e-142 |
| Kdm5b | 2202.514358 | -5.510307 | 0.220788 | -24.957441 | 1.772980e-137 | 8.028642e-134 |
| Milr1 | 1647.365126 | 1.244237 | 0.050208 | 24.781408 | 1.422471e-135 | 4.831069e-132 |
| Adgre5 | 10346.763067 | -3.102891 | 0.125360 | -24.751915 | 2.956557e-135 | 8.032965e-132 |
Contrast 2: Leukemic_SB vs Normal#
# Obtain genes that pass the thresholds
dc.pp.filter_by_expr(
adata2,
group="condition",
min_count=10,
min_total_count=15,
large_n=10,
min_prop=0.8,
)
adata2
AnnData object with n_obs × n_vars = 46 × 13404
obs: 'condition'
# Estimation of differential expression
inference2 = DefaultInference()
dds2 = DeseqDataSet(
adata=adata2,
design_factors="condition",
refit_cooks=True,
inference=inference2,
)
dds2.deseq2()
Using None as control genes, passed at DeseqDataSet initialization
stat_res_leu2 = DeseqStats(dds2, contrast=["condition", "Leukemic_SB", "Normal"], inference=inference2)
stat_res_leu2.summary()
Log2 fold change & Wald test p-value: condition Leukemic_SB vs Normal
baseMean log2FoldChange lfcSE stat pvalue \
ID
Itm2a 336.316323 0.710476 0.224228 3.168537 1.532084e-03
Sergef 692.078769 1.218302 0.070644 17.245745 1.204506e-66
Fam109a 382.268647 0.329745 0.097073 3.396873 6.816056e-04
Dhx9 27794.304719 -0.031645 0.064081 -0.493828 6.214275e-01
Ssu72 2994.457694 0.004410 0.061412 0.071816 9.427485e-01
... ... ... ... ... ...
Ctsw 151.303280 -1.150780 0.374313 -3.074381 2.109398e-03
Ryk 50.267450 -2.400835 0.309995 -7.744751 9.576951e-15
Rhd 259.540718 -4.740190 0.424503 -11.166439 5.951916e-29
Pxmp4 783.587183 -0.148554 0.095118 -1.561793 1.183368e-01
Gm25500 8.618646 0.386326 0.302098 1.278810 2.009639e-01
padj
ID
Itm2a 2.906320e-03
Sergef 2.409731e-64
Fam109a 1.365860e-03
Dhx9 6.748999e-01
Ssu72 9.529865e-01
... ...
Ctsw 3.916660e-03
Ryk 6.703365e-14
Rhd 1.323043e-27
Pxmp4 1.573129e-01
Gm25500 2.530503e-01
[13404 rows x 6 columns]
results_df2 = stat_res_leu2.results_df
results_df2.sort_values(by="padj", ascending=True, inplace=False).head()
| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ID | ||||||
| Mfsd2b | 2252.744179 | -9.412063 | 0.308793 | -30.480193 | 4.769906e-204 | 6.393582e-200 |
| Itga2b | 8434.243135 | -11.914356 | 0.401437 | -29.679285 | 1.421308e-193 | 9.525603e-190 |
| Adgrl4 | 148.349960 | -9.117488 | 0.322229 | -28.295029 | 3.978579e-176 | 1.777629e-172 |
| Gp1ba | 2396.209744 | -8.226071 | 0.323265 | -25.446856 | 7.648486e-143 | 2.563008e-139 |
| Ppbp | 17321.675345 | -15.613716 | 0.621843 | -25.108767 | 3.989417e-139 | 1.069483e-135 |
Prior knowledge with Decoupler and Omnipath#
# Retrieve CollecTRI gene regulatory network (through Omnipath)
collectri = dc.op.collectri(organism="mouse")
collectri.head()
| source | target | weight | resources | references | sign_decision | |
|---|---|---|---|---|---|---|
| 0 | Myc | Tert | 1.0 | DoRothEA-A;ExTRI;HTRI;NTNU.Curated;Pavlidis202... | 10022128;10491298;10606235;10637317;10723141;1... | PMID |
| 1 | Spi1 | Bglap3 | 1.0 | ExTRI | 10022617 | default activation |
| 2 | Spi1 | Bglap | 1.0 | ExTRI | 10022617 | default activation |
| 3 | Spi1 | Bglap2 | 1.0 | ExTRI | 10022617 | default activation |
| 4 | Smad3 | Jun | 1.0 | ExTRI;NTNU.Curated;TFactS;TRRUST | 10022869;12374795 | PMID |
TF activity inference per condition#
Contrast 1: Leukemic_spontaneous vs Normal#
mat1 = results_df1[["stat"]].T.rename(index={"stat": "Leukemic_spontaneous.vs.Normal"})
mat1
| ID | Itm2a | Sergef | Fam109a | Dhx9 | Ssu72 | Eif2b2 | Mks1 | Vps28 | Setd6 | Psma4 | ... | Cd37 | Rag2 | Itgb1bp2 | Sec23ip | Gosr1 | Ctsw | Ryk | Rhd | Pxmp4 | Gm25500 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Leukemic_spontaneous.vs.Normal | 6.016588 | 11.387579 | 2.723045 | -0.163747 | -1.056967 | 7.738055 | 3.636604 | 2.671039 | 6.008118 | 3.939383 | ... | 3.891274 | -2.875231 | -0.371605 | -11.321515 | 2.340763 | -1.979131 | -8.266294 | -10.982454 | -5.376307 | 2.049048 |
1 rows × 13585 columns
tf_acts1, tf_pvals1 = dc.mt.ulm(data=mat1, net=collectri, verbose=True)
tf_acts1
| Abl1 | Ahr | Aire | Apex1 | Ar | Arid1a | Arid1b | Arid3a | Arid3b | Arid4a | ... | Zfpm1 | Zfpm2 | Zglp1 | Zgpat | Zhx2 | Zic1 | Zic2 | Zkscan3 | Zkscan4 | Zkscan7 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Leukemic_spontaneous.vs.Normal | 1.452427 | -0.114799 | -2.245622 | 0.882775 | -0.297864 | -0.176459 | -0.632378 | -0.453538 | -1.821644 | 0.260495 | ... | -1.439184 | 0.110701 | -1.527737 | -0.94774 | 1.098365 | 0.530819 | -0.030262 | 1.943526 | 1.943526 | 1.024536 |
1 rows × 640 columns
dc.pl.barplot(
data=tf_acts1,
name="Leukemic_spontaneous.vs.Normal",
top=25,
figsize=(3, 6),
)
Contrast 2: Leukemic_SB vs Normal#
mat2 = results_df2[["stat"]].T.rename(index={"stat": "Leukemic_SB.vs.Normal"})
mat2
| ID | Itm2a | Sergef | Fam109a | Dhx9 | Ssu72 | Eif2b2 | Mks1 | Vps28 | Setd6 | Psma4 | ... | Cd37 | Rag2 | Itgb1bp2 | Sec23ip | Gosr1 | Ctsw | Ryk | Rhd | Pxmp4 | Gm25500 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Leukemic_SB.vs.Normal | 3.168537 | 17.245745 | 3.396873 | -0.493828 | 0.071816 | 6.812841 | 6.314824 | 3.452622 | 5.558162 | 5.259825 | ... | 1.184029 | -2.941684 | -2.355503 | -7.619795 | 2.099912 | -3.074381 | -7.744751 | -11.166439 | -1.561793 | 1.27881 |
1 rows × 13404 columns
tf_acts2, tf_pvals2 = dc.mt.ulm(data=mat2, net=collectri, verbose=True)
tf_acts2
| Abl1 | Ahr | Aire | Apex1 | Ar | Arid1a | Arid1b | Arid3a | Arid3b | Arid4a | ... | Zfpm1 | Zfpm2 | Zglp1 | Zgpat | Zhx2 | Zic1 | Zic2 | Zkscan3 | Zkscan4 | Zkscan7 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Leukemic_SB.vs.Normal | 1.198672 | 0.664946 | -0.932051 | 1.41717 | 0.441607 | 0.676726 | -1.148908 | -0.229291 | -1.489206 | 0.929682 | ... | -1.618189 | -0.165162 | -0.623841 | -1.382759 | 1.355504 | 0.659013 | -0.049629 | 1.390874 | 1.390874 | 1.316169 |
1 rows × 625 columns
dc.pl.barplot(
data=tf_acts2,
name="Leukemic_SB.vs.Normal",
top=25,
figsize=(3, 6),
)
Retrieving potential receptors per condition#
# We obtain ligand-receptor interactions from Omnipath, and we keep only the receptors
# This is our list of a prior potential receptors from which we will infer the network
unique_receptors = set(
op.interactions.LigRecExtra.get(organisms="mouse", genesymbols=True)["target_genesymbol"].values.tolist()
)
len(unique_receptors)
849
Contrast 1: Leukemic_spontaneous vs Normal#
df_de_receptors1 = results_df1.loc[results_df1.index.intersection(unique_receptors)]
df_de_receptors1 = df_de_receptors1.sort_values(by="stat", ascending=False)
# We will take the top 20 receptors that increased the expression after treatment
df_top_receptors1 = df_de_receptors1.head(30)
df_top_receptors1.head()
| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ID | ||||||
| Mrc1 | 314.324326 | 3.312083 | 0.179667 | 18.434602 | 6.932705e-76 | 2.943150e-73 |
| St14 | 1073.819725 | 1.904158 | 0.116095 | 16.401736 | 1.858495e-60 | 3.078982e-58 |
| Cd48 | 5533.470898 | 0.664655 | 0.054604 | 12.172306 | 4.365699e-34 | 1.212843e-32 |
| Ifngr1 | 6551.268902 | 1.418280 | 0.120961 | 11.725146 | 9.473804e-32 | 2.174014e-30 |
| Lrrc4c | 82.859784 | 10.019258 | 0.895652 | 11.186547 | 4.745541e-29 | 8.584311e-28 |
Contrast 2: Leukemic_SB vs Normal#
df_de_receptors2 = results_df2.loc[results_df2.index.intersection(unique_receptors)]
df_de_receptors2 = df_de_receptors2.sort_values(by="stat", ascending=False)
# We will take the top 20 receptors that increased the expression after treatment
df_top_receptors2 = df_de_receptors2.head(30)
df_top_receptors2.head()
| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ID | ||||||
| Lsr | 770.417182 | 4.360408 | 0.346097 | 12.598814 | 2.143445e-36 | 8.425435e-35 |
| Ptprs | 5164.424536 | 1.303309 | 0.109463 | 11.906393 | 1.096177e-32 | 3.257906e-31 |
| Nlgn2 | 1309.701665 | 1.950302 | 0.165263 | 11.801208 | 3.847470e-32 | 1.106684e-30 |
| Cd244 | 2529.473998 | 2.992844 | 0.255033 | 11.735143 | 8.418380e-32 | 2.345945e-30 |
| Lrrc4c | 267.125842 | 10.650499 | 0.916744 | 11.617743 | 3.348651e-31 | 8.801043e-30 |
Inferring intracellular signalling network with CORNETO#
cn.info()
|
|
|
from corneto.methods import CarnivalFlow
# CarnivalFlow.show_citations()
Setting prior knowledge graph#
pkn = op.interactions.OmniPath.get(organisms="mouse", databases=["SIGNOR"], genesymbols=True)
pkn = pkn[pkn.consensus_direction == True]
pkn.head()
| source | target | source_genesymbol | target_genesymbol | is_directed | is_stimulation | is_inhibition | consensus_direction | consensus_stimulation | consensus_inhibition | curation_effort | references | sources | n_sources | n_primary_sources | n_references | references_stripped | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | P0C605 | Q9QZC1 | Prkg1 | Trpc3 | True | False | True | True | False | True | 9 | HPRD:14983059;KEA:14983059;ProtMapper:14983059... | HPRD;HPRD_KEA;HPRD_MIMP;KEA;MIMP;PhosphoPoint;... | 15 | 8 | 2 | 14983059;16331690 |
| 1 | P0C605 | Q9WVC5 | Prkg1 | Trpc7 | True | True | False | True | True | False | 3 | SIGNOR:21402151;TRIP:21402151;iPTMnet:21402151 | SIGNOR;TRIP;iPTMnet | 3 | 3 | 1 | 21402151 |
| 2 | Q8K2C7 | Q9EPK8 | Os9 | Trpv4 | True | True | True | True | True | True | 3 | HPRD:17932042;SIGNOR:17932042;TRIP:17932042 | HPRD;SIGNOR;TRIP | 3 | 3 | 1 | 17932042 |
| 3 | P35821 | Q91WD2 | Ptpn1 | Trpv6 | True | False | True | True | False | True | 11 | DEPOD:15894168;DEPOD:17197020;HPRD:15894168;In... | DEPOD;HPRD;IntAct;Lit-BM-17;SIGNOR;SPIKE_LC;TRIP | 7 | 6 | 2 | 15894168;17197020 |
| 4 | P68040 | Q8CIR4 | Rack1 | Trpm6 | True | False | True | True | False | True | 2 | SIGNOR:18258429;TRIP:18258429 | SIGNOR;TRIP | 2 | 2 | 1 | 18258429 |
pkn["interaction"] = pkn["is_stimulation"].astype(int) - pkn["is_inhibition"].astype(int)
sel_pkn = pkn[["source_genesymbol", "interaction", "target_genesymbol"]]
sel_pkn.head()
| source_genesymbol | interaction | target_genesymbol | |
|---|---|---|---|
| 0 | Prkg1 | -1 | Trpc3 |
| 1 | Prkg1 | 1 | Trpc7 |
| 2 | Os9 | 0 | Trpv4 |
| 3 | Ptpn1 | -1 | Trpv6 |
| 4 | Rack1 | -1 | Trpm6 |
# We create the CORNETO graph by importing the edges and interaction
G = cn.Graph.from_sif_tuples([(r[0], r[1], r[2]) for _, r in sel_pkn.iterrows() if r[1] != 0])
G.shape # nodes, edges
(4304, 9505)
Identifying target TFs per condition#
max_pval = 0.01
Contrast 1: Leukemic_spontaneous vs Normal#
# As measurements, we take the estimated TFs, we will filter out TFs with p-val > 0.001
significant_tfs1 = (
tf_acts1[tf_pvals1 <= max_pval].T.dropna().sort_values(by="Leukemic_spontaneous.vs.Normal", ascending=False)
)
significant_tfs1.head()
| Leukemic_spontaneous.vs.Normal | |
|---|---|
| Myc | 5.927294 |
| Nr3c1 | -3.807611 |
| Klf1 | -3.827762 |
| Pax2 | -3.833372 |
| Pknox1 | -3.854072 |
# We keep only the ones in the PKN graph
measurements1 = significant_tfs1.loc[significant_tfs1.index.intersection(G.V)].to_dict()[
"Leukemic_spontaneous.vs.Normal"
]
measurements1
{'Myc': 5.9272938726510445,
'Nr3c1': -3.8076109957400055,
'Pax2': -3.833372272814152,
'Pknox1': -3.8540720496124936,
'Spi1': -3.870956617433766,
'Cebpb': -4.041315700038792,
'Srf': -4.189456236895589,
'Ets1': -4.531838015984738,
'Stat6': -4.762013996738547,
'Fli1': -5.337254315751466,
'Sp3': -5.344906237393298,
'Runx1': -5.966766274845339,
'Gata1': -6.145405587377315,
'Sp1': -6.180554655220554}
Contrast 2: Leukemic_SB vs Normal#
# As measurements, we take the estimated TFs, we will filter out TFs with p-val > 0.001
significant_tfs2 = tf_acts2[tf_pvals2 <= max_pval].T.dropna().sort_values(by="Leukemic_SB.vs.Normal", ascending=False)
significant_tfs2.head()
| Leukemic_SB.vs.Normal | |
|---|---|
| Myc | 5.456814 |
| Phox2a | -4.102671 |
| Runx1 | -5.481601 |
| Klf1 | -5.749658 |
| Pknox1 | -5.823616 |
# We keep only the ones in the PKN graph
measurements2 = significant_tfs2.loc[significant_tfs2.index.intersection(G.V)].to_dict()["Leukemic_SB.vs.Normal"]
measurements2
{'Myc': 5.4568137702067805,
'Phox2a': -4.102671408127372,
'Runx1': -5.481600720105362,
'Pknox1': -5.823615535106708,
'Fli1': -7.269774022313037,
'Gata1': -8.17857460262512}
def balance_dicts_by_top_abs(dict1, dict2, discretize=False):
# Decide which is smaller and which is larger
if len(dict1) <= len(dict2):
small, large = dict1, dict2
else:
small, large = dict2, dict1
n = len(small)
# Take the n keys from large with the largest abs(values)
top_large_items = sorted(large.items(), key=lambda kv: abs(kv[1]), reverse=True)[:n]
# Prepare output dicts
small_balanced = dict(small) # already of size n
large_balanced = dict(top_large_items)
if discretize:
small_balanced = {k: np.sign(v) for k, v in small_balanced.items()}
large_balanced = {k: np.sign(v) for k, v in large_balanced.items()}
return small_balanced, large_balanced
# Discretize not needed, but simplifies multi-cond analysis since both conditions have the same cost
d_measurements1, d_measurements2 = balance_dicts_by_top_abs(measurements1, measurements2, discretize=True)
d_measurements1
{'Myc': np.float64(1.0),
'Phox2a': np.float64(-1.0),
'Runx1': np.float64(-1.0),
'Pknox1': np.float64(-1.0),
'Fli1': np.float64(-1.0),
'Gata1': np.float64(-1.0)}
d_measurements2
{'Sp1': np.float64(-1.0),
'Gata1': np.float64(-1.0),
'Runx1': np.float64(-1.0),
'Myc': np.float64(1.0),
'Sp3': np.float64(-1.0),
'Fli1': np.float64(-1.0)}
Creating a CARNIVAL problem per condition#
Contrast 1: Leukemic_spontaneous vs Normal#
# We will infer the direction, so for the inputs, we use a value of 0 (=unknown direction)
inputs1 = {k: 0 for k in df_top_receptors1.index.intersection(G.V).values}
inputs1
{'Ifngr1': 0,
'Lrrc4c': 0,
'Tnfrsf11a': 0,
'Cdon': 0,
'Notch1': 0,
'Ptprs': 0,
'Nlgn2': 0,
'Flt3': 0,
'Ephb4': 0,
'Il2rg': 0,
'Marco': 0,
'Axl': 0,
'Tlr4': 0,
'Ifngr2': 0,
'Znrf3': 0}
# Create the dataset in standard format
carnival_data1 = dict()
for inp, v in inputs1.items():
carnival_data1[inp] = dict(value=v, role="input", mapping="vertex")
for out, v in d_measurements1.items():
carnival_data1[out] = dict(value=v, role="output", mapping="vertex")
data1 = cn.Data.from_cdict({"sample1": carnival_data1})
data1
Data(n_samples=1, n_feats=[21])
Contrast 2: Leukemic_SB vs Normal#
# We will infer the direction, so for the inputs, we use a value of 0 (=unknown direction)
inputs2 = {k: 0 for k in df_top_receptors2.index.intersection(G.V).values}
inputs2
{'Ptprs': 0,
'Nlgn2': 0,
'Lrrc4c': 0,
'Tnfrsf11a': 0,
'Ifngr1': 0,
'Marco': 0,
'Ifngr2': 0,
'Ctla4': 0,
'Znrf3': 0,
'Flt3': 0,
'Ephb4': 0,
'Il9r': 0,
'Icam1': 0,
'Axl': 0}
# Create the dataset in standard format
carnival_data2 = dict()
for inp, v in inputs2.items():
carnival_data2[inp] = dict(value=v, role="input", mapping="vertex")
for out, v in d_measurements2.items():
carnival_data2[out] = dict(value=v, role="output", mapping="vertex")
data2 = cn.Data.from_cdict({"sample2": carnival_data2})
data2
Data(n_samples=1, n_feats=[20])
Solving multi-condition CARNIVAL problem with CORNETO#
data = cn.Data.from_cdict({"sample1": carnival_data1, "sample2": carnival_data2})
data
Data(n_samples=2, n_feats=[21 20])
from corneto.utils import check_gurobi
check_gurobi()
Gurobipy successfully imported.
Gurobi environment started successfully.
Starting optimization of the test model...
Test optimization was successful.
Gurobi environment disposed.
Gurobi is correctly installed and working.
True
def subgraph(G, inp, out):
inp_s = set(G.V).intersection(inp)
out_s = set(G.V).intersection(out)
Gs = G.prune(inp_s, out_s)
tot_inp = set(Gs.V).intersection(inp)
tot_out = set(Gs.V).intersection(out)
print(f"Inputs: ({len(tot_inp)}/{len(inp)}), Outputs: ({len(tot_out)}/{len(out)})")
return Gs
Gs1 = subgraph(G, df_top_receptors1.index.tolist(), list(d_measurements1.keys()))
Gs2 = subgraph(G, df_top_receptors2.index.tolist(), list(d_measurements2.keys()))
Inputs: (7/30), Outputs: (5/6)
Inputs: (4/30), Outputs: (6/6)
c = CarnivalFlow(lambda_reg=0, indirect_rule_penalty=1)
P = c.build(G, data)
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
len(set(c.processed_graph.V).intersection(d_measurements1.keys()))
5
len(set(c.processed_graph.V).intersection(d_measurements2.keys()))
6
# How many TFs are in common between conditions?
s1 = set(c.processed_graph.V).intersection(d_measurements1.keys())
s2 = set(c.processed_graph.V).intersection(d_measurements2.keys())
for g in s1.intersection(s2):
print(g, d_measurements1[g], d_measurements2[g])
Fli1 -1.0 -1.0
Runx1 -1.0 -1.0
Gata1 -1.0 -1.0
Myc 1.0 1.0
Which is the best lambda to choose?#
To choose a robust value of \(\lambda\) we sample multiple solutions, as in the tutorial https://saezlab.github.io/corneto/dev/tutorials/network-sampler.html#vertex-based-perturbation.
We define 3 metrics to choose \(\lambda\):
\(sign\_agreement\_ratio = \frac{n.\:of\:nodes\:matching\:sign\:with\:stat\:value}{total\:n.\:of\:nodes}\)
\(de\_overlap\_ratio = \frac{n.\:of\:nodes\:also\:differentially\:expressed}{total\:n.\:of\:nodes}\)
\(overall\_agreement = \frac{sign\_agreement\_ratio\:+\:de\_overlap\_ratio}{2}\)
# retrieving alternative network solutions. This snippet will take a while to run
from corneto.methods.sampler import sample_alternative_solutions
lambda_val = [0, 0.01, 0.1, 0.2, 0.3, 0.5, 0.7, 0.9, 0.95]
other_optimal_v = dict() # collecting the vertex sampled solutions for each lambda
index = dict()
for i in lambda_val:
print("lambda: ", i)
c = CarnivalFlow(lambda_reg=i, indirect_rule_penalty=1)
P = c.build(G, data)
vertex_results = sample_alternative_solutions(
P,
"vertex_value",
percentage=0.03,
scale=0.03,
rel_opt_tol=0.05,
max_samples=30, # number of alternative solutions to sample
solver_kwargs=dict(solver="gurobi", max_seconds=max_time, mip_gap=0.01, seed=seed),
)
other_optimal_v["lambda" + str(i)] = vertex_results["vertex_value"]
index["lambda" + str(i)] = c.processed_graph.V
lambda: 0
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
Set parameter Username
Set parameter LicenseID to value 2775543
Academic license - for non-commercial use only - expires 2027-02-05
lambda: 0.01
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
lambda: 0.1
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
lambda: 0.2
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
lambda: 0.3
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
lambda: 0.5
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
lambda: 0.7
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
lambda: 0.9
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
lambda: 0.95
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
def compute_nodes_metrics(nodes, inputs, measurements, results_df, column):
"""Computes the metrics for the nodes in the network based on the provided inputs, measurements, and results dataframe.
Args:
nodes (pd.DataFrame): data frame containing the node scores per condition.
inputs (dict): Dictionary of input nodes with their values.
measurements (dict): Dictionary of output nodes with their values.
results_df (pd.DataFrame): DataFrame containing the results of differential expression analysis.
column (str): Column name to be used for filtering conditions.
Returns:
sign_agreement_ratio (float): Ratio of nodes with the same sign as the DE analysis.
de_overlap_ratio (float): Ratio of DE genes that overlap with the nodes.
overall_agreement (float): Overall agreement between DE genes and nodes.
"""
nodes = nodes.loc[(nodes[column] != 0.0)][column] # nodes with a score different from 0
de = results_df.loc[results_df["padj"] < 0.05].index.tolist() # de genes
in_out_nodes = set(measurements.keys()).union(
set(nodes.index[nodes.index.isin(inputs.keys())])
) # collecting input and output nodes which we don´t want to evaluate
middle_nodes = set(nodes.index).difference(in_out_nodes)
commmon_genes = list(
set(results_df.index).intersection(middle_nodes)
) # common genes between network and DE analysis genes
stat = results_df.loc[commmon_genes, "stat"]
stat_sign = np.sign(stat)
middle_nodes_sign = np.sign(nodes.loc[commmon_genes])
sign_agreement_ratio = 0
if len(middle_nodes_sign) > 0:
sign_agreement_ratio = (middle_nodes_sign == stat_sign).sum() / len(middle_nodes_sign)
de_overlap_ratio = 0
if len(middle_nodes) > 0:
de_overlap_ratio = len(set(de).intersection(middle_nodes)) / len(middle_nodes)
overall_agreement = (np.array(de_overlap_ratio) + np.array(sign_agreement_ratio)) / 2
return sign_agreement_ratio, de_overlap_ratio, overall_agreement
# Initialize metric containers
metrics_avg_1_per_lambda = []
metrics_avg_2_per_lambda = []
metrics_sd_1_per_lambda = []
metrics_sd_2_per_lambda = []
for k in other_optimal_v.keys():
metrics1_per_sample = []
metrics2_per_sample = []
for i in range(len(other_optimal_v[k])):
n = pd.DataFrame(
other_optimal_v[k][i],
index=index[k],
columns=["vertex_activity1", "vertex_activity2"],
) # retrieving node scores
sign_agreement1, de_overlap1, overall_agreement1 = compute_nodes_metrics(
n, inputs1, d_measurements1, results_df1, "vertex_activity1"
)
sign_agreement2, de_overlap2, overall_agreement2 = compute_nodes_metrics(
n, inputs2, d_measurements2, results_df2, "vertex_activity2"
)
metrics1_per_sample.append(
{
"sign_agreement": sign_agreement1,
"de_overlap": de_overlap1,
"overall": overall_agreement1,
}
)
metrics2_per_sample.append(
{
"sign_agreement": sign_agreement2,
"de_overlap": de_overlap2,
"overall": overall_agreement2,
}
)
metrics1_per_sample_df = pd.DataFrame(metrics1_per_sample)
metrics2_per_sample_df = pd.DataFrame(metrics2_per_sample)
metrics_avg_1_per_lambda.append(metrics1_per_sample_df.mean().to_dict())
metrics_avg_2_per_lambda.append(metrics2_per_sample_df.mean().to_dict())
metrics_sd_1_per_lambda.append(metrics1_per_sample_df.std().to_dict())
metrics_sd_2_per_lambda.append(metrics2_per_sample_df.std().to_dict())
# Convert lists of dictionaries to DataFrames
m_avg_df1 = pd.DataFrame(metrics_avg_1_per_lambda)
m_avg_df2 = pd.DataFrame(metrics_avg_2_per_lambda)
m_sd_df1 = pd.DataFrame(metrics_sd_1_per_lambda)
m_sd_df2 = pd.DataFrame(metrics_sd_2_per_lambda)
##plotting metrics
import matplotlib.pyplot as plt
metric_name = ["sign_agreement", "de_overlap", "overall"]
metric_label = ["Sign agreement ratio", "DE overlap ratio", "Overall agreement"]
fig, axes = plt.subplots(1, 3, figsize=(25, 5), sharex=True)
for i in range(len(metric_name)):
axes[i].errorbar(
lambda_val,
m_avg_df1[metric_name[i]],
yerr=m_sd_df1[metric_name[i]],
label="Leukemic_spontaneous",
marker="o",
)
axes[i].errorbar(
lambda_val,
m_avg_df2[metric_name[i]],
yerr=m_sd_df2[metric_name[i]],
label="Leukemic_SB",
marker="o",
)
axes[i].grid()
axes[i].set_xscale("linear") # adjust scale according to your results
axes[i].set_xlabel("Lambda")
axes[i].set_ylabel(metric_label[i])
axes[i].set_title(metric_label[i] + " vs Lambda")
if i == 0:
axes[i].legend(loc="upper left")
We choose the lambda value basing on the Overall agreement metric. In this case a reasonable value of \(\lambda\) is 2
c = CarnivalFlow(lambda_reg=0.5, indirect_rule_penalty=1)
P = c.build(G, data)
# Optimal solutions are found when Gap is close to 0
P.solve(solver="GUROBI", verbosity=1)
Unreachable vertices for sample: 12
Unreachable vertices for sample: 7
===============================================================================
CVXPY
v1.9.2
===============================================================================
-------------------------------------------------------------------------------
Compilation
-------------------------------------------------------------------------------
-------------------------------------------------------------------------------
Numerical solver
-------------------------------------------------------------------------------
Set parameter OutputFlag to value 1
Set parameter QCPDual to value 1
Gurobi Optimizer version 12.0.2 build v12.0.2rc0 (mac64[arm] - Darwin 25.5.0 25F84)
CPU model: Apple M4 Pro
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads
Non-default parameters:
QCPDual 1
Optimize a model with 45059 rows, 15944 columns and 159976 nonzeros
Model fingerprint: 0xad65e9b6
Variable types: 3844 continuous, 12100 integer (12100 binary)
Coefficient statistics:
Matrix range [5e-01, 7e+02]
Objective range [5e-01, 2e+00]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+03]
Found heuristic solution: objective 0.0000000
Presolve removed 21164 rows and 1239 columns
Presolve time: 0.36s
Presolved: 23895 rows, 14705 columns, 137721 nonzeros
Variable types: 3431 continuous, 11274 integer (11261 binary)
Deterministic concurrent LP optimizer: primal and dual simplex
Showing primal log only...
Concurrent spin time: 0.01s
Solved with dual simplex
Root relaxation: objective -1.434375e+00, 2089 iterations, 0.14 seconds (0.31 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 -1.43438 0 93 0.00000 -1.43438 - - 0s
0 0 -0.93977 0 113 0.00000 -0.93977 - - 0s
0 0 -0.93977 0 78 0.00000 -0.93977 - - 0s
0 0 -0.92708 0 82 0.00000 -0.92708 - - 0s
0 0 -0.92045 0 76 0.00000 -0.92045 - - 0s
H 0 0 -0.5000000 -0.89063 78.1% - 0s
0 0 -0.89063 0 76 -0.50000 -0.89063 78.1% - 0s
Cutting planes:
Gomory: 3
Implied bound: 16
Zero half: 2
RLT: 2
BQP: 1
Explored 1 nodes (5098 simplex iterations) in 0.88 seconds (1.64 work units)
Thread count was 12 (of 12 available processors)
Solution count 2: -0.5 0
Optimal solution found (tolerance 1.00e-04)
Best objective -5.000000000000e-01, best bound -5.000000000000e-01, gap 0.0000%
-------------------------------------------------------------------------------
Summary
-------------------------------------------------------------------------------
Problem(Minimize(Expression(AFFINE, UNKNOWN, (1,))), [Inequality(Constant(CONSTANT, ZERO, (2420,))), Inequality(Variable((2420,), _flow)), Inequality(Constant(CONSTANT, ZERO, (712, 2))), Inequality(Variable((712, 2), _dag_layer)), Equality(Expression(AFFINE, UNKNOWN, (712,)), Constant(CONSTANT, ZERO, ())), Inequality(Expression(AFFINE, NONNEGATIVE, (2420, 2))), Equality(Expression(AFFINE, UNKNOWN, (12,)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, NONNEGATIVE, (12,)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, NONNEGATIVE, (12,)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, UNKNOWN, (7,)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, NONNEGATIVE, (7,)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, NONNEGATIVE, (7,)), Constant(CONSTANT, ZERO, ())), Inequality(Expression(AFFINE, NONNEGATIVE, (712, 2))), Inequality(Expression(AFFINE, UNKNOWN, (2406,))), Inequality(Expression(AFFINE, UNKNOWN, (2406,))), Inequality(Constant(CONSTANT, NONNEGATIVE, (688,))), Inequality(Expression(AFFINE, UNKNOWN, (2406,))), Inequality(Expression(AFFINE, UNKNOWN, (2406,))), Inequality(Constant(CONSTANT, NONNEGATIVE, (688,))), Inequality(Expression(AFFINE, NONNEGATIVE, (2420, 2))), Inequality(Expression(AFFINE, NONNEGATIVE, (2413, 2))), Inequality(Expression(AFFINE, NONNEGATIVE, (2413, 2))), Equality(Expression(AFFINE, NONNEGATIVE, (3,)), Constant(CONSTANT, ZERO, ())), Equality(Expression(AFFINE, NONNEGATIVE, (3,)), Constant(CONSTANT, ZERO, ())), Inequality(Expression(AFFINE, NONNEGATIVE, (2420,))), Inequality(Variable((2420,), edge_has_signal_OR, boolean=True))])
# optimization metrics
for o in P.objectives:
print(o.name, o.value)
error_sample1_0 4.0
penalty_indirect_rules_0 [0.]
error_sample2_1 5.0
penalty_indirect_rules_1 [0.]
regularization_edge_has_signal_OR 3.0
# having a look to edge values per condition
pd.DataFrame(
P.expr.edge_value.value,
index=c.processed_graph.E,
columns=["edge_activity1", "edge_activity2"],
).head(5)
| edge_activity1 | edge_activity2 | ||
|---|---|---|---|
| (Mapk14) | (Mapkapk2) | 0.0 | 0.0 |
| (Map2k3) | (Dyrk1b) | 0.0 | 0.0 |
| (Akt1) | (Chuk) | 0.0 | 0.0 |
| (Cdkn1a) | (Cdk1) | 0.0 | 0.0 |
| (Mdc1) | (Rnf8) | 0.0 | 0.0 |
# having a look to node values per condition
pd.DataFrame(
P.expr.vertex_value.value,
index=c.processed_graph.V,
columns=["vertex_activity1", "vertex_activity2"],
).head(5)
| vertex_activity1 | vertex_activity2 | |
|---|---|---|
| Fus | 0.0 | 0.0 |
| Pstpip1 | 0.0 | 0.0 |
| Ppp1r1b | 0.0 | 0.0 |
| Arhgap33 | 0.0 | 0.0 |
| Lrrk2 | 0.0 | 0.0 |
Inferred network for contrast 1: Leukemic_spontaneous vs Normal#
sol_edges1 = np.flatnonzero(np.abs(P.expr.edge_value.value[:, 0]) > 0.5)
c.processed_graph.plot_values(
vertex_values=P.expr.vertex_value.value[:, 0],
edge_values=P.expr.edge_value.value[:, 0],
edge_indexes=sol_edges1,
)
Inferred network for contrast 2: Leukemic_SB vs Normal#
sol_edges2 = np.flatnonzero(np.abs(P.expr.edge_value.value[:, 1]) > 0.5)
c.processed_graph.plot_values(
vertex_values=P.expr.vertex_value.value[:, 1],
edge_values=P.expr.edge_value.value[:, 1],
edge_indexes=sol_edges2,
)