Multi-resolution deconvolution of spatial transcriptomics

In this tutorial, we through the steps of applying DestVI for deconvolution of 10x Visium spatial transcriptomics profiles using an accompanying single-cell RNA sequencing data.

Background:

Spatial transcriptomics technologies are currently limited, because their resolution is limited to niches (spots) of sizes well beyond that of a single cell. Although several pipelines proposed joint analysis with single-cell RNA-sequencing (scRNA-seq) to alleviate this problem they are limited to a discrete view of cell type proportion inside every spot. This limitation becomes critical in the common case where, even within a cell type, there is a continuum of cell states. We present Deconvolution of Spatial Transcriptomics profiles using Variational Inference (DestVI), a probabilistic method for multi-resolution analysis for spatial transcriptomics that explicitly models continuous variation within cell types.

Plan for this tutorial:

1. Loading the data

2. Training the single-cell model (scLVM) to learn a basis of gene expression on the scRNA-seq data

3. Training the spatial model (stLVM) to perform the deconvolution

4. Visualize the learned cell type proportions

5. Run our automated pipeline

6. Dig into the intra cell type information

7. Run cell-type specific differential expression

!pip install --quiet scvi-colab
from scvi_colab import install

install()
!pip install --quiet git+https://github.com/yoseflab/destvi_utils.git@main

import destvi_utils
import matplotlib.pyplot as plt
import numpy as np
import scanpy as sc
import scvi
from scvi.model import CondSCVI, DestVI

scvi.settings.seed = 0
print("Last run with scvi-tools version:", scvi.__version__)

Let’s download our data from a comparative study of murine lymph nodes, comparing wild-type with a stimulation after injection of a mycobacteria. We have at disposal a 10x Visium dataset as well as a matching scRNA-seq dataset from the same tissue.

sc.set_figure_params(figsize=(4, 4), frameon=False)

%config InlineBackend.print_figure_kwargs={'facecolor' : "w"}
%config InlineBackend.figure_format='retina'

url1 = "https://github.com/romain-lopez/DestVI-reproducibility/blob/master/lymph_node/deconvolution/ST-LN-compressed.h5ad?raw=true"
url2 = "https://github.com/romain-lopez/DestVI-reproducibility/blob/master/lymph_node/deconvolution/scRNA-LN-compressed.h5ad?raw=true"
out1 = "data/ST-LN-compressed.h5ad"
out2 = "data/scRNA-LN-compressed.h5ad"

Data loading & processing

First, let’s load the single-cell data. We profiled immune cells from murine lymph nodes with 10x Chromium, as a control / case study to study the immune response to exposure to a mycobacteria (refer to our paper for more info). We provide the preprocessed data from our reproducibility repository: it contains the raw counts (DestVI always takes raw counts as input).

sc_adata = sc.read(out2, backup_url=url2)

We clustered the single-cell data by major immune cell types. DestVI can resolve beyond discrete clusters, but need to work with an existing level of clustering. A rule of thumb to keep in mind while clsutering is that DestVI assumes only a single state from each cell type exists in each spot. For example, resting and inflammed monocytes cannot co-exist in one unique spot according to our assumption. Users may cluster their data so that this modeling assumption is as accurate as possible.

sc.pl.umap(sc_adata, color="broad_cell_types")

# let us filter some genes
G = 2000
sc.pp.filter_genes(sc_adata, min_counts=10)

sc_adata.layers["counts"] = sc_adata.X.copy()

sc.pp.highly_variable_genes(
    sc_adata, n_top_genes=G, subset=True, layer="counts", flavor="seurat_v3"
)

sc.pp.normalize_total(sc_adata, target_sum=10e4)
sc.pp.log1p(sc_adata)
sc_adata.raw = sc_adata

Now, let’s load the spatial data and choose a common gene subset. Users will note that having a common gene set is a prerequisite of the method.

st_adata = sc.read(out1, backup_url=url1)

st_adata.layers["counts"] = st_adata.X.copy()
st_adata.obsm["spatial"] = st_adata.obsm["location"]

