Benchmarking pipeline

Biological activities can be estimated by coupling a statistical method with a prior knowledge network. However, it is not clear which method or network performs better. When perturbation experiments are available, one can benchmark these methods by evaluating how good are they at recovering perturbed regulators. In our original decoupler publication, we observed that wsum, mlm, ulm and consensus where the top performer methods.

In this notebook we showcase how to use decoupler for the benchmarking of methods and prior knowledge networks. The data used here consist of single-gene perturbation experiments from the KnockTF2 database, an extensive dataset of literature curated perturbation experiments. If you use these data for your resarch, please cite the original publication of the resource:

Chenchen Feng, Chao Song, Yuejuan Liu, Fengcui Qian, Yu Gao, Ziyu Ning, Qiuyu Wang, Yong Jiang, Yanyu Li, Meng Li, Jiaxin Chen, Jian Zhang, Chunquan Li, KnockTF: a comprehensive human gene expression profile database with knockdown/knockout of transcription factors, Nucleic Acids Research, Volume 48, Issue D1, 08 January 2020, Pages D93–D100, https://doi.org/10.1093/nar/gkz881

Note

This tutorial assumes that you already know the basics of decoupler. Else, check out the Usage tutorial first.

Loading packages

First, we need to load the relevant packages.

import pandas as pd
import numpy as np

import decoupler as dc

Loading the data

We can download the data easily from Zenodo:

!wget 'https://zenodo.org/record/7035528/files/knockTF_expr.csv?download=1' -O knockTF_expr.csv
!wget 'https://zenodo.org/record/7035528/files/knockTF_meta.csv?download=1' -O knockTF_meta.csv

--2023-06-01 11:09:30--  https://zenodo.org/record/7035528/files/knockTF_expr.csv?download=1
Resolving zenodo.org (zenodo.org)... 188.185.124.72
Connecting to zenodo.org (zenodo.org)|188.185.124.72|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 146086808 (139M) [text/plain]
Saving to: ‘knockTF_expr.csv’

knockTF_expr.csv    100%[===================>] 139,32M  18,8MB/s    in 47s

2023-06-01 11:10:17 (2,97 MB/s) - ‘knockTF_expr.csv’ saved [146086808/146086808]

--2023-06-01 11:10:18--  https://zenodo.org/record/7035528/files/knockTF_meta.csv?download=1
Resolving zenodo.org (zenodo.org)... 188.185.124.72
Connecting to zenodo.org (zenodo.org)|188.185.124.72|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 144861 (141K) [text/plain]
Saving to: ‘knockTF_meta.csv’

knockTF_meta.csv    100%[===================>] 141,47K  --.-KB/s    in 0,03s

2023-06-01 11:10:18 (4,15 MB/s) - ‘knockTF_meta.csv’ saved [144861/144861]

We can then read it using pandas:

# Read data
mat = pd.read_csv('knockTF_expr.csv', index_col=0)
obs = pd.read_csv('knockTF_meta.csv', index_col=0)

mat consist of a matrix of logFCs between the perturbed and the control samples. Positive values mean that these genes were overexpressed after the perturbation, negative values mean the opposite.

mat

	A1BG	A1BG-AS1	A1CF	A2LD1	A2M	A2ML1	A4GALT	A4GNT	AA06	AAA1	...	ZWINT	ZXDA	ZXDB	ZXDC	ZYG11A	ZYG11B	ZYX	ZZEF1	ZZZ3	psiTPTE22
