Note

This page was generated from DestVI_tutorial.ipynb. Interactive online version: Colab badge. Some tutorial content may look better in light mode.

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

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.

[ ]:
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt

from scvi.model import CondSCVI, DestVI

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")
../../_images/tutorials_notebooks_DestVI_tutorial_9_0.png
[ ]:
# 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)
../../_images/tutorials_notebooks_DestVI_tutorial_15_0.png

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()
GPU available: True, used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
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()
../../_images/tutorials_notebooks_DestVI_tutorial_22_0.png

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)
GPU available: True, used: True
TPU available: False, using: 0 TPU cores
IPU available: False, using: 0 IPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
/usr/local/lib/python3.7/dist-packages/pytorch_lightning/trainer/data_loading.py:433: UserWarning: The number of training samples (9) is smaller than the logging interval Trainer(log_every_n_steps=10). Set a lower value for log_every_n_steps if you want to see logs for the training epoch.
  f"The number of training samples ({self.num_training_batches}) is smaller than the logging interval"
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()
../../_images/tutorials_notebooks_DestVI_tutorial_30_0.png

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)
../../_images/tutorials_notebooks_DestVI_tutorial_36_0.png

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).

[ ]:
import destvi_utils
Creating directory /root/.config/bioservices
/usr/local/lib/python3.7/dist-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
  import pandas.util.testing as tm
[ ]:
ct_thresholds = destvi_utils.automatic_proportion_threshold(st_adata, ct_list=ct_list, kind_threshold='secondary')
/usr/local/lib/python3.7/dist-packages/numba/np/ufunc/parallel.py:363: NumbaWarning: The TBB threading layer requires TBB version 2019.5 or later i.e., TBB_INTERFACE_VERSION >= 11005. Found TBB_INTERFACE_VERSION = 9107. The TBB threading layer is disabled.
  warnings.warn(problem)
100%|██████████| 100/100 [00:00<00:00, 232.41it/s]
../../_images/tutorials_notebooks_DestVI_tutorial_39_1.png
100%|██████████| 100/100 [00:00<00:00, 12786.73it/s]
../../_images/tutorials_notebooks_DestVI_tutorial_39_3.png
100%|██████████| 100/100 [00:00<00:00, 6621.78it/s]
../../_images/tutorials_notebooks_DestVI_tutorial_39_5.png

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["{}_gamma".format(ct)] = 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)
/usr/local/lib/python3.7/dist-packages/cmap2d/util.py:32: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  bary_coords = la.lstsq(self._a,b)[0]
../../_images/tutorials_notebooks_DestVI_tutorial_44_1.png
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*



../../_images/tutorials_notebooks_DestVI_tutorial_44_30.png
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*



../../_images/tutorials_notebooks_DestVI_tutorial_44_59.png
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()
../../_images/tutorials_notebooks_DestVI_tutorial_48_0.png

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
../../_images/tutorials_notebooks_DestVI_tutorial_50_1.png