<jupyter_start><jupyter_text>Annotation with CellAssign  Assigning single-cell RNA-seq data to known cell types CellAssign is a probabilistic model that uses prior knowledge of cell-type marker genes to annotate scRNA data into predefined cell types. Unlike other methods for assigning cell types, CellAssign does not require labeled single cell data and only needs to know whether or not each given gene is a marker of each cell type. The original paper and R code are linked below.Paper: [Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling, *Nature Methods 2019*](https://www.nature.com/articles/s41592-019-0529-1 )Code: https://github.com/Irrationone/cellassign This notebook will demonstrate how to use CellAssign on follicular lymphoma and HGSC scRNA data.<jupyter_code>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]<jupyter_output>[K     |████████████████████████████████| 184kB 18.5MB/s 
[K     |████████████████████████████████| 122kB 48.2MB/s 
[K     |████████████████████████████████| 829kB 58.6MB/s 
[K     |████████████████████████████████| 634kB 51.3MB/s 
[K     |████████████████████████████████| 81kB 11.7MB/s 
[K     |████████████████████████████████| 204kB 58.3MB/s 
[K     |████████████████████████████████| 245kB 51.3MB/s 
[K     |████████████████████████████████| 8.7MB 29.6MB/s 
[K     |████████████████████████████████| 2.4MB 41.4MB/s 
[K     |████████████████████████████████| 51kB 7.9MB/s 
[K     |████████████████████████████████| 10.3MB 25.8MB/s 
[K     |████████████████████████████████| 3.1MB 55.3MB/s 
[K     |████████████████████████████████| 276kB 61.3MB/s 
[K     |████████████████████████████████| 112kB 56.2MB/s 
[K     |████████████████████████████████| 829kB 53.2MB/s 
[K     |████████████████████████████████| 51kB 9.1MB/s 
[K     |████████████████████████████████| 112kB 63.8MB/s 
[...]<jupyter_text>To demonstrate CellAssign, we use the data from the original publication, which we converted into h5ad format. The data are originally available from here:https://zenodo.org/record/3372746<jupyter_code>import gdown
url = 'https://drive.google.com/uc?id=10l6m2KKKioCZnQlRHomheappHh-jTFmx'
output = 'sce_follicular_annotated_final.h5ad'
gdown.download(url, output, quiet=False)

url = 'https://drive.google.com/uc?id=1Pae7VEcoZbKRvtllGAEWG4SOLWSjjtCO'
output = 'sce_hgsc_annotated_final.h5ad'
gdown.download(url, output, quiet=False)

url = 'https://drive.google.com/uc?id=1Mk5uPdnPC4IMRnuG5N4uFvypT8hPdJ74'
output = 'HGSC_celltype.csv'
gdown.download(url, output, quiet=False)

url = 'https://drive.google.com/uc?id=1tJSOI9ve0i78WmszMLx2ul8F8tGycBTd'
output = 'FL_celltype.csv'
gdown.download(url, output, quiet=False)
import scvi
import scanpy as sc
import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Follicular Lymphoma DataLoad follicular lymphoma data and marker gene matrix (see Supplementary Table 2 from the original paper).<jupyter_code>adata = sc.read("sce_follicular_annotated_final.h5ad")
adata.var_names_make_unique()
adata.obs_names_make_unique()
marker_gene_mat = pd.read_csv('FL_celltype.csv', index_col=0)<jupyter_output><empty_output><jupyter_text>Create and fit CellAssign model The anndata object and cell type marker matrix should contain the same genes, so we index into `adata` to include only the genes from `marker_gene_mat`.<jupyter_code>bdata = adata[:, marker_gene_mat.index].copy()
scvi.data.setup_anndata(bdata)<jupyter_output>[34mINFO    [0m No batch_key inputted, assuming all cells are same batch                            
[34mINFO    [0m No label_key inputted, assuming all cells have same label                           
[34mINFO    [0m Using data from adata.X                                                             
[34mINFO    [0m Computing library size prior per batch                                              
[31mWARNING [0m This dataset has some empty cells, this might fail inference.Data should be filtered
         with `scanpy.pp.filter_cells[1m([0m[1m)[0m`                                                     
[34mINFO    [0m Successfully registered anndata object containing [1;34m9156[0m cells, [1;34m24[0m vars, [1;34m1[0m batches, [1;34m1[0m 
         labels, and [1;34m0[0m proteins. Also registered [1;34m0[0m extra categorical covariates and [1;34m0[0m extra  
         continuous covariates.                                                              
[34mIN[...]<jupyter_text>Then we setup anndata and initialize a `CellAssign` model. Here we set the `size_factor_key` to "size_factor", which is a column in `bdata.obs`. A size factor may be defined manually as UMI count and should not be placed on the log scale, as the model will do this manually.For example,```pythonadata.obs["size_factor"] = adata.X.sum(1)```<jupyter_code>from scvi.external import CellAssign
model = CellAssign(bdata, marker_gene_mat, size_factor_key="size_factor")
model.train()<jupyter_output>GPU available: True, used: True
TPU available: None, using: 0 TPU cores<jupyter_text>Inspecting the convergence:<jupyter_code>model.history["elbo_validation"].plot()<jupyter_output><empty_output><jupyter_text>Predict and plot assigned cell types Predict the soft cell type assignment probability for each cell.<jupyter_code>predictions = model.predict()
predictions.head()<jupyter_output><empty_output><jupyter_text>We can visualize the probabilities of assignment with a heatmap that returns the probability matrix for each cell and cell type.<jupyter_code>sns.clustermap(predictions, cmap="viridis")<jupyter_output>/usr/local/lib/python3.7/dist-packages/seaborn/matrix.py:649: UserWarning: Clustering large matrix with scipy. Installing `fastcluster` may give better performance.
  warnings.warn(msg)<jupyter_text>We then create a UMAP plot labeled by maximum probability assignments from the CellAssign model. The left plot contains the true cell types and the right plot contains our model's predictions.<jupyter_code>bdata.obs["scvi-tools predictions"] = predictions.idxmax(axis=1).values
# celltype is the original CellAssign prediction
sc.pl.umap(
    bdata,
    color=["celltype", "scvi-tools predictions"],
    frameon=False,
    ncols=1
)<jupyter_output><empty_output><jupyter_text>Model reproducibility We see that the scvi-tools implementation highly reproduces the original implementation's predictions.<jupyter_code>df = bdata.obs
confusion_matrix = pd.crosstab(
    df["scvi-tools predictions"],
    df["celltype"],
    rownames=["scvi-tools predictions"],
    colnames=["Original predictions"],
)
confusion_matrix /= confusion_matrix.sum(1).ravel().reshape(-1, 1)
fig, ax = plt.subplots(figsize=(5, 4))
sns.heatmap(
    confusion_matrix,
    cmap=sns.diverging_palette(245, 320, s=60, as_cmap=True),
    ax=ax,
    square=True,
    cbar_kws=dict(shrink=0.4, aspect=12),
)<jupyter_output><empty_output><jupyter_text>HGSC Data We can repeat the same process for HGSC data.<jupyter_code>hgsc_adata = scvi.data.read_h5ad("sce_hgsc_annotated_final.h5ad")
hgsc_adata.var_names_make_unique()
hgsc_adata.obs_names_make_unique()
marker_gene_mat_hgsc = pd.read_csv('HGSC_celltype.csv', index_col=0)<jupyter_output><empty_output><jupyter_text>Create and fit CellAssign model<jupyter_code>hgsc_bdata = hgsc_adata[:, marker_gene_mat_hgsc.index].copy()
scvi.data.setup_anndata(hgsc_bdata)
from scvi.external import CellAssign
model_hgsc = CellAssign(hgsc_bdata, marker_gene_mat_hgsc, "size_factor")
model_hgsc.train()
model.history["elbo_validation"].plot()<jupyter_output><empty_output><jupyter_text>Predict and plot assigned cell types<jupyter_code>predictions_hgsc = model_hgsc.predict()
sns.clustermap(predictions_hgsc, cmap="viridis")
hgsc_bdata.obs["scvi-tools predictions"] = predictions_hgsc.idxmax(axis=1).values
sc.pl.umap(
    hgsc_bdata,
    color=["celltype", "scvi-tools predictions"],
    ncols=1,
    frameon=False,
)<jupyter_output>... storing 'scvi-tools predictions' as categorical<jupyter_text>Model reproducibility<jupyter_code>df = hgsc_bdata.obs
confusion_matrix = pd.crosstab(
    df["scvi-tools predictions"],
    df["celltype"],
    rownames=["scvi-tools predictions"],
    colnames=["Original predictions"],
)
confusion_matrix /= confusion_matrix.sum(1).ravel().reshape(-1, 1)
fig, ax = plt.subplots(figsize=(5, 4))
sns.heatmap(
    confusion_matrix,
    cmap=sns.diverging_palette(245, 320, s=60, as_cmap=True),
    ax=ax,
    square=True,
    cbar_kws=dict(shrink=0.4, aspect=12),
)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1637661758143686, 'fraction_numerical': 0.028781793842034806, 'mean_word_length': 4.890275952693824, 'pattern_counts': {'":': 2, '<': 52, '<?xml version=': 0, '>': 52, 'https://': 8, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '198270', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3713, 'n_tokens_neox': 3129, 'n_words': 698, 'path': 'cellassign_tutorial.ipynb', 'repo_name': 'PierreBoyeau/scvi-tutorials', 'size': '8964'}