Transcription factor activity inference

scRNA-seq yield many molecular readouts that are hard to interpret by themselves. One way of summarizing this information is by infering transcription factor (TF) activities from prior knowledge.

In this notebook we showcase how to use decoupler for TF activity inference with the 3k PBMCs 10X data-set. The data consists of 3k PBMCs from a Healthy Donor and is freely available from 10x Genomics here from this webpage

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 scRNA-seq data and decoupler to use statistical methods.

[1]:
import scanpy as sc
import decoupler as dc

# Only needed for visualization:
import matplotlib.pyplot as plt
import seaborn as sns

Loading the data

We can download the data easily using scanpy:

[2]:
adata = sc.datasets.pbmc3k_processed()
adata
[2]:
AnnData object with n_obs × n_vars = 2638 × 1838
    obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain'
    var: 'n_cells'
    uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'

We can visualize the different cell types in it:

[3]:
sc.pl.umap(adata, color='louvain')
../_images/notebooks_dorothea_7_0.png

DoRothEA network

DoRothEA is a comprehensive resource containing a curated collection of TFs and their transcriptional targets. Since these regulons were gathered from different types of evidence, interactions in DoRothEA are classified in different confidence levels, ranging from A (highest confidence) to D (lowest confidence). Moreover, each interaction is weighted by its confidence level and the sign of its mode of regulation (activation or inhibition).

For this example we will use the human version (mouse is also available) and we will use the confidence levels ABC. To access it we can use decoupler.

[4]:
net = dc.get_dorothea(organism='human', levels=['A','B','C'])
net
[4]:
source confidence target weight
0 ETS1 A IL12B 1.000000
1 RELA A IL6 1.000000
2 MITF A BCL2A1 -1.000000
3 E2F1 A TRERF1 1.000000
4 MITF A BCL2 1.000000
... ... ... ... ...
32272 IKZF1 C PTK2B 0.333333
32273 IKZF1 C PRKCB 0.333333
32274 IKZF1 C PREX1 0.333333
32275 IRF4 C SLAMF7 0.333333
32276 ZNF83 C ZNF331 0.333333

32277 rows × 4 columns

Activity inference with Multivariate Linear Model

To infer activities we will run the Multivariate Linear Model method (mlm). It models the observed gene expression by using a regulatory adjacency matrix (target genes x TFs) as covariates of a linear model. The values of this matrix are the associated interaction weights. The obtained t-values of the fitted model are the activity scores.

To run decoupler methods, we need an input matrix (mat), an input prior knowledge network/resource (net), and the name of the columns of net that we want to use.

[5]:
dc.run_mlm(mat=adata, net=net, source='source', target='target', weight='weight', verbose=True)
58 features of mat are empty in 2635 samples, they will be ignored.
Running mlm on mat with 2638 samples and 13656 targets for 281 sources.
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1/1 [00:02<00:00,  2.11s/it]

The obtained scores (t-values)(mlm_estimate) and p-values (mlm_pvals) are stored in the .obsm key:

[6]:
adata.obsm['mlm_estimate']
[6]:
AHR AR ARID2 ARID3A ARNT ARNTL ASCL1 ATF1 ATF2 ATF3 ... ZNF217 ZNF24 ZNF263 ZNF274 ZNF384 ZNF584 ZNF592 ZNF639 ZNF644 ZNF740
AAACATACAACCAC-1 1.378150 2.264298 -0.873331 -0.393862 0.138301 -0.672725 -0.517485 -0.248989 3.503303 2.835994 ... 0.184589 -0.368378 -2.079983 -0.819110 -1.070388 -0.513945 -0.168149 1.016380 -0.426697 0.321187
AAACATTGAGCTAC-1 -0.771588 1.791014 -0.843765 -2.423640 1.121624 -1.710007 -0.098785 -0.437353 1.153923 0.325599 ... 0.975830 -1.015501 -1.298656 -0.695448 -0.717212 -0.854894 -0.634031 -0.014624 0.189120 0.082494
AAACATTGATCAGC-1 1.692741 2.992088 -0.532625 -1.705129 0.845746 -1.167037 -0.788174 1.828790 3.575234 2.071524 ... 0.001384 -0.424550 -2.065619 -0.535101 -0.751276 0.057936 0.295797 1.095647 1.089456 1.264060
AAACCGTGCTTCCG-1 0.950014 2.362606 -0.064996 -2.908076 0.925004 -1.022275 -0.007307 0.157358 2.789861 0.088665 ... 0.729626 -0.744912 -1.341288 -1.054425 -1.698280 -0.830635 -0.641135 -0.139098 -0.632303 -0.595420
AAACCGTGTATGCG-1 1.064510 1.690478 -0.638209 -1.559541 0.021607 -1.932919 -0.375869 -1.908948 -2.098662 2.609694 ... 1.399450 -0.787965 -1.913958 0.275981 -0.631440 -0.596082 0.775437 0.925863 -0.350852 0.218543
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
TTTCGAACTCTCAT-1 0.318669 3.333388 -0.411639 -1.804543 -0.024563 -1.427912 -0.363295 -0.074911 1.548207 -1.499464 ... 0.924315 0.169201 -1.506775 -1.248257 -0.712768 -0.823545 -0.452053 -0.257264 0.458856 -0.125243
TTTCTACTGAGGCA-1 0.030769 3.491299 -0.003791 -1.533298 1.548915 -1.666713 -0.371389 0.731881 2.436971 -0.304424 ... 0.801170 -0.535308 -1.734625 -0.584285 -1.229417 -1.455677 0.495567 1.191070 -0.316076 -0.160704
TTTCTACTTCCTCG-1 0.832451 1.216455 -0.077293 -2.480380 1.021595 -1.059853 -0.544001 0.110235 4.672344 0.357117 ... 0.063307 -0.711507 -0.494090 -0.123923 0.307083 -0.587551 -0.336608 -0.489183 -0.105758 0.136212
TTTGCATGAGAGGC-1 1.639109 -1.020677 -0.213234 -1.343528 -1.097390 -1.058503 -0.473822 -1.477212 3.461018 -0.166352 ... 0.090235 -0.554029 -1.828934 -0.717332 -1.356414 0.781463 -0.567849 1.205468 -0.513302 -0.485558
TTTGCATGCCTCAC-1 0.603507 1.357398 -0.001079 -0.808254 0.574953 -0.918306 0.057768 -0.138383 3.882593 0.035894 ... 0.123599 -0.818726 -1.486402 -0.168509 -0.349769 -0.609659 -0.524839 0.367849 -0.639759 -1.074101

2638 rows × 281 columns

Note: Each run of run_mlm overwrites what is inside of mlm_estimate and mlm_pvals. if you want to run mlm with other resources and still keep the activities inside the same AnnData object, you can store the results in any other key in .obsm with different names, for example:

[7]:
adata.obsm['dorothea_mlm_estimate'] = adata.obsm['mlm_estimate'].copy()
adata.obsm['dorothea_mlm_pvals'] = adata.obsm['mlm_pvals'].copy()
adata
[7]:
AnnData object with n_obs × n_vars = 2638 × 1838
    obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain'
    var: 'n_cells'
    uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr', 'mlm_estimate', 'mlm_pvals', 'dorothea_mlm_estimate', 'dorothea_mlm_pvals'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'

Visualization

To visualize the obtained scores, we can re-use many of scanpy’s plotting functions. First though, we need to extract them from the adata object.

[8]:
acts = dc.get_acts(adata, obsm_key='mlm_estimate')
acts
[8]:
AnnData object with n_obs × n_vars = 2638 × 281
    obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain'
    uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr', 'mlm_estimate', 'mlm_pvals', 'dorothea_mlm_estimate', 'dorothea_mlm_pvals'