DataSet_01_001	0.000000	0.000000	0.000000	0.0000	0.000000	0.00000	0.00000	0.00000	0.0	0.0	...	0.000000	0.00000	0.000000	0.000000	0.00000	0.000000	-0.033220	-0.046900	-0.143830	0.00000
DataSet_01_002	-0.173750	0.000000	-0.537790	-0.4424	-1.596640	-0.00010	-0.08260	-3.64590	0.0	0.0	...	0.000000	0.39579	0.000000	0.086110	0.28574	-0.130210	-0.276640	-0.881580	0.812950	0.47896
DataSet_01_003	-0.216130	0.000000	-0.220270	-0.0008	0.000000	0.15580	-0.35802	-0.32025	0.0	0.0	...	0.000000	-0.82964	0.000000	0.064470	0.59323	0.313950	-0.232600	0.065510	-0.147140	0.00000
DataSet_01_004	-0.255680	0.111050	-0.285270	0.0000	-0.035860	-0.46970	0.18559	-0.25601	0.0	0.0	...	0.000000	-0.39888	0.000000	0.104440	-0.16434	0.191460	0.415610	0.393840	0.127900	0.00000
DataSet_01_005	0.478500	-0.375710	-0.847180	0.0000	3.354450	0.17104	-0.34852	-0.95517	0.0	0.0	...	0.000000	0.24849	0.000000	-0.286920	-0.01815	0.119410	0.077850	0.234740	0.228690	0.00000
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
DataSet_04_039	0.668721	0.000000	0.120893	0.0000	0.263564	0.00000	0.00000	0.00000	0.0	0.0	...	-0.324808	0.00000	-0.370491	0.092874	0.00000	0.014075	-0.215169	0.141374	-0.076054	0.00000
DataSet_04_040	1.816027	0.000000	0.139721	0.0000	0.459498	0.00000	0.00000	0.00000	0.0	0.0	...	-0.321726	0.00000	-0.144934	-0.238343	0.00000	-0.134995	0.577406	0.082028	-0.171424	0.00000
DataSet_04_041	-0.174250	0.000000	-0.084499	0.0000	0.112492	0.00000	0.00000	0.00000	0.0	0.0	...	0.405369	0.00000	0.272526	-0.125475	0.00000	0.019763	-0.098584	0.072675	0.199813	0.00000
DataSet_04_042	0.976681	0.000000	0.099845	0.0000	0.366973	0.00000	0.00000	0.00000	0.0	0.0	...	0.397698	0.00000	-0.398350	0.037134	0.00000	0.409704	-0.056978	0.308122	-0.058794	0.00000
DataSet_04_043	0.677944	0.584963	0.271215	0.0000	0.104090	0.00000	0.00000	0.00000	0.0	0.0	...	-0.357371	0.00000	-0.152622	-0.076498	0.00000	0.047949	-0.169766	0.178053	0.199003	0.00000

907 rows × 21985 columns

On the other hand, obs contains the meta-data information of each perturbation experiment. Here one can find which TF was perturbed (TF column), and many other useful information.

obs

	TF	Species	Knock.Method	Biosample.Name	Profile.ID	Platform	TF.Class	TF.Superclass	Tissue.Type	Biosample.Type	Data.Source	Pubmed.ID	logFC
DataSet_01_001	ESR1	Homo sapiens	siRNA	MCF7	GSE10061	GPL3921	Nuclear receptors with C4 zinc fingers	Zinc-coordinating DNA-binding domains	Mammary_gland	Cell line	GEO	18631401	-0.713850
DataSet_01_002	HNF1A	Homo sapiens	shRNA	HuH7	GSE103128	GPL18180	Homeo domain factors	Helix-turn-helix domains	Liver	Cell line	GEO	29466992	0.164280
DataSet_01_003	MLXIP	Homo sapiens	shRNA	HA1ER	GSE11242	GPL4133	Basic helix-loop-helix factors (bHLH)	Basic domains	Embryo_kidney	Stem cell	GEO	18458340	0.262150
DataSet_01_004	CREB1	Homo sapiens	shRNA	K562	GSE12056	GPL570	Basic leucine zipper factors (bZIP)	Basic domains	Haematopoietic_and_lymphoid_tissue	Cell line	GEO	18801183	-0.950180
DataSet_01_005	POU5F1	Homo sapiens	siRNA	GBS6	GSE12320	GPL570	Homeo domain factors	Helix-turn-helix domains	Bone_marrow	Cell line	GEO	20203285	0.000000
...	...	...	...	...	...	...	...	...	...	...	...	...	...
DataSet_04_039	SRPK2	Homo sapiens	CRISPR	HepG2	ENCSR929EIP	-	Others	ENCODE_TF	Liver	Cell line	ENCODE	22955616	-1.392551
DataSet_04_040	WRN	Homo sapiens	shRNA	HepG2	ENCSR093FHC	-	Others	ENCODE_TF	Liver	Cell line	ENCODE	22955616	-0.173964
DataSet_04_041	YBX1	Homo sapiens	CRISPR	HepG2	ENCSR548OTL	-	Cold-shock domain factors	beta-Barrel DNA-binding domains	Liver	Cell line	ENCODE	22955616	-2.025170
DataSet_04_042	ZC3H8	Homo sapiens	shRNA	HepG2	ENCSR184YDW	-	C3H zinc finger factors	Zinc-coordinating DNA-binding domains	Liver	Cell line	ENCODE	22955616	-0.027152
DataSet_04_043	ZNF574	Homo sapiens	CRISPR	HepG2	ENCSR821GQN	-	C2H2 zinc finger factors	Zinc-coordinating DNA-binding domains	Liver	Cell line	ENCODE	22955616	0.990298