sc.pp.normalize_total(st_adata, target_sum=10e4)
sc.pp.log1p(st_adata)
st_adata.raw = st_adata

# filter genes to be the same on the spatial data
intersect = np.intersect1d(sc_adata.var_names, st_adata.var_names)
st_adata = st_adata[:, intersect].copy()
sc_adata = sc_adata[:, intersect].copy()
G = len(intersect)

sc.pl.embedding(st_adata, basis="spatial", color="lymph_node", s=80)

Fit the scLVM

In order to learn cell state specific gene expression patterns, we will fit the single-cell Latent Variable Model (scLVM) to single-cell RNA sequencing data from the same tissue.

CondSCVI.setup_anndata(sc_adata, layer="counts", labels_key="broad_cell_types")

As a first step, we embed our data using a cell type conditional VAE. We pass the layer containing the raw counts and the labels key. We train this model without reweighting the loss by the cell type abundance. Training will take around 5 minutes in a Colab GPU session.

sc_model = CondSCVI(sc_adata, weight_obs=False)
sc_model.view_anndata_setup()

Anndata setup with scvi-tools version 0.16.0.

Setup via `CondSCVI.setup_anndata` with arguments:

{'labels_key': 'broad_cell_types', 'layer': 'counts'}

     Summary Statistics     
┏━━━━━━━━━━━━━━━━━━┳━━━━━━━┓
┃ Summary Stat Key ┃ Value ┃
┡━━━━━━━━━━━━━━━━━━╇━━━━━━━┩
│     n_cells      │ 14989 │
│      n_vars      │ 1888  │
│     n_labels     │  12   │
└──────────────────┴───────┘

               Data Registry                
┏━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Registry Key ┃    scvi-tools Location    ┃
┡━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│      X       │  adata.layers['counts']   │
│    labels    │ adata.obs['_scvi_labels'] │
└──────────────┴───────────────────────────┘

                         labels State Registry                         
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃        Source Location        ┃  Categories   ┃ scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['broad_cell_types'] │    B cells    │          0          │
│                               │  CD4 T cells  │          1          │
│                               │  CD8 T cells  │          2          │
│                               │  GD T cells   │          3          │
│                               │  Macrophages  │          4          │
│                               │ Migratory DCs │          5          │
│                               │   Monocytes   │          6          │
│                               │   NK cells    │          7          │
│                               │     Tregs     │          8          │
│                               │     cDC1s     │          9          │
│                               │     cDC2s     │         10          │
│                               │     pDCs      │         11          │
└───────────────────────────────┴───────────────┴─────────────────────┘

sc_model.train()

Epoch 300/300: 100%|██████████| 300/300 [05:42<00:00,  1.14s/it, loss=814, v_num=1]

sc_model.history["elbo_train"].iloc[5:].plot()
plt.show()

Note that model converges quickly. Over experimentation with the model drastically reducing the number of epochs leads to decreased performance and performance deteriorates as max_epochs<200.

Deconvolution with stLVM

As a second step, we train our deconvolution model: spatial transcriptomics Latent Variable Model (stLVM). We setup the DestVI model using the counts layer in st_adata that contains the raw counts. We then pass the trained CondSCVI model and generate a new model based on st_adata and sc_model using DestVI.from_rna_model.

The decoder network architecture will be generated from sc_model. Two neural networks are initiated for cell type proportions and gamma value amortization. Training will take around 5 minutes in a Colab GPU session.

Potential adaptations of DestVI.from_rna_model are:

1. increasing vamp_prior_p leads to less gradual changes in gamma values

2. more discretized values. Increasing l1_sparsity will lead to sparser results for cell type proportions.

3. Although we recommend using similar sequencing technology for both assays, consider changing beta_weighting_prior otherwise.

Technical Note: During inference, we adopt a variational mixture of posterior as a prior to enforce gamma values in stLVM match scLVM (see details in original publication). This empirical prior is based on cell type specific subclustering (using k-means to find vamp_prior_p clusters) of the posterior distribution in latent space for every cell.

DestVI.setup_anndata(st_adata, layer="counts")

st_model = DestVI.from_rna_model(st_adata, sc_model)
st_model.view_anndata_setup()

