Note

This page was generated from totalVI.ipynb. Interactive online version: Colab badge.

CITE-seq analysis with totalVI

With totalVI, we can produce a joint latent representation of cells, denoised data for both protein and RNA, integrate datasets, and compute differential expression of RNA and protein. Here we demonstrate this functionality with an integrated analysis of PBMC10k and PBMC5k, datasets of peripheral blood mononuclear cells publicly available from 10X Genomics subset to the 14 shared proteins between them. The same pipeline would generally be used to analyze a single CITE-seq dataset.

If you use totalVI, please consider citing:

  • Gayoso, A., Steier, Z., Lopez, R., Regier, J., Nazor, K. L., Streets, A., & Yosef, N. (2021). Joint probabilistic modeling of single-cell multi-omic data with totalVI. Nature Methods, 18(3), 272-282.

[1]:
import sys

#if branch is stable, will install via pypi, else will install from source
branch = "stable"
IN_COLAB = "google.colab" in sys.modules

if IN_COLAB and branch == "stable":
    !pip install --quiet scvi-tools[tutorials]
elif IN_COLAB and branch != "stable":
    !pip install --quiet --upgrade jsonschema
    !pip install --quiet git+https://github.com/yoseflab/scvi-tools@$branch#egg=scvi-tools[tutorials]
     |████████████████████████████████| 61kB 6.3MB/s
ERROR: nbclient 0.5.3 has requirement jupyter-client>=6.1.5, but you'll have jupyter-client 5.3.5 which is incompatible.
  Installing build dependencies ... done
  Getting requirements to build wheel ... done
    Preparing wheel metadata ... done
     |████████████████████████████████| 204kB 15.8MB/s
     |████████████████████████████████| 133kB 23.5MB/s
     |████████████████████████████████| 245kB 29.9MB/s
     |████████████████████████████████| 81kB 10.8MB/s
     |████████████████████████████████| 634kB 30.0MB/s
     |████████████████████████████████| 808kB 64.3MB/s
     |████████████████████████████████| 1.4MB 48.8MB/s
     |████████████████████████████████| 3.2MB 56.4MB/s
     |████████████████████████████████| 8.8MB 39.9MB/s
     |████████████████████████████████| 51kB 974kB/s
     |████████████████████████████████| 10.3MB 22.2MB/s
     |████████████████████████████████| 51kB 7.7MB/s
     |████████████████████████████████| 112kB 53.0MB/s
     |████████████████████████████████| 645kB 53.8MB/s
     |████████████████████████████████| 276kB 53.5MB/s
     |████████████████████████████████| 829kB 50.0MB/s
     |████████████████████████████████| 112kB 62.7MB/s
     |████████████████████████████████| 81kB 10.7MB/s
     |████████████████████████████████| 1.3MB 45.1MB/s
     |████████████████████████████████| 71kB 10.4MB/s
     |████████████████████████████████| 1.2MB 40.3MB/s
     |████████████████████████████████| 143kB 60.6MB/s
     |████████████████████████████████| 296kB 66.2MB/s
  Building wheel for scvi-tools (PEP 517) ... done
  Building wheel for loompy (setup.py) ... done
  Building wheel for future (setup.py) ... done
  Building wheel for numpy-groupies (setup.py) ... done
  Building wheel for sinfo (setup.py) ... done
  Building wheel for umap-learn (setup.py) ... done
  Building wheel for pynndescent (setup.py) ... done

Imports and data loading

[2]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

import scvi
import scanpy as sc

sc.set_figure_params(figsize=(4, 4))
/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)

This dataset was filtered as described in the totalVI manuscript (low quality cells, doublets, lowly expressed genes, etc.)