907 rows × 13 columns

Filtering

Since some of the experiments collected inside KnockTF are of low quality, it is best to remove them. We can do this by filtering perturbation experiments where the knocked TF does not show a strong negative logFC (meaning that the knock method did not work).

# Filter out low quality experiments (positive logFCs in knockout experiments)
msk = obs['logFC'] < -1
mat = mat[msk]
obs = obs[msk]
mat.shape, obs.shape

((388, 21985), (388, 13))

Evaluation

Single network

As an example, we will evaluate the gene regulatory network CollecTRI, for more information you can visit this other vignette.

We can retrieve this network by running:

# Get collectri
collectri = dc.get_collectri(organism='human', split_complexes=False)
collectri

	source	target	weight
0	MYC	TERT	1
1	SPI1	BGLAP	1
2	SMAD3	JUN	1
3	SMAD4	JUN	1
4	STAT5A	IL2	1
...	...	...	...
43173	NFKB	hsa-miR-143-3p	1
43174	AP1	hsa-miR-206	1
43175	NFKB	hsa-miR-21-5p	1
43176	NFKB	hsa-miR-224-5p	1
43177	AP1	hsa-miR-144	1

43178 rows × 3 columns

We will run the benchmark pipeline using the top performer methods from decoupler. If we would like, we could use any other, and we can optionaly change their parametters too (for more information check the decouple function).

To run the benchmark pipeline, we need input molecular data (mat), its associated metadata (obs) with the name of the column where the perturbed regulators are (perturb), and which direction the pertrubations are (sign, positive or negative).

# Example on how to set up decouple arguments
decouple_kws={
    'methods' : ['wsum', 'ulm', 'mlm'],
    'consensus': True,
    'args' : {
        'wsum' : {'times': 100}
    }
}

# Run benchmark pipeline
df = dc.benchmark(mat, obs, collectri, perturb='TF', sign=-1, verbose=True, decouple_kws=decouple_kws)

Extracting inputs...
Formating net...
109 experiments without sources in net, they will be removed.
Running methods...
52 features of mat are empty, they will be removed.
Running wsum on mat with 279 samples and 21933 targets for 766 sources.
Infering activities on 1 batches.

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:17<00:00, 17.22s/it]

52 features of mat are empty, they will be removed.
Running ulm on mat with 279 samples and 21933 targets for 766 sources.
52 features of mat are empty, they will be removed.
Running mlm on mat with 279 samples and 21933 targets for 766 sources.

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:02<00:00,  2.33s/it]

Calculating metrics...
Computing metrics...
Done.

The result of running the pipeline is a long-format dataframe containing different performance metrics across methods:

• auroc: Area Under the Receiver Operating characteristic Curve (AUROC)

• auprc: Area Under the Precision-Recall Curve (AUPRC), which by default uses a calibrated version where 0.5 indicates a random classifier (adapted from this publication).

• mauroc: Monte-Carlo Area Under the Precision-Recall Curve (AUPRC)

• mauprc: Monte-Carlo Area Under the Receiver Operating characteristic Curve (AUROC)