Anndata setup with scvi-tools version 0.16.0.

Setup via `DestVI.setup_anndata` with arguments:

{'layer': 'counts'}

     Summary Statistics     
┏━━━━━━━━━━━━━━━━━━┳━━━━━━━┓
┃ Summary Stat Key ┃ Value ┃
┡━━━━━━━━━━━━━━━━━━╇━━━━━━━┩
│     n_cells      │ 1092  │
│      n_vars      │ 1888  │
└──────────────────┴───────┘

              Data Registry              
┏━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Registry Key ┃  scvi-tools Location   ┃
┡━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━┩
│      X       │ adata.layers['counts'] │
│    ind_x     │ adata.obs['_indices']  │
└──────────────┴────────────────────────┘

st_model.train(max_epochs=2500)

Epoch 2500/2500: 100%|██████████| 2500/2500 [04:42<00:00,  8.84it/s, loss=1.95e+06, v_num=1]

Note that model converges quickly. Over experimentation with the model drastically reducing the number of epochs leads to decreased performance and we advocate against max_epochs<1000.

st_model.history["elbo_train"].iloc[10:].plot()
plt.show()

The output of DestVI has two resolution. At the broader resolution, DestVI returns the cell type proportion in every spot. At the more granular resolution, DestVI can impute cell type specific gene expression in every spot.

Cell type proportions

We extract the computed cell type proportions and display them in spatial embedding. These values are directly calculated by normalized the spot-level parameters from the stLVM model.

st_adata.obsm["proportions"] = st_model.get_proportions()

st_adata.obsm["proportions"].head(5)

	B cells	CD4 T cells	CD8 T cells	GD T cells	Macrophages	Migratory DCs	Monocytes	NK cells	Tregs	cDC1s	cDC2s	pDCs
AAACCGGGTAGGTACC-1-0	0.642821	0.009133	0.048631	0.001883	0.028526	0.028299	0.041506	0.068736	0.042653	0.056685	0.012831	0.018297
AAACCTCATGAAGTTG-1-0	0.599079	0.024904	0.073012	0.007396	0.022820	0.053057	0.055390	0.064917	0.025092	0.026119	0.037690	0.010524
AAAGACTGGGCGCTTT-1-0	0.553710	0.024247	0.098830	0.004738	0.020433	0.077447	0.031703	0.074599	0.029511	0.047068	0.036233	0.001481
AAAGGGCAGCTTGAAT-1-0	0.079860	0.135610	0.347425	0.003620	0.009415	0.097961	0.037403	0.065563	0.139114	0.052310	0.029629	0.002090
AAAGTCGACCCTCAGT-1-0	0.774980	0.018126	0.010346	0.002952	0.020493	0.030376	0.027595	0.029145	0.014909	0.037448	0.029934	0.003698

ct_list = ["B cells", "CD8 T cells", "Monocytes"]
for ct in ct_list:
    data = st_adata.obsm["proportions"][ct].values
    st_adata.obs[ct] = np.clip(data, 0, np.quantile(data, 0.99))

sc.pl.embedding(st_adata, basis="spatial", color=ct_list, cmap="Reds", s=80)

Because the inference of cell type specific gene expression is prone to error when the cell type is not present in a spot, and because the cell type proportion values estimated by stLVM are not sparse, we provide an automated way of thresholding them. For follow-up analysis we recommend checking these threshold values and adjust them for each cell type.

This part of the software is not directly available in scvi-tools, but instead in the util package destvi_utils (installable from GitHub; refer to the top of this tutorial).

ct_thresholds = destvi_utils.automatic_proportion_threshold(
    st_adata, ct_list=ct_list, kind_threshold="secondary"
)

In terms of cell type location, we observe a strong compartimentalization of the cell types in the lymph node (B cells / T cells), as expected. We also observe a differential localization of the monocytes (refer to the paper for further details).

Intra cell type information

At the heart of DestVI is a multitude of latent variables (5 per cell type per spots). We refer to them as “gamma”, and we may manually examine them for downstream analysis