[3]:
adata = scvi.data.pbmcs_10x_cite_seq(run_setup_anndata=False)
adata.layers["counts"] = adata.X.copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata
INFO     Downloading file at data/pbmc_10k_protein_v3.h5ad
Downloading...: 24938it [00:00, 117848.96it/s]
INFO     Downloading file at data/pbmc_5k_protein_v3.h5ad
Downloading...: 100%|██████████| 18295/18295.0 [00:00<00:00, 84501.31it/s]
Observation names are not unique. To make them unique, call `.obs_names_make_unique`.
[4]:
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=4000,
    flavor="seurat_v3",
    batch_key="batch",
    subset=True,
    layer="counts"
)
Observation names are not unique. To make them unique, call `.obs_names_make_unique`.
Observation names are not unique. To make them unique, call `.obs_names_make_unique`.
[5]:
scvi.data.setup_anndata(
    adata,
    layer="counts",
    batch_key="batch",
    protein_expression_obsm_key="protein_expression"
)
INFO     Using batches from adata.obs["batch"]
INFO     No label_key inputted, assuming all cells have same label
INFO     Using data from adata.layers["counts"]
INFO     Computing library size prior per batch
INFO     Using protein expression from adata.obsm['protein_expression']
INFO     Using protein names from columns of adata.obsm['protein_expression']
INFO     Successfully registered anndata object containing 10849 cells, 4000 vars, 2 batches,
         1 labels, and 14 proteins. Also registered 0 extra categorical covariates and 0
         extra continuous covariates.
INFO     Please do not further modify adata until model is trained.

Prepare and run model

[6]:
vae = scvi.model.TOTALVI(adata, latent_distribution="normal")
[7]:
vae.train()
GPU available: True, used: True
TPU available: False, using: 0 TPU cores
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]
Epoch 262/400:  66%|██████▌   | 262/400 [04:02<01:59,  1.15it/s, loss=1.22e+03, v_num=1]Epoch   262: reducing learning rate of group 0 to 2.4000e-03.
Epoch 298/400:  74%|███████▍  | 298/400 [04:34<01:27,  1.17it/s, loss=1.22e+03, v_num=1]Epoch   298: reducing learning rate of group 0 to 1.4400e-03.
Epoch 364/400:  91%|█████████ | 364/400 [05:32<00:33,  1.08it/s, loss=1.21e+03, v_num=1]Epoch   364: reducing learning rate of group 0 to 8.6400e-04.
Epoch 378/400:  94%|█████████▍| 378/400 [05:45<00:20,  1.09it/s, loss=1.22e+03, v_num=1]
[8]:
plt.plot(vae.history["elbo_train"], label="train")
plt.plot(vae.history["elbo_validation"], label="validation")
plt.title("Negative ELBO over training epochs")
plt.ylim(1200, 1400)
plt.legend()
[8]:
<matplotlib.legend.Legend at 0x7fbb8119f990>
../../_images/user_guide_notebooks_totalVI_11_1.png

Analyze outputs

We use Scanpy for clustering and visualization after running totalVI. It’s also possible to save totalVI outputs for an R-based workflow. First, we store the totalVI outputs in the appropriate slots in AnnData.

[9]:
adata.obsm["X_totalVI"] = vae.get_latent_representation()

rna, protein = vae.get_normalized_expression(
    n_samples=25,
    return_mean=True,
    transform_batch=["PBMC10k", "PBMC5k"]
)

adata.layers["denoised_rna"], adata.obsm["denoised_protein"] = rna, protein

adata.obsm["protein_foreground_prob"] = vae.get_protein_foreground_probability(
    n_samples=25,
    return_mean=True,
    transform_batch=["PBMC10k", "PBMC5k"]
)
parsed_protein_names = [p.split("_")[0] for p in adata.obsm["protein_expression"].columns]
adata.obsm["protein_foreground_prob"].columns = parsed_protein_names

Now we can compute clusters and visualize the latent space.