The Monte-Carlo metrics perform random permutations where each one samples the same number of true positives and true negatives to have balanced sets, returning one value per iteration done (10k by default).

df

	groupby	group	source	method	metric	score	ci
0	None	None	None	wsum_estimate	auroc	0.709890	0.001305
1	None	None	None	wsum_estimate	auprc	0.735267	0.001305
2	None	None	None	wsum_estimate	mcauroc	0.718978	0.001305
3	None	None	None	wsum_estimate	mcauroc	0.730129	0.001305
4	None	None	None	wsum_estimate	mcauroc	0.715317	0.001305
...	...	...	...	...	...	...	...
12007	None	None	None	consensus_estimate	mcauprc	0.771847	0.001305
12008	None	None	None	consensus_estimate	mcauprc	0.791121	0.001305
12009	None	None	None	consensus_estimate	mcauprc	0.768232	0.001305
12010	None	None	None	consensus_estimate	mcauprc	0.782516	0.001305
12011	None	None	None	consensus_estimate	mcauprc	0.783199	0.001305

12012 rows × 7 columns

We can easily visualize the results in a scatterplot:

dc.plot_metrics_scatter(df, x='mcauroc', y='mcauprc')

It looks like mlm this time had worse performance than other methods like ulm at recovering perturbed regulators, this might be caused because there are many TFs that share too many targets. Whenever there are multiple sources (TFs) and they share a lot of targets we recommend to use ulm. On the other hand, when there are few sources and they are orthogonal (they share few targets), we recommend to use mlm. Since consensus remains one of the top performer methods we also recommend its use but at a higher computational cost.

We can also visualise the distirbutions of the obtained Monte-Carlo metrics, for example mcauroc:

dc.plot_metrics_boxplot(df, metric='mcauroc')

Multiple networks

This pipeline can evaluate multiple networks at the same time. Now we will evaluate CollecTRI and DoRothEA, a literature-based network that we developed in the past. We will also add a randomized version as control.

# Get dorothea
dorothea = dc.get_dorothea(organism='human')

# Randomize nets
ctri_rand = dc.shuffle_net(collectri, target='target', weight='weight').drop_duplicates(['source', 'target'])
doro_rand = dc.shuffle_net(dorothea, target='target', weight='weight').drop_duplicates(['source', 'target'])

We can run the same pipeline but now using a dictionary of networks (and of extra arguments if needed).

# Build dictionary of networks to test
nets = {
    'r_collectri': ctri_rand,
    'r_dorothea': doro_rand,
    'collectri': collectri,
    'dorothea': dorothea,
}

# Example extra arguments
decouple_kws = {
    'r_collectri': {
        'methods' : ['ulm', 'mlm'],
        'consensus': True
    },
    'r_dorothea': {
        'methods' : ['ulm', 'mlm'],
        'consensus': True
    },
    'collectri': {
        'methods' : ['ulm', 'mlm'],
        'consensus': True
    },
    'dorothea': {
        'methods' : ['ulm', 'mlm'],
        'consensus': True
    }
}

# Run benchmark pipeline
df = dc.benchmark(mat, obs, nets, perturb='TF', sign=-1, verbose=True, decouple_kws=decouple_kws)

Using r_collectri network...
Extracting inputs...
Formating net...
109 experiments without sources in net, they will be removed.
Running methods...
52 features of mat are empty, they will be removed.
Running ulm on mat with 279 samples and 21933 targets for 762 sources.
52 features of mat are empty, they will be removed.
Running mlm on mat with 279 samples and 21933 targets for 762 sources.

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:01<00:00,  1.62s/it]

Calculating metrics...
Computing metrics...
Done.
Using r_dorothea network...
Extracting inputs...
Formating net...
175 experiments without sources in net, they will be removed.
Running methods...
55 features of mat are empty, they will be removed.
Running ulm on mat with 213 samples and 21930 targets for 295 sources.
55 features of mat are empty, they will be removed.
Running mlm on mat with 213 samples and 21930 targets for 295 sources.

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00,  2.03it/s]