# more globally, the values of the gamma are all summarized in this dictionary of data frames
for ct, g in st_model.get_gamma().items():
    st_adata.obsm[f"{ct}_gamma"] = g

st_adata.obsm["B cells_gamma"].head(5)

	0	1	2	3	4
AAACCGGGTAGGTACC-1-0	0.027521	-0.526543	0.491092	-0.180955	0.300883
AAACCTCATGAAGTTG-1-0	-0.601435	-0.306528	-2.216609	-0.012317	0.413721
AAAGACTGGGCGCTTT-1-0	-1.254444	-0.088996	-1.609583	0.299078	0.915245
AAAGGGCAGCTTGAAT-1-0	-0.082905	-0.035406	2.039065	0.107000	0.499636
AAAGTCGACCCTCAGT-1-0	0.350596	-0.700346	-0.607777	0.057218	-0.230876

Because those values may be hard to examine manually for end-users, we presented in the manuscript several methods for prioritizing the study of different cell types (based on spatially-weighted PCA and Hotspot). Below we provide the result of our automated pipeline with the spatially-weighted PCA.

More precisely, for de novo detection of spatial patterns, we study the gamma space and use a spatially-informed PCA to find the spatial axis of variation in this gamma space. We use EnrichR to functionally annotate these genes. In particular, we recover enrichment of IFN genes across monocytes as well as B cells

The function explore_gamma_space operates as follow, for each cell type individually:

1. Select all the spots with proportions beyond the magnitude threshold,

2. Calculate the spot-specific cell-type-specific embeddings gamma,

3. Calculate the first two principal vectors of those gamma values, weighted by the spatial coordinates,

4. Project all the embeddings (considered spots, and single-cell profiles) onto this 2D space,

5. Map each spot (or cell) to a specific color via its 2d coordinate, using the cmap2d package

6. Plot (A) the color of every spot in spatial coordinates (B) the color of every spot in sPC space (C) the color of every single cell in sPC space

7. Calculate genes enriched in each direction and group into pathways with EnrichR

destvi_utils.explore_gamma_space(
    st_model, sc_model, ct_list=ct_list, ct_thresholds=ct_thresholds
)

Genes associated with SpatialPC1

Positively

---------------------------------------------------------------------------------------

Birc5, Ccna2, Ube2c, Ccnb2, Hmmr, Cdca3, Cdk1, Tpx2, Cdkn3, Rrm2, Kif11, Ccnb1, Mki67, Bub1, 
Plk1, Ncaph, Cdca8, Kif2c, Kif22, Gm26532, Ncapg, Melk, Spc24, Pbk, Ska1, Neil3, Kif20a, 
Ckap2l, Cenpf, Kif4, Kif18b, Esco2, Nuf2, Cks1b, Cep55, Tyms, Uhrf1, Prc1, Hist1h3c, Stmn1, 
Tmsb4x, Cdca5, Cenpe, Top2a, Cdc20, Mxd3, Cdca2, E2f8, Fignl1, Ttk

---------------------------------------------------------------------------------------

Cell cycle*, M phase pathway*, Polo-like kinase 1 (PLK1) pathway*, DNA replication*, Aurora B
signaling*, FOXM1 transcription factor network*, Mitotic prometaphase*, Cyclin 
A/B1-associated events during G2/M transition*, Phosphorylation of Emi1*, APC/C activator 
regulation between G1/S and early anaphase*

Negatively

---------------------------------------------------------------------------------------

Ifit3, Stat1, Slfn5, Irf7, Rtp4, Zbp1, Ifi47, Usp18, Ifit3b, Oasl2, Trim30a, Ifit1, Ifit2, 
Isg15, Rsad2, Serpina3g, Parp14, Igtp, Ifi27l2a, Bst2, Malat1, Rnf213, Oas3, Gbp4, Jun, Gbp7,
Pkib, Phf11b, Ms4a4c, Isg20, Hspa1b, Ifitm3, Irf1, Ly6a, Eif2ak2, Lgals3bp, Cmpk2, Serpina3f,
Socs1, Herc6, Sdc3, Ighm, Trafd1, B2m, Hsp90aa1, Phf11a, Psmb9, Cybb, Gbp2, Gbp5

