Bulk functional analysis

Bulk RNA-seq yields many molecular readouts that are hard to interpret by themselves. One way of summarizing this information is by inferring pathway and transcription factor activities from prior knowledge.

In this notebook we showcase how to use decoupler for transcription factor (TF) and pathway activity inference from a human data-set. The data consists of 6 samples of hepatic stellate cells (HSC) where three of them were activated by the cytokine Transforming growth factor (TGF-β), it is available at GEO here.

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, scanpy to handle RNA-seq data and decoupler to use statistical methods.

import scanpy as sc
import decoupler as dc

# Only needed for processing
import numpy as np
import pandas as pd
from anndata import AnnData

Loading the data

We can download the data easily from GEO:

!wget 'https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE151251&format=file&file=GSE151251%5FHSCs%5FCtrl%2Evs%2EHSCs%5FTGFb%2Ecounts%2Etsv%2Egz' -O counts.txt.gz
!gzip -d -f counts.txt.gz

--2023-06-01 11:29:28--  https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE151251&format=file&file=GSE151251%5FHSCs%5FCtrl%2Evs%2EHSCs%5FTGFb%2Ecounts%2Etsv%2Egz
Resolving www.ncbi.nlm.nih.gov (www.ncbi.nlm.nih.gov)... 130.14.29.110, 2607:f220:41e:4290::110
Connecting to www.ncbi.nlm.nih.gov (www.ncbi.nlm.nih.gov)|130.14.29.110|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1578642 (1,5M) [application/octet-stream]
Saving to: ‘counts.txt.gz’

counts.txt.gz       100%[===================>]   1,50M  2,14MB/s    in 0,7s

2023-06-01 11:29:29 (2,14 MB/s) - ‘counts.txt.gz’ saved [1578642/1578642]

We can then read it using pandas:

# Read raw data and process it
adata = pd.read_csv('counts.txt', index_col=2, sep='\t').iloc[:, 5:].T
adata

GeneName	DDX11L1	WASH7P	MIR6859-1	MIR1302-11	MIR1302-9	FAM138A	OR4G4P	OR4G11P	OR4F5	RP11-34P13.7	...	MT-ND4	MT-TH	MT-TS2	MT-TL2	MT-ND5	MT-ND6	MT-TE	MT-CYB	MT-TT	MT-TP
25_HSCs-Ctrl1	0	9	10	1	0	0	0	0	0	33	...	93192	342	476	493	54466	17184	1302	54099	258	475
26_HSCs-Ctrl2	0	12	14	0	0	0	0	0	0	66	...	114914	355	388	436	64698	21106	1492	62679	253	396
27_HSCs-Ctrl3	0	14	10	0	0	0	0	0	0	52	...	155365	377	438	480	85650	31860	2033	89559	282	448
31_HSCs-TGFb1	0	11	16	0	0	0	0	0	0	54	...	110866	373	441	481	60325	19496	1447	66283	172	341
32_HSCs-TGFb2	0	5	8	0	0	0	0	0	0	44	...	45488	239	331	343	27442	9054	624	27535	96	216
33_HSCs-TGFb3	0	12	5	0	0	0	0	0	0	32	...	70704	344	453	497	45443	13796	1077	43415	192	243

6 rows × 64253 columns

Note

In case your data does not contain gene symbols but rather gene ids like ENSMBL, you can use the function sc.queries.biomart_annotations to retrieve them. In this other vignette, ENSMBL ids are transformed into gene symbols.

Note

In case your data is not from human but rather a model organism, you can find their human orthologs using pypath. Check this GitHub issue where mouse genes are transformed into human.

The obtained data consist of raw read counts for six different samples (three controls, three treatments) for ~60k genes. Before continuing, we will transform our expression data into an AnnData object. It handles annotated data matrices efficiently in memory and on disk and is used in the scverse framework. You can read more about it here and here.

# Transform to AnnData object
adata = AnnData(adata, dtype=np.float32)
adata.var_names_make_unique()
adata