Calculating metrics...
Computing metrics...
Done.
Using collectri network...
Extracting inputs...
Formating net...
109 experiments without sources in net, they will be removed.
Running methods...
52 features of mat are empty, they will be removed.
Running ulm on mat with 279 samples and 21933 targets for 766 sources.
52 features of mat are empty, they will be removed.
Running mlm on mat with 279 samples and 21933 targets for 766 sources.

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:01<00:00,  1.49s/it]

Calculating metrics...
Computing metrics...
Done.
Using dorothea network...
Extracting inputs...
Formating net...
174 experiments without sources in net, they will be removed.
Running methods...
55 features of mat are empty, they will be removed.
Running ulm on mat with 214 samples and 21930 targets for 297 sources.
55 features of mat are empty, they will be removed.
Running mlm on mat with 214 samples and 21930 targets for 297 sources.

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00,  1.50it/s]

Calculating metrics...
Computing metrics...
Done.

Like before we can try to visualize the results in a scatterplot, this time grouping by the network used:

dc.plot_metrics_scatter(df, groupby='net', x='mcauroc', y='mcauprc')

Since this plot is too crowded, we can plot it again separating the methods by columns:

dc.plot_metrics_scatter_cols(df, col='method', figsize=(9, 5), groupby='net')

Or by selecting only one of the methods, in this case consensus since it is the most reliable method:

dc.plot_metrics_scatter(df[df['method'] == 'consensus_estimate'], groupby='net', x='mcauroc', y='mcauprc', show_text=False)

As expected, the random network has a performance close to 0.5 while the rest show higher predictive performance. As previously reported, the collectri network has better performance than dorothea.

If needed, we can also plot the distributions of the Monte-Carlo metrics, now grouped by network:

dc.plot_metrics_boxplot(df, metric='mcauroc', groupby='net')

By source

The benchmark pipeline also allows to test the performance of individual regulators using the argument by='source'. For simplicity let us just use the best performing network:

# Example on how to set up decouple arguments
decouple_kws={
    'methods' : ['ulm', 'mlm'],
    'consensus': True,
}

# Run benchmark pipeline
df = dc.benchmark(mat, obs, collectri, perturb='TF', sign=-1, by='source', verbose=True, decouple_kws=decouple_kws)

Extracting inputs...
Formating net...
109 experiments without sources in net, they will be removed.
Running methods...
52 features of mat are empty, they will be removed.
Running ulm on mat with 279 samples and 21933 targets for 766 sources.
52 features of mat are empty, they will be removed.
Running mlm on mat with 279 samples and 21933 targets for 766 sources.

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:01<00:00,  1.55s/it]

Calculating metrics...
Computing metrics...
Done.

We can plot the results with:

dc.plot_metrics_scatter_cols(df, col='source', figsize=(9, 5), groupby='method')

It looks like we get good predictions for HIF1A, RELA and SOX2 among others. On the other side, we worse predictions for STAT1.

By meta-data groups

The benchmarking pipeline can also be run by grouping the input experiments using the groupby argument. Multiple groups can be selected at the same time but for simplicity, we will just test which knockout method seems to perform the best inhibitions:

# Example on how to set up decouple arguments
decouple_kws={
    'methods' : ['ulm', 'mlm'],
    'consensus': True,
}

# Run benchmark pipeline
df = dc.benchmark(mat, obs, collectri, perturb='TF', sign=-1, groupby='Knock.Method', verbose=True, decouple_kws=decouple_kws)

Extracting inputs...
Formating net...
109 experiments without sources in net, they will be removed.
Running methods...
52 features of mat are empty, they will be removed.
Running ulm on mat with 279 samples and 21933 targets for 766 sources.
52 features of mat are empty, they will be removed.
Running mlm on mat with 279 samples and 21933 targets for 766 sources.

100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:01<00:00,  1.39s/it]

Calculating metrics...
Computing metrics for groupby Knock.Method...
Done.

Again, we can visualize the results in a collection of scatterplots:

dc.plot_metrics_scatter_cols(df, col='method', figsize=(9, 5), groupby='group')

From our limited data-set, we can see that apparently shRNA works better than other knockout protocols.