---------------------------------------------------------------------------------------

Interferon signaling*, Interferon alpha/beta signaling*, Immune system signaling by 
interferons, interleukins, prolactin, and growth hormones*, Immune system*, Interferon-gamma 
signaling pathway*, Type II interferon signaling (interferon-gamma)*, Antiviral mechanism by 
interferon-stimulated genes*, Toll-like receptor signaling pathway regulation*, Interferon 
alpha signaling regulation*, Prolactin activation of MAPK signaling*

Genes associated with SpatialPC2

Positively

---------------------------------------------------------------------------------------

Rps2, Rpsa, Rplp0, Rpl41, Rplp1, Rps12, Nme2, Ppia, Npm1, Gapdh, Mif, Hspe1, Ybx1, Impdh2, 
Actg1, Hnrnpa1, Tubb5, Pgk1, Hsp90ab1, Nme1, Atp5g1, C1qbp, Plac8, Dbi, Tuba1b, Eif5a, Prdx1,
Phgdh, Pkm, Cdk4, Ran, Ppa1, Cct8, Ranbp1, Pomp, Sec61b, Aprt, Psmb9, Snrpd1, Tkt, Nhp2, 
Psme2, Slc25a5, Npm3, Mettl1, Psma7, Fbl, Cox5a, Psma5, Ftl1

---------------------------------------------------------------------------------------

T cell receptor regulation of apoptosis*, HIV factor interactions with host*, Disease*, Gene 
expression*, HIV infection*, Protein metabolism*, Influenza infection*, Translation*, 
Cytoplasmic ribosomal proteins*, Destabilization of mRNA by AUF1 (hnRNP D0)*

Negatively

---------------------------------------------------------------------------------------

Malat1, Gm42418, Fos, Junb, Dennd4a, Fosb, Kdm6b, Pim1, Nr4a1, H2-Ab1, Dusp1, Birc3, Gm26532,
Klf4, Myadm, Ighm, Neat1, Cd69, H2-Eb1, Ahnak, Bcl11a, H2-Aa, Lyst, Rel, Zbtb20, Dusp5, 
Kcnq1ot1, Nfkbid, Klf6, Tnfaip3, Vps37b, Atp2b1, Flna, Egr1, Fcer2a, Ccr2, Gm15987, St8sia4, 
Cd79a, Zfp36, Gpr183, Vim, Srgn, Zeb2, Osbpl3, Cd74, Ddhd1, Swap70, Fam43a, Irf2bp2

---------------------------------------------------------------------------------------

Interleukin-2 signaling pathway*, T cell receptor regulation of apoptosis*, BDNF signaling 
pathway*, Interleukin-5 regulation of apoptosis*, Interleukin-1 regulation of extracellular 
matrix*, AP-1 transcription factor network*, TSH regulation of gene expression*, ATF2 
transcription factor network*, Interleukin-1 signaling pathway*, TGF-beta signaling pathway*

Genes associated with SpatialPC1

Positively

---------------------------------------------------------------------------------------

Ptma, Ppia, Npm1, Lad1, Bcat1, Rps2, Hsp90ab1, Cdca7, Utf1, Dctpp1, Tnfrsf9, Eif4a1, Set, 
Ppp1r14b, Mcm10, Fabp5, Mrpl12, Cdc45, Ran, Rpl41, Actg1, Nefh, Mpp6, Shmt1, Il2ra, Dnph1, 
Eif5a, Marcksl1, Atad3a, Cd83, Hnrnpa1, Rgs16, Ybx1, Hivep3, Dut, Ube2s, Gm26532, Eif4ebp1, 
Ncl, Ftl1, Chchd10, Mcm7, Spire1, Cct3, Timm8a1, Hmgb1, Pabpc4, Ftsj3, Stmn1, Lyar

---------------------------------------------------------------------------------------

Myc active pathway*, T cell receptor regulation of apoptosis*, Activation of the 
pre-replicative complex*, Aurora B signaling*, G2/M checkpoints*, Translation factors, 
Unwinding of DNA, Activation of mRNA upon binding of the cap-binding complex and eIFs, and 
subsequent binding to 43S, Translation, Protein metabolism