[10]:
sc.pp.neighbors(adata, use_rep="X_totalVI")
sc.tl.umap(adata, min_dist=0.4)
sc.tl.leiden(adata, key_added="leiden_totalVI")
[11]:
sc.pl.umap(
    adata,
    color=["leiden_totalVI", "batch"],
    frameon=False,
    ncols=1,
)
... storing 'batch' as categorical
../../_images/user_guide_notebooks_totalVI_16_1.png

To visualize protein values on the umap, we make a temporary protein adata object. We have to copy over the umap from the original adata object.

[12]:
pro_adata = sc.AnnData(adata.obsm["protein_expression"].copy(), obs=adata.obs)
sc.pp.log1p(pro_adata)
# Keep log normalized data in raw
pro_adata.raw = pro_adata
pro_adata.X = adata.obsm["denoised_protein"]
# these are cleaner protein names -- "_TotalSeqB" removed
pro_adata.var["protein_names"] = parsed_protein_names
pro_adata.obsm["X_umap"] = adata.obsm["X_umap"]
pro_adata.obsm["X_totalVI"] = adata.obsm["X_totalVI"]

Observation names are not unique. To make them unique, call `.obs_names_make_unique`.
[13]:
names = adata.obsm["protein_foreground_prob"].columns
for p in names:
    pro_adata.obs["{}_fore_prob".format(p)] = adata.obsm["protein_foreground_prob"].loc[:, p]

Visualize denoised protein values

[14]:
sc.pl.umap(
    pro_adata,
    color=pro_adata.var_names,
    gene_symbols="protein_names",
    ncols=3,
    vmax="p99",
    use_raw=False,
    frameon=False,
    wspace=0.1
)
../../_images/user_guide_notebooks_totalVI_21_0.png

Visualize probability of foreground

Here we visualize the probability of foreground for each protein and cell (projected on UMAP). Some proteins are easier to disentangle than others. Some proteins end up being “all background”. For example, CD15 does not appear to be captured well, when looking at the denoised values above we see little localization in the monocytes.

Note

While the foreground probability could theoretically be used to identify cell populations, we recommend using the denoised protein expression, which accounts for the foreground/background probability, but preserves the dynamic range of the protein measurements. Consequently, the denoised values are on the same scale as the raw data and it may be desirable to take a transformation like log or square root.

By viewing the foreground probability, we can get a feel for the types of cells in our dataset. For example, it’s very easy to see a population of monocytes based on the CD14 foregroud probability.

[15]:
sc.pl.umap(
    pro_adata,
    color=["{}_fore_prob".format(p) for p in parsed_protein_names],
    ncols=3,
    color_map="cividis",
    frameon=False,
    wspace=0.1
)
../../_images/user_guide_notebooks_totalVI_25_0.png

Differential expression

Here we do a one-vs-all DE test, where each cluster is tested against all cells not in that cluster. The results for each of the one-vs-all tests is concatenated into one DataFrame object. Inividual tests can be sliced using the “comparison” column. Genes and proteins are included in the same DataFrame.

Important

We do not recommend using totalVI denoised values in other differential expression tools, as denoised values are a summary of a random quantity. The totalVI DE test takes into account the full uncertainty of the denoised quantities.