/home/badi/miniconda3/envs/dcp39/lib/python3.9/site-packages/anndata/_core/anndata.py:1830: UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`.
  utils.warn_names_duplicates("var")
/home/badi/miniconda3/envs/dcp39/lib/python3.9/site-packages/anndata/utils.py:111: UserWarning: Suffix used (-[0-9]+) to deduplicate index values may make index values difficult to interpret. There values with a similar suffixes in the index. Consider using a different delimiter by passing `join={delimiter}`Example key collisions generated by the make_index_unique algorithm: ['SNORD116-1', 'SNORD116-2', 'SNORD116-3', 'SNORD116-4', 'SNORD116-5']
  warnings.warn(

AnnData object with n_obs × n_vars = 6 × 64253

Inside an AnnData object, there is the .obs attribute where we can store the metadata of our samples. Here we will infer the metadata by processing the sample ids, but this could also be added from a separate dataframe:

# Process treatment information
adata.obs['condition'] = ['control' if '-Ctrl' in sample_id else 'treatment' for sample_id in adata.obs.index]

# Process sample information
adata.obs['sample_id'] = [sample_id.split('_')[0] for sample_id in adata.obs.index]

# Visualize metadata
adata.obs

	condition	sample_id
25_HSCs-Ctrl1	control	25
26_HSCs-Ctrl2	control	26
27_HSCs-Ctrl3	control	27
31_HSCs-TGFb1	treatment	31
32_HSCs-TGFb2	treatment	32
33_HSCs-TGFb3	treatment	33

Quality control

Before doing anything we need to ensure that our data passes some quality control thresholds. In transcriptomics it can happen that some genes were not properly profiled and thus need to be removed.

To filter genes, we will follow the strategy implemented in the function filterByExpr from edgeR. It keeps genes that have a minimum total number of reads across samples (min_total_count), and that have a minimum number of counts in a number of samples (min_count).

We can plot how many genes do we keep, you can play with the min_count and min_total_count to check how many genes would be kept when changed:

dc.plot_filter_by_expr(adata, group=None, min_count=10, min_total_count=15, large_n=1, min_prop=1)

Here we can observe the frequency of genes with different total sum of counts and number of samples. The dashed lines indicate the current thresholds, meaning that only the genes in the upper-right corner are going to be kept. Filtering parameters is completely arbitrary, but a good rule of thumb is to identify bimodal distributions and split them modifying the available thresholds. In this example, with the default values we would keep a good quantity of genes while filtering potential noisy genes.

Note

Changing the value of min_count will drastically change the distribution of “Number of samples”, not change its threshold. In case you want to lower or increase it, you need to play with the group, large_n and min_prop parameters.

Note

This thresholds can vary a lot between datasets, manual assessment of them needs to be considered. For example, it might be the case that many genes are not expressed in just one sample which they would get removed by the current setting. For this specific dataset it is fine.

Once we are content with the threshold parameters, we can perform the actual filtering:

# Obtain genes that pass the thresholds
genes = dc.filter_by_expr(adata, group=None, min_count=10, min_total_count=15, large_n=1, min_prop=1)

# Filter by these genes
adata = adata[:, genes].copy()
adata

AnnData object with n_obs × n_vars = 6 × 17575
    obs: 'condition', 'sample_id'

Differential expression analysis

In order to identify which are the genes that are changing the most between treatment and control we can perform differential expression analysis (DEA). For this example, we will perform a simple experimental design where we compare the gene expression of treated cells against controls. We will use the python implementation of the framework DESeq2, but we could have used any other one (limma or edgeR for example). For a better understanding how it works, check DESeq2’s documentation. Note that more complex experimental designs can be used by adding more factors to the design_factors argument.

# Import DESeq2
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

# Build DESeq2 object
dds = DeseqDataSet(
    adata=adata,
    design_factors='condition',
    refit_cooks=True,
    n_cpus=8,
)

# Compute LFCs
dds.deseq2()

Fitting size factors...
... done in 0.01 seconds.

Fitting dispersions...
... done in 39.97 seconds.

Fitting dispersion trend curve...
... done in 7.35 seconds.

Fitting MAP dispersions...
... done in 43.56 seconds.

Fitting LFCs...
... done in 3.16 seconds.

Refitting 0 outliers.

# Extract contrast between COVID-19 vs normal
stat_res = DeseqStats(dds, contrast=["condition", 'treatment', 'control'], n_cpus=8)

# Compute Wald test
stat_res.summary()

Running Wald tests...
... done in 1.53 seconds.

Log2 fold change & Wald test p-value: condition treatment vs control

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
GeneName
RP11-34P13.7	45.820190	0.067094	0.294257	0.228010	0.819639	0.871654
RP11-34P13.8	29.560568	-0.086227	0.249617	-0.345438	0.729765	0.801351
CICP27	106.257057	0.144798	0.177248	0.816925	0.413971	0.521994
FO538757.2	33.114796	-0.634443	0.286075	-2.217747	0.026572	0.050794
AP006222.2	67.154701	0.579098	0.270626	2.139848	0.032367	0.060587
...	...	...	...	...	...	...
MT-ND6	17962.507812	-0.445816	0.278378	-1.601475	0.109272	0.175000
MT-TE	1284.633179	-0.343227	0.287943	-1.191996	0.233263	0.330373
MT-CYB	55097.121094	-0.323882	0.302438	-1.070903	0.284213	0.387995
MT-TT	205.288528	-0.495090	0.332233	-1.490190	0.136174	0.210545
MT-TP	346.101318	-0.473558	0.200422	-2.362806	0.018137	0.036124

17575 rows × 6 columns

# Shrink LFCs
stat_res.lfc_shrink(coeff='condition_treatment_vs_control')

Fitting MAP LFCs...
... done in 6.70 seconds.

Shrunk Log2 fold change & Wald test p-value: condition treatment vs control

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
GeneName
RP11-34P13.7	45.820190	0.050407	0.271447	0.228010	0.819639	0.871654
RP11-34P13.8	29.560568	-0.068949	0.234608	-0.345438	0.729765	0.801351
CICP27	106.257057	0.129368	0.172303	0.816925	0.413971	0.521994
FO538757.2	33.114796	-0.523892	0.280505	-2.217747	0.026572	0.050794
AP006222.2	67.154701	0.483006	0.266258	2.139848	0.032367	0.060587
...	...	...	...	...	...	...
MT-ND6	17962.507812	-0.240886	0.260519	-1.601475	0.109272	0.175000
MT-TE	1284.633179	-0.265368	0.273098	-1.191996	0.233263	0.330373
MT-CYB	55097.121094	-0.523154	0.297741	-1.070903	0.284213	0.387995
MT-TT	205.288528	-0.366886	0.316293	-1.490190	0.136174	0.210545
MT-TP	346.101318	-0.424278	0.197644	-2.362806	0.018137	0.036124

17575 rows × 6 columns

# Extract results
results_df = stat_res.results_df
results_df

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
GeneName
RP11-34P13.7	45.820190	0.050407	0.271447	0.228010	0.819639	0.871654
RP11-34P13.8	29.560568	-0.068949	0.234608	-0.345438	0.729765	0.801351
CICP27	106.257057	0.129368	0.172303	0.816925	0.413971	0.521994
FO538757.2	33.114796	-0.523892	0.280505	-2.217747	0.026572	0.050794
AP006222.2	67.154701	0.483006	0.266258	2.139848	0.032367	0.060587
...	...	...	...	...	...	...
MT-ND6	17962.507812	-0.240886	0.260519	-1.601475	0.109272	0.175000
MT-TE	1284.633179	-0.265368	0.273098	-1.191996	0.233263	0.330373
MT-CYB	55097.121094	-0.523154	0.297741	-1.070903	0.284213	0.387995
MT-TT	205.288528	-0.366886	0.316293	-1.490190	0.136174	0.210545
MT-TP	346.101318	-0.424278	0.197644	-2.362806	0.018137	0.036124

17575 rows × 6 columns

We can plot the obtained results in a volcano plot:

dc.plot_volcano_df(results_df, x='log2FoldChange', y='padj', top=20)

/home/badi/miniconda3/envs/dcp39/lib/python3.9/site-packages/pandas/core/arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)

After performing DEA, we can use the obtained gene level statistics to perform enrichment analysis. Any statistic can be used, but we recommend using the t-values instead of logFCs since t-values incorporate the significance of change in their value. We will transform the obtained t-values stored in stats to a wide matrix so that it can be used by decoupler:

mat = results_df[['stat']].T.rename(index={'stat': 'treatment.vs.control'})
mat

GeneName	RP11-34P13.7	RP11-34P13.8	CICP27	FO538757.2	AP006222.2	RP4-669L17.10	MTND1P23	MTND2P28	AC114498.1	MIR6723	...	MT-ND4	MT-TH	MT-TS2	MT-TL2	MT-ND5	MT-ND6	MT-TE	MT-CYB	MT-TT	MT-TP
treatment.vs.control	0.22801	-0.345438	0.816925	-2.217747	2.139848	-0.422515	-1.373661	-0.53866	-2.845663	-3.862065	...	-1.472827	0.550156	0.704443	1.049607	-1.283784	-1.601475	-1.191996	-1.070903	-1.49019	-2.362806

1 rows × 17575 columns

Transcription factor activity inference

The first functional analysis we can perform is to infer transcription factor (TF) activities from our transcriptomics data. We will need a gene regulatory network (GRN) and a statistical method.

CollecTRI network

CollecTRI is a comprehensive resource containing a curated collection of TFs and their transcriptional targets compiled from 12 different resources. This collection provides an increased coverage of transcription factors and a superior performance in identifying perturbed TFs compared to our previous DoRothEA network and other literature based GRNs. Similar to DoRothEA, interactions are weighted by their mode of regulation (activation or inhibition).

For this example we will use the human version (mouse and rat are also available). We can use decoupler to retrieve it from omnipath. The argument split_complexes keeps complexes or splits them into subunits, by default we recommend to keep complexes together.

Note

In this tutorial we use the network CollecTRI, but we could use any other GRN coming from an inference method such as CellOracle, pySCENIC or SCENIC+.

# Retrieve CollecTRI gene regulatory network
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

Activity inference with Univariate Linear Model (ULM)

To infer TF enrichment scores we will run the Univariate Linear Model (ulm) method. For each sample in our dataset (mat) and each TF in our network (net), it fits a linear model that predicts the observed gene expression based solely on the TF’s TF-Gene interaction weights. Once fitted, the obtained t-value of the slope is the score. If it is positive, we interpret that the TF is active and if it is negative we interpret that it is inactive.

We can run ulm with a one-liner:

# Infer TF activities with ulm
tf_acts, tf_pvals = dc.run_ulm(mat=mat, net=collectri, verbose=True)
tf_acts

Running ulm on mat with 1 samples and 17575 targets for 629 sources.

	ABL1	AHR	AIRE	AP1	APEX1	AR	ARID1A	ARID3A	ARID3B	ARID4A	...	ZNF354C	ZNF362	ZNF382	ZNF384	ZNF395	ZNF436	ZNF699	ZNF76	ZNF804A	ZNF91
treatment.vs.control	-1.757384	-2.041384	-1.317279	-2.718626	-2.338566	0.178862	-6.440889	1.410645	2.23343	-0.704304	...	1.164067	-0.933561	2.687689	1.49982	-1.249087	1.407577	0.820753	0.485726	0.782225	0.86518

1 rows × 629 columns

Let us plot the obtained scores for the top active/inactive transcription factors:

dc.plot_barplot(tf_acts, 'treatment.vs.control', top=25, vertical=True)

SRF, MYOCD and SMAD3 seem to be the most activated in this treatment while IRF1, KLF11 AND ARID1A seem to be inactivated.

As before, we can manually inspect the downstream targets of each transcription factor:

dc.plot_targets(results_df, stat='stat', source_name='SRF', net=collectri, top=15)

SRF seems to be active in treated cells since their positive targets are up-regulated.

If needed, we can also look at this at the volcano plot:

# Extract logFCs and pvals
logFCs = results_df[['log2FoldChange']].T.rename(index={'log2FoldChange': 'treatment.vs.control'})
pvals = results_df[['padj']].T.rename(index={'padj': 'treatment.vs.control'})

# Plot
dc.plot_volcano(logFCs, pvals, 'treatment.vs.control', name='SRF', net=collectri, top=10, sign_thr=0.05, lFCs_thr=0.5)

/home/badi/miniconda3/envs/dcp39/lib/python3.9/site-packages/pandas/core/internals/blocks.py:351: RuntimeWarning: divide by zero encountered in log10
  result = func(self.values, **kwargs)

Pathway activity inference

Another analysis we can perform is to infer pathway activities from our transcriptomics data.

PROGENy model

PROGENy is a comprehensive resource containing a curated collection of pathways and their target genes, with weights for each interaction. For this example we will use the human weights (other organisms are available) and we will use the top 500 responsive genes ranked by p-value. Here is a brief description of each pathway:

• Androgen: involved in the growth and development of the male reproductive organs.

• EGFR: regulates growth, survival, migration, apoptosis, proliferation, and differentiation in mammalian cells

• Estrogen: promotes the growth and development of the female reproductive organs.

• Hypoxia: promotes angiogenesis and metabolic reprogramming when O2 levels are low.

• JAK-STAT: involved in immunity, cell division, cell death, and tumor formation.

• MAPK: integrates external signals and promotes cell growth and proliferation.

• NFkB: regulates immune response, cytokine production and cell survival.

• p53: regulates cell cycle, apoptosis, DNA repair and tumor suppression.

• PI3K: promotes growth and proliferation.

• TGFb: involved in development, homeostasis, and repair of most tissues.

• TNFa: mediates haematopoiesis, immune surveillance, tumour regression and protection from infection.

• Trail: induces apoptosis.

• VEGF: mediates angiogenesis, vascular permeability, and cell migration.

• WNT: regulates organ morphogenesis during development and tissue repair.

To access it we can use decoupler.

# Retrieve PROGENy model weights
progeny = dc.get_progeny(top=500)
progeny

	source	target	weight	p_value
0	Androgen	TMPRSS2	11.490631	0.000000e+00
1	Androgen	NKX3-1	10.622551	2.242078e-44
2	Androgen	MBOAT2	10.472733	4.624285e-44
3	Androgen	KLK2	10.176186	1.944414e-40
4	Androgen	SARG	11.386852	2.790209e-40
...	...	...	...	...
6995	p53	ZMYM4	-2.325752	1.522388e-06
6996	p53	CFDP1	-1.628168	1.526045e-06
6997	p53	VPS37D	2.309503	1.537098e-06
6998	p53	TEDC1	-2.274823	1.547037e-06
6999	p53	CCDC138	-3.205113	1.568160e-06

7000 rows × 4 columns

Activity inference with Multivariate Linear Model (MLM)

To infer pathway enrichment scores we will run the Multivariate Linear Model (mlm) method. For each sample in our dataset (adata), it fits a linear model that predicts the observed gene expression based on all pathways’ Pathway-Gene interactions weights. Once fitted, the obtained t-values of the slopes are the scores. If it is positive, we interpret that the pathway is active and if it is negative we interpret that it is inactive.

We can run mlm with a one-liner:

# Infer pathway activities with mlm
pathway_acts, pathway_pvals = dc.run_mlm(mat=mat, net=progeny, verbose=True)
pathway_acts

Running mlm on mat with 1 samples and 17575 targets for 14 sources.

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

	Androgen	EGFR	Estrogen	Hypoxia	JAK-STAT	MAPK	NFkB	PI3K	TGFb	TNFa	Trail	VEGF	WNT	p53
treatment.vs.control	1.92976	-1.739555	1.506251	6.360106	-11.535424	7.474121	0.744243	7.730807	37.363144	-5.540495	-2.149369	3.378371	-0.239079	-5.815882

Let us plot the obtained scores:

dc.plot_barplot(pathway_acts, 'treatment.vs.control', top=25, vertical=False)

As expected, after treating cells with the cytokine TGFb we see an increase of activity for this pathway.

On the other hand, it seems that this treatment has decreased the activity of other pathways like JAK-STAT or TNFa.

We can visualize the targets of TFGb in a scatter plot:

dc.plot_targets(results_df, stat='stat', source_name='TGFb', net=progeny, top=15)

The observed activation of TGFb is due to the fact that majority of its target genes with positive weights have positive t-values (1st quadrant), and the majority of the ones with negative weights have negative t-values (3d quadrant).

Functional enrichment of biological terms

Finally, we can also infer activities for general biological terms or processes.

MSigDB gene sets

The Molecular Signatures Database (MSigDB) is a resource containing a collection of gene sets annotated to different biological processes.

msigdb = dc.get_resource('MSigDB')
msigdb

	genesymbol	collection	geneset
0	MAFF	chemical_and_genetic_perturbations	BOYAULT_LIVER_CANCER_SUBCLASS_G56_DN
1	MAFF	chemical_and_genetic_perturbations	ELVIDGE_HYPOXIA_UP
2	MAFF	chemical_and_genetic_perturbations	NUYTTEN_NIPP1_TARGETS_DN
3	MAFF	immunesigdb	GSE17721_POLYIC_VS_GARDIQUIMOD_4H_BMDC_DN
4	MAFF	chemical_and_genetic_perturbations	SCHAEFFER_PROSTATE_DEVELOPMENT_12HR_UP
...	...	...	...
3838543	PRAMEF22	go_biological_process	GOBP_POSITIVE_REGULATION_OF_CELL_POPULATION_PR...
3838544	PRAMEF22	go_biological_process	GOBP_APOPTOTIC_PROCESS
3838545	PRAMEF22	go_biological_process	GOBP_REGULATION_OF_CELL_DEATH
3838546	PRAMEF22	go_biological_process	GOBP_NEGATIVE_REGULATION_OF_DEVELOPMENTAL_PROCESS
3838547	PRAMEF22	go_biological_process	GOBP_NEGATIVE_REGULATION_OF_CELL_DEATH

3838548 rows × 3 columns

As an example, we will use the hallmark gene sets, but we could have used any other.

Note

To see what other collections are available in MSigDB, type: msigdb['collection'].unique().

We can filter by for hallmark:

# Filter by hallmark
msigdb = msigdb[msigdb['collection']=='hallmark']

# Remove duplicated entries
msigdb = msigdb[~msigdb.duplicated(['geneset', 'genesymbol'])]

# Rename
msigdb.loc[:, 'geneset'] = [name.split('HALLMARK_')[1] for name in msigdb['geneset']]

msigdb

	genesymbol	collection	geneset
233	MAFF	hallmark	IL2_STAT5_SIGNALING
250	MAFF	hallmark	COAGULATION
270	MAFF	hallmark	HYPOXIA
373	MAFF	hallmark	TNFA_SIGNALING_VIA_NFKB
377	MAFF	hallmark	COMPLEMENT
...	...	...	...
1449668	STXBP1	hallmark	PANCREAS_BETA_CELLS
1450315	ELP4	hallmark	PANCREAS_BETA_CELLS
1450526	GCG	hallmark	PANCREAS_BETA_CELLS
1450731	PCSK2	hallmark	PANCREAS_BETA_CELLS
1450916	PAX6	hallmark	PANCREAS_BETA_CELLS

7318 rows × 3 columns

Enrichment with Over Representation Analysis (ORA)

To infer functional enrichment scores we will run the Over Representation Analysis (ora) method. As input data it accepts an expression matrix (decoupler.run_ora) or the results of differential expression analysis (decoupler.run_ora_df). For the former, by default the top 5% of expressed genes by sample are selected as the set of interest (S*), and for the latter a user-defined significance filtering can be used. Once we have S*, it builds a contingency table using set operations for each set stored in the gene set resource being used (net). Using the contingency table, ora performs a one-sided Fisher exact test to test for significance of overlap between sets. The final score is obtained by log-transforming the obtained p-values, meaning that higher values are more significant.

We can run ora with a simple one-liner:

# Infer enrichment with ora using significant deg
top_genes = results_df[results_df['padj'] < 0.05].reset_index()
top_genes['group'] = 'treatment.vs.control'

# Run ora
enr_pvals = dc.get_ora_df(
    df=top_genes,
    net=msigdb,
    groupby='group',
    features='GeneName',
    source='geneset',
    target='genesymbol'
)

enr_pvals

	ADIPOGENESIS	ALLOGRAFT_REJECTION	ANDROGEN_RESPONSE	ANGIOGENESIS	APICAL_JUNCTION	APICAL_SURFACE	APOPTOSIS	BILE_ACID_METABOLISM	CHOLESTEROL_HOMEOSTASIS	COAGULATION	...	PROTEIN_SECRETION	REACTIVE_OXYGEN_SPECIES_PATHWAY	SPERMATOGENESIS	TGF_BETA_SIGNALING	TNFA_SIGNALING_VIA_NFKB	UNFOLDED_PROTEIN_RESPONSE	UV_RESPONSE_DN	UV_RESPONSE_UP	WNT_BETA_CATENIN_SIGNALING	XENOBIOTIC_METABOLISM
treatment.vs.control	7.529597e-41	3.818284e-24	2.010532e-21	1.677850e-07	9.410098e-37	7.686361e-08	2.064548e-36	1.054623e-18	2.105839e-20	1.744101e-24	...	2.010532e-21	6.268350e-13	1.153058e-16	6.008191e-14	0.0	7.966153e-25	4.288777e-37	2.835032e-30	3.249205e-10	9.936055e-36

1 rows × 50 columns

We can then transform these p-values to their -log10 (so that the higher the value, the more significant the p-value). We will also set 0s to a minimum p-value so that we do not get infinites.

# Set 0s to min p-value
enr_pvals.values[enr_pvals.values == 0] = np.min(enr_pvals.values[enr_pvals.values != 0])

# Log-transform
enr_pvals = -np.log10(enr_pvals)

Then we can visualize the most enriched terms:

dc.plot_barplot(enr_pvals, 'treatment.vs.control', top=25, vertical=True)

TNFa and interferons response (JAK-STAT) processes seem to be enriched. We previously observed a similar result with the PROGENy pathways, were they were significantly downregulated. Therefore, one of the limitations of using a prior knowledge resource without weights is that it doesn’t provide direction.

Enrichment of ligand-receptor interactions

Recently, study of interactions between ligands and receptors have gained significant traction, notably pushed by the democratization of single cell sequencing technologies. While most methods (such as the one described in LIANA) are developed for single cell datasets, they rely on a relatively simple assumption of co-expression or co-regulation of two (or more in the context of complexes) genes acting as ligand and receptors to propose hypothetical ligand-receptor interaction events. This concept can seamlessly be applied to a bulk RNA dataset, where the assumption is that sender and receiver cells are convoluted in a single dataset, but the observation of significant co-regulation of ligand and receptors should still correspond to hypothetical ligand-receptor interaction events.

At the core, most current ligand-receptor interaction methods rely on averaging the measurements obtained for ligand and receptor and standardizing them against a background distribution. Thus, an enrichment method based on a weighted mean can emulate this, where the sets are simply the members of a ligand receptor pair (or more in the context of complexes).

Thus, we can extract ligand-receptor interaction ressources from the LIANA package (available both in R and python), and use it as a prior knowledge network with decoupler to find the most significant pairs of ligand-receptors in a given bulk dataset.

While more work is required to fully understand the functional relevance of the highlighted ligand-receptor interactions with such an approach, this represents a very straightforward and intuitive approach to embed a bulk RNA dataset with ligand-receptor interaction prior knowledge.

First, we extract ligand-receptor interactions from liana, and decomplexify them to format them into an appropriate decoupleR input.

import liana as ln

liana_lr = ln.resource.select_resource()
liana_lr = ln.resource.explode_complexes(liana_lr)

# Create two new DataFrames, each containing one of the pairs of columns to be concatenated
df1 = liana_lr[['interaction', 'ligand']]
df2 = liana_lr[['interaction', 'receptor']]

# Rename the columns in each new DataFrame
df1.columns = ['interaction', 'genes']
df2.columns = ['interaction', 'genes']

# Concatenate the two new DataFrames
liana_lr = pd.concat([df1, df2], axis=0)
liana_lr['weight'] = 1

# Find duplicated rows
duplicates = liana_lr.duplicated()

# Remove duplicated rows
liana_lr = liana_lr[~duplicates]

liana_lr

	interaction	genes	weight
0	LGALS9|PTPRC	LGALS9	1
1	LGALS9|MET	LGALS9	1
2	LGALS9|CD44	LGALS9	1
3	LGALS9|LRP1	LGALS9	1
4	LGALS9|CD47	LGALS9	1
...	...	...	...
5849	BMP2|ACTR2	ACTR2	1
5850	BMP15|ACTR2	ACTR2	1
5851	CSF1|CSF3R	CSF3R	1
5852	IL36G|IFNAR1	IFNAR1	1
5853	IL36G|IFNAR2	IFNAR2	1

10532 rows × 3 columns

Then we can use ulm to find the significant co-regulated pairs of ligand and receptors.

# Infer lr activities with ulm
lr_score, lr_pvalue = dc.run_ulm(
    mat=mat,
    net=liana_lr,
    source='interaction',
    target='genes',
    min_n=2,
    verbose=True
)

Running ulm on mat with 1 samples and 17575 targets for 1854 sources.

dc.plot_barplot(lr_score, 'treatment.vs.control', top=25, vertical=True)

Interactions between colagens and the ITG_ITG complexes seems to be quite enrichned. That is especially relevant since the EMT pathway was significantly enriched in the functional pathway analysis.