Negatively

---------------------------------------------------------------------------------------

Ifit3, Ifit1, Isg15, Rsad2, Ifit3b, Rtp4, Oasl2, Zbp1, Usp18, Irf7, Slfn5, Igtp, Isg20, Gbp7,
Oas3, Trim30a, Slfn1, Bst2, Parp14, Phf11b, Iigp1, Stat1, Gbp2, Herc6, Ms4a4b, Cmpk2, Ddx60, 
Rnf213, Ly6a, Lgals3bp, Ifi47, Gbp4, Eif2ak2, Ifih1, Epsti1, Trafd1, AW112010, Plac8, B2m, 
Ifi27l2a, Gbp8, Ifit2, Psmb9, Parp12, Ly6c2, Cd274, Ms4a6b, Phf11a, Gbp5, Ms4a4c

---------------------------------------------------------------------------------------

Interferon signaling*, Immune system signaling by interferons, interleukins, prolactin, and 
growth hormones*, Interferon alpha/beta signaling*, Immune system*, Interferon-gamma 
signaling pathway*, Type II interferon signaling (interferon-gamma)*, Antiviral mechanism by 
interferon-stimulated genes*, TRAF3-dependent IRF activation pathway*, RIG-I-like receptor 
signaling pathway, RIG-I/MDA5-mediated induction of interferon-alpha/beta pathways

Genes associated with SpatialPC2

Positively

---------------------------------------------------------------------------------------

Tcrg-C4, Tcrg-C2, Ikzf2, Trdc, Anxa2, Malat1, Tcrg-C1, S100a6, Lgals1, Actb, Ms4a4b, Klrc2, 
Ahnak, Sox4, Cd160, Rbpms, Klrc1, Ccl5, Fcer1g, Ctla2a, Spry2, Klrk1, Zyx, Ahr, Ctla2b, 
Slfn5, Ar, Il2rb, Prex1, Xcl1, Klre1, Itgb1, Klra1, Tmsb4x, Cd7, Sept11, Slamf7, Myo1f, 
Klra7, S100a4, Serpina3g, Fyn, Itgal, Spn, Bcl2, Klf6, Nkg7, Dok2, Rgs2, Klrb1c

---------------------------------------------------------------------------------------

Interleukin-2 signaling pathway*, Cell surface interactions at the vascular wall*, T cell 
receptor regulation of apoptosis*, Natural killer cell-mediated cytotoxicity*, Platelet 
adhesion to exposed collagen*, Immunoregulatory interactions between a lymphoid and a 
non-lymphoid cell*, Focal adhesion*, Ras-independent pathway in NK cell-mediated 
cytotoxicity*, Myc repressed pathway*, Integrin-mediated cell adhesion*

Negatively

---------------------------------------------------------------------------------------

Rplp0, Rpsa, Rplp1, Rps2, Rpl41, Rps12, Npm1, Cd8b1, Nme2, Cd8a, Gbp2, Dapl1, Ly6a, Igtp, 
Npc2, Stat1, Plac8, Rgs10, AW112010, Psme2, Hspe1, Impdh2, Hspa8, Gbp4, Fth1, Igfbp4, Iigp1, 
Ifi47, Cd274, Trac, Psmb9, Hnrnpa1, Lef1, Il7r, Ppa1, B2m, Sec61b, Fbl, Klk8, Irf1, Socs3, 
Ppia, Tubb5, Gbp8, Sp140, Lcn4, Dnajc15, Art2b, Ctss, Gbp5

---------------------------------------------------------------------------------------

T cell receptor regulation of apoptosis*, Interferon-gamma signaling pathway*, Disease*, 
Immune system*, Antigen processing and presentation*, Translation*, Interferon signaling*, 
Cytoplasmic ribosomal proteins*, Immune system signaling by interferons, interleukins, 
prolactin, and growth hormones*, Influenza viral RNA transcription and replication*

Genes associated with SpatialPC1

Positively