[16]:
de_df = vae.differential_expression(
        groupby="leiden_totalVI",
        delta=0.5,
        batch_correction=True
)
de_df.head(5)
DE...: 100%|██████████| 18/18 [00:36<00:00,  2.02s/it]
[16]:
proba_de proba_not_de bayes_factor scale1 scale2 pseudocounts delta lfc_mean lfc_median lfc_std lfc_min lfc_max raw_mean1 raw_mean2 non_zeros_proportion1 non_zeros_proportion2 raw_normalized_mean1 raw_normalized_mean2 is_de_fdr_0.05 comparison
S100A12 0.9978 0.0022 6.117091 0.003423 5.052456e-05 0.0 0.5 8.857993 9.102179 2.531390 -2.200211 17.580215 9.887289 0.151105 0.965858 0.046923 37.031631 0.452420 True 0 vs Rest
MARC1 0.9974 0.0026 5.949637 0.000124 1.637846e-06 0.0 0.5 8.603086 8.766019 2.684580 -5.585446 17.082640 0.323324 0.004871 0.248457 0.004633 1.257425 0.015061 True 0 vs Rest
CLEC4D 0.9968 0.0032 5.741396 0.000099 1.667533e-06 0.0 0.5 8.163298 8.302477 2.470065 -2.996992 17.291874 0.268614 0.004395 0.221308 0.003920 1.017459 0.014606 True 0 vs Rest
LIN7A 0.9966 0.0034 5.680571 0.000046 9.869541e-07 0.0 0.5 7.891823 8.094858 2.475345 -4.837774 18.115129 0.133690 0.002376 0.116002 0.002376 0.456346 0.006435 True 0 vs Rest
PADI2 0.9964 0.0036 5.623212 0.000084 1.486880e-06 0.0 0.5 7.693967 7.775435 2.694114 -3.209021 15.885664 0.220897 0.004277 0.178939 0.004158 0.803408 0.012067 True 0 vs Rest

Now we filter the results such that we retain features above a certain Bayes factor (which here is on the natural log scale) and genes with greater than 10% non-zero entries in the cluster of interest.

[17]:
filtered_pro = {}
filtered_rna = {}
cats = adata.obs.leiden_totalVI.cat.categories
for i, c in enumerate(cats):
    cid = "{} vs Rest".format(c)
    cell_type_df = de_df.loc[de_df.comparison == cid]
    cell_type_df = cell_type_df.sort_values("lfc_median", ascending=False)

    cell_type_df = cell_type_df[cell_type_df.lfc_median > 0]

    pro_rows = cell_type_df.index.str.contains('TotalSeqB')
    data_pro = cell_type_df.iloc[pro_rows]
    data_pro = data_pro[data_pro["bayes_factor"] > 0.7]

    data_rna = cell_type_df.iloc[~pro_rows]
    data_rna = data_rna[data_rna["bayes_factor"] > 3]
    data_rna = data_rna[data_rna["non_zeros_proportion1"] > 0.1]

    filtered_pro[c] = data_pro.index.tolist()[:3]
    filtered_rna[c] = data_rna.index.tolist()[:2]

We can also use general scanpy visualization functions

[18]:
sc.tl.dendrogram(adata, groupby="leiden_totalVI", use_rep="X_totalVI")
sc.tl.dendrogram(pro_adata, groupby="leiden_totalVI", use_rep="X_totalVI")
/usr/lib/python3.7/importlib/_bootstrap.py:219: RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject
  return f(*args, **kwds)
[19]:
sc.pl.dotplot(
    adata,
    filtered_rna,
    groupby="leiden_totalVI",
    dendrogram=True,
    standard_scale="var",
    swap_axes=True
)
../../_images/user_guide_notebooks_totalVI_34_0.png

Matrix plot displays totalVI denoised protein expression per leiden cluster.

[20]:
sc.pl.matrixplot(
    pro_adata,
    pro_adata.var["protein_names"],
    groupby="leiden_totalVI",
    gene_symbols="protein_names",
    dendrogram=True,
    swap_axes=True,
    use_raw=False, # use totalVI denoised
    cmap="Greens",
    standard_scale="var"
)
../../_images/user_guide_notebooks_totalVI_36_0.png

This is a selection of some of the markers that turned up in the RNA DE test.

[21]:
sc.pl.umap(
    adata,
    color=[
           "leiden_totalVI",
           "IGHD",
           "FCER1A",
           "SCT",
           "GZMH",
           "NOG",
           "FOXP3",
           "CD8B",
           "C1QA",
           "SIGLEC1",
           "XCL2",
           "GZMK",
           ],
    legend_loc="on data",
    frameon=False,
    ncols=3,
    layer="denoised_rna",
    wspace=0.1
)
../../_images/user_guide_notebooks_totalVI_38_0.png