dc.get_acts returns a new AnnData object which holds the obtained activities in its .X attribute, allowing us to re-use many scanpy functions, for example:

[9]:
sc.pl.umap(acts, color=['PAX5', 'louvain'], cmap='coolwarm', vcenter=0)
../_images/notebooks_dorothea_19_0.png

Here we observe the activity infered for PAX5 across cells, which it is particulary active in B cells. Interestingly, PAX5 is a known TF crucial for B cell identity and function. The inference of activities from “foot-prints” of target genes is more informative than just looking at the molecular readouts of a given TF, as an example here is the gene expression of PAX5, which is not very informative by itself:

[10]:
sc.pl.umap(adata, color='PAX5')
../_images/notebooks_dorothea_21_0.png

Exploration

With decoupler we can also see what is the mean activity per group (to show more TFs, increase the min_std parameter):

[11]:
mean_acts = dc.summarize_acts(acts, groupby='louvain', min_std=0.75)
mean_acts
[11]:
ATF2 ATF3 BACH2 CEBPB CREB1 CREM E2F7 ELK1 ESR2 FLI1 ... REL RFX5 SPI1 SRF STAT3 STAT4 TWIST1 USF2 VDR ZEB2
B cells 2.862964 0.019416 1.465129 -1.228271 0.766851 -1.464808 2.403027 2.052576 3.010636 -1.225934 ... -0.235808 13.664436 1.983893 1.009390 0.446209 -0.634329 -0.379898 1.135776 3.559168 1.004977
CD14+ Monocytes 1.999010 0.171994 2.853401 0.873312 -1.068735 -1.807254 3.214971 1.698302 5.127343 -0.795039 ... -0.484025 6.790131 3.596484 2.012941 1.400468 2.490372 -0.830063 0.463730 3.634529 2.253693
CD4 T cells 2.803933 0.874976 1.300079 -1.456408 -0.922026 -1.801688 2.631588 2.954416 3.755872 -0.997119 ... -0.048044 1.297521 1.478310 1.341562 1.228791 0.390741 -0.088309 1.126816 3.673022 1.261324
CD8 T cells 1.982264 2.649553 1.488398 -1.810813 -0.213379 -1.688924 2.406077 2.480685 4.022677 -0.927006 ... 1.185045 2.509188 2.047695 1.809774 0.565882 1.710762 -0.155667 1.347121 3.513260 0.990922
Dendritic cells 2.866896 0.180116 1.889387 -0.742983 -0.417720 -1.442676 3.151718 1.870627 3.726094 -0.932661 ... 0.012765 15.064636 3.585419 2.635822 1.193597 1.011139 -0.950940 0.472480 4.095103 2.617512
FCGR3A+ Monocytes 2.293198 0.446685 2.491488 -0.908603 -1.255402 -1.755676 3.188233 1.914761 4.890978 -0.538166 ... -0.813959 7.569202 3.998404 2.627059 2.216238 2.493336 -1.909961 0.895054 3.707708 1.937340
Megakaryocytes 0.040962 2.745303 4.102540 -2.279636 0.885389 0.852653 0.671578 -0.539585 5.077672 2.363797 ... 2.611817 2.147261 3.802499 7.186048 -1.160673 -1.038153 -2.471840 3.480007 1.433540 0.918447
NK cells 1.108534 2.085616 1.788463 -1.807453 -0.826893 -2.283551 2.114203 1.565360 4.122890 -1.014566 ... -0.035408 2.060028 2.794987 2.362190 0.249298 3.375603 -0.687441 1.471111 3.432891 0.413902

8 rows × 29 columns

We can visualize which TF is more active across cell types using seaborn:

[12]:
sns.clustermap(mean_acts, xticklabels=mean_acts.columns, vmin=-5, vmax=5, cmap='coolwarm')
plt.show()
../_images/notebooks_dorothea_25_0.png

Here we can observe other known marker TFs appering, GATA1 for megakaryocytes, RFX5 for the myeloid lineage and JUND for the lymphoid.