---------------------------------------------------------------------------------------

Rplp1, Spdl1, Mcm10, Ppia, Chaf1a, Tuba1c, Hnrnpll, Nr4a3, S100a4, Rps2, Hnrnpab, Tarm1, 
Tuba1a, Dhfr, Rpl41, Rps12, Espl1, Odc1, Kntc1, Dbn1, Exo1, Ube2t, Rad51ap1, Dscc1, Cenpp, 
Il1b, Mif, Fabp5, Atox1, Klf4, Eif1ax, Cdk14, Fen1, Cbx5, Rad21, Fosl2, Smc1a, Shmt1, Cdk4, 
Nhp2, Itgax, Eif4a1, Eif4ebp1, Fam149a, Avpi1, Procr, Plxnc1, Sema7a, Ncapg, Neat1

---------------------------------------------------------------------------------------

Cell cycle*, Translation*, Activation of mRNA upon binding of the cap-binding complex and 
eIFs, and subsequent binding to 43S*, DNA replication*, Protein metabolism*, M phase 
pathway*, Chromosome maintenance*, Cytoplasmic ribosomal proteins*, Cap-dependent translation
initiation*, Myc active pathway*

Negatively

---------------------------------------------------------------------------------------

Ifi204, Ms4a6b, Oasl2, Igtp, Ifi47, Fcgr1, Cxcl10, Ifit3, Rtp4, Apoe, Ifit1, Ifit2, Oas3, 
Gbp7, Slfn5, Isg15, Ifit3b, Rnf213, Stat1, Sp140, Gbp2, Irf8, Mpeg1, Irf7, Zbp1, Scimp, Bst2,
Usp18, Slfn1, Gbp3, Ifi205, Cybb, Cmpk2, Gbp5, Ms4a4c, Ctss, Trafd1, Ccl12, Ifih1, Ms4a7, 
Tiam1, Ccr2, Cd7, Vcan, Pbx1, Flt1, Mertk, Lpl, Pgap2, Irf1

---------------------------------------------------------------------------------------

Interferon signaling*, Interferon alpha/beta signaling*, Immune system signaling by 
interferons, interleukins, prolactin, and growth hormones*, Immune system*, Interferon-gamma 
signaling pathway*, Type II interferon signaling (interferon-gamma)*, Innate immune system*, 
Antiviral mechanism by interferon-stimulated genes*, RIG-I-like receptor signaling pathway*, 
RIG-I/MDA5-mediated induction of interferon-alpha/beta pathways*

Genes associated with SpatialPC2

Positively

---------------------------------------------------------------------------------------

Gfpt1, Nedd4, Dusp5, Fosb, Hsph1, Nuak1, Uhrf1, AA467197, Ccnd2, Ddx60, Gm42418, Kif23, 
Sirpa, Ifit2, Egr1, Dusp1, Zc3h12c, Rnf213, Kdr, B2m, Pmaip1, Slc6a6, Apol7c, Cd83, Tnfrsf1b,
Hells, Fscn1, Ifit3, Rasal2, Rp2, Kctd14, Igf1, Rel, Dab2, Hpn, Eaf2, Cxcr3, Synj1, Klf4, 
Pim1, Zfp366, Rell1, Prkcd, Sgpl1, Arl5a, Ear2, Stmn1, Mical3, Tmem26, Marcksl1

---------------------------------------------------------------------------------------

T cell receptor regulation of apoptosis*, Interferon signaling*, Interleukin-2 signaling 
pathway*, BDNF signaling pathway*, Immune system signaling by interferons, interleukins, 
prolactin, and growth hormones*, p53 activity regulation, Apoptosis, SRF and microRNAs in 
smooth muscle differentiation and proliferation, EGF/EGFR signaling pathway, Interleukin-5 
regulation of apoptosis

Negatively

---------------------------------------------------------------------------------------

Atp6v0c, Sec61b, Nme2, Smdt1, Taldo1, Atp5g1, Tyrobp, Capg, Snrpd1, Hmgb1, Nme1, Ly86, Prdx4,
Bax, S100a6, Ppia, Fcer1g, Lyz2, Anxa5, Hexa, Slc25a5, Cd68, Eif4ebp1, Aprt, Tmed3, Guca1a, 
Ndufab1, Rplp0, Ramp1, Ndufc2, Hint1, Tmem256, Cnih4, Anxa2, Ifitm2, Cd302, Uqcc2, Tifab, 
Gapdh, Tpd52, Cyb5a, Etfb, Ifi30, Psme2, Pgam1, Dbi, Tspo, Pgk1, Snrpa1, Tmem14c

---------------------------------------------------------------------------------------

T cell receptor regulation of apoptosis*, Metabolism*, Glycolysis*, Respiratory electron 
transport, ATP biosynthesis by chemiosmotic coupling, and heat production by uncoupling 
proteins*, Gluconeogenesis*, Prostaglandin biosynthesis and regulation*, Huntington's 
disease*, Electron transport chain*, Tricarboxylic acid (TCA) cycle and respiratory electron 
transport*, Carbohydrate metabolism*

We anticipate this to be a valuable ressource for formulating scientific hypotheses from ST data.

Example with B cells; and differential expression

First, we display the genes identified via the pipeline as well as Hotspot (see manuscript), using the B-cell-specific gene expression values imputed by DestVI.

plt.figure(figsize=(8, 8))

ct_name = "B cells"
gene_name = ["Ifit3", "Ifit3b", "Ifit1", "Isg15", "Oas3", "Usp18", "Isg20"]

# we must filter spots with low abundance (consult the paper for an automatic procedure)
indices = np.where(st_adata.obsm["proportions"][ct_name].values > 0.2)[0]

# impute genes and combine them
specific_expression = np.sum(
    st_model.get_scale_for_ct(ct_name, indices=indices)[gene_name], 1
)
specific_expression = np.log(1 + 1e4 * specific_expression)

# plot (i) background (ii) g
plt.scatter(
    st_adata.obsm["location"][:, 0], st_adata.obsm["location"][:, 1], alpha=0.05
)
plt.scatter(
    st_adata.obsm["location"][indices][:, 0],
    st_adata.obsm["location"][indices][:, 1],
    c=specific_expression,
    s=10,
    cmap="Reds",
)
plt.colorbar()
plt.title(f"Imputation of {gene_name} in {ct_name}")
plt.show()

Second, we apply a Kolmogorov-Smirnov test on the generated counts to study the differential expression of B cells in the exposed lymph nodes, between the interfollicular area (IFA) and the rest. We display the identified IFN genes in a Volcano plot and see significant upregulation in the IFA area of exposed lymph nodes.

ct = "B cells"
imputation = st_model.get_scale_for_ct(ct)
color = np.log(1 + 1e5 * imputation["Ifit3"].values)
threshold = 4

mask = np.logical_and(
    np.logical_or(st_adata.obs["LN"] == "TC", st_adata.obs["LN"] == "BD"),
    color > threshold,
).values

mask2 = np.logical_and(
    np.logical_or(st_adata.obs["LN"] == "TC", st_adata.obs["LN"] == "BD"),
    color < threshold,
).values

_ = destvi_utils.de_genes(
    st_model, mask=mask, mask2=mask2, threshold=ct_thresholds[ct], ct=ct, key="IFN_rich"
)

display(st_adata.uns["IFN_rich"]["de_results"].head(10))

destvi_utils.plot_de_genes(
    st_adata,
    interesting_genes=["Ifit3", "Ifit3b", "Ifit1", "Isg15", "Oas3", "Usp18", "Isg20"],
    key="IFN_rich",
)

	log2FC	score	pval
Irf7	2.901810	0.933216	0.0
Zbp1	2.736697	0.923558	0.0
Rnf213	1.877064	0.905463	0.0
Trim30a	1.883341	0.896659	0.0
Trafd1	1.855225	0.895122	0.0
Ifit3	4.828797	0.887100	0.0
Usp18	3.743639	0.887058	0.0
Parp14	1.977952	0.881785	0.0
Rtp4	2.848005	0.876900	0.0
Oas3	3.433762	0.869146	0.0