Data loading and preparation

Here we walk through the necessary steps to get your data into ready for scvi-tools.

Open In Colab

[1]:
import sys

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

if IN_COLAB and stable:
    !pip install --quiet scvi-tools[tutorials]
elif IN_COLAB and not stable:
    !pip install --quiet --upgrade jsonschema
    !pip install --quiet git+https://github.com/yoseflab/scvi-tools@master#egg=scvi-tools[tutorials]

     |████████████████████████████████| 61kB 2.1MB/s
ERROR: nbclient 0.5.0 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
     |████████████████████████████████| 112kB 3.5MB/s
     |████████████████████████████████| 122kB 11.1MB/s
     |████████████████████████████████| 163kB 11.3MB/s
     |████████████████████████████████| 8.7MB 5.2MB/s
     |████████████████████████████████| 3.2MB 46.9MB/s
     |████████████████████████████████| 51kB 6.2MB/s
     |████████████████████████████████| 2.4MB 41.3MB/s
     |████████████████████████████████| 7.7MB 33.3MB/s
     |████████████████████████████████| 51kB 5.9MB/s
     |████████████████████████████████| 112kB 51.8MB/s
     |████████████████████████████████| 61kB 7.3MB/s
     |████████████████████████████████| 51kB 6.1MB/s
  Building wheel for scvi-tools (PEP 517) ... done
  Building wheel for loompy (setup.py) ... done
  Building wheel for numpy-groupies (setup.py) ... done
  Building wheel for sinfo (setup.py) ... done
[2]:
import scvi
import scanpy as sc

Loading data

scvi-tools supports the AnnData data format, which also underlies Scanpy. AnnData is quite similar to other popular single cell objects like that of Seurat and SingleCellExperiment. In particular, it allows cell-level and feature-level metadata to coexist in the same data structure as the molecular counts.

It’s also now possible to automatically convert these R-based objects to AnnData within a Jupyter notebook. See the following tutorial for more information.

scvi-tools has a number of convenience methods for loading data from .csv, .loom, and .h5ad formats. To load ouputs from Cell Ranger, please use Scanpy’s reading functionality.

Let us now download an AnnData object (.h5ad format) and load it using scvi-tools.

PBMC3k

[3]:
!wget 'http://falexwolf.de/data/pbmc3k_raw.h5ad'
--2020-10-14 03:13:20--  http://falexwolf.de/data/pbmc3k_raw.h5ad
Resolving falexwolf.de (falexwolf.de)... 85.13.135.70
Connecting to falexwolf.de (falexwolf.de)|85.13.135.70|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 5855727 (5.6M)
Saving to: ‘pbmc3k_raw.h5ad’

pbmc3k_raw.h5ad     100%[===================>]   5.58M  4.94MB/s    in 1.1s

2020-10-14 03:13:22 (4.94 MB/s) - ‘pbmc3k_raw.h5ad’ saved [5855727/5855727]

[4]:
pbmc3k = scvi.data.read_h5ad("pbmc3k_raw.h5ad")
[5]:
pbmc3k
[5]:
AnnData object with n_obs × n_vars = 2700 × 32738
    var: 'gene_ids'

This is a fairly simple object, it just contains the count data and the ENSEMBL ids for the genes.

[6]:
pbmc3k.var.head()
[6]:
gene_ids
index
MIR1302-10 ENSG00000243485
FAM138A ENSG00000237613
OR4F5 ENSG00000186092
RP11-34P13.7 ENSG00000238009
RP11-34P13.8 ENSG00000239945

PBMC5k

As another example, let’s download a dataset from 10x Genomics. This data was obtained from a CITE-seq experiment, so it also contains protein count data.

[7]:
!wget https://cf.10xgenomics.com/samples/cell-exp/3.0.2/5k_pbmc_protein_v3/5k_pbmc_protein_v3_filtered_feature_bc_matrix.h5
--2020-10-14 03:13:22--  https://cf.10xgenomics.com/samples/cell-exp/3.0.2/5k_pbmc_protein_v3/5k_pbmc_protein_v3_filtered_feature_bc_matrix.h5
Resolving cf.10xgenomics.com (cf.10xgenomics.com)... 104.18.0.173, 104.18.1.173, 2606:4700::6812:ad, ...
Connecting to cf.10xgenomics.com (cf.10xgenomics.com)|104.18.0.173|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 17129253 (16M) [binary/octet-stream]
Saving to: ‘5k_pbmc_protein_v3_filtered_feature_bc_matrix.h5’

5k_pbmc_protein_v3_ 100%[===================>]  16.33M  21.3MB/s    in 0.8s

2020-10-14 03:13:24 (21.3 MB/s) - ‘5k_pbmc_protein_v3_filtered_feature_bc_matrix.h5’ saved [17129253/17129253]

[8]:
pbmc5k = sc.read_10x_h5(
    "5k_pbmc_protein_v3_filtered_feature_bc_matrix.h5",
    gex_only=False
)
Variable names are not unique. To make them unique, call `.var_names_make_unique`.

It’s often helpful to give the gene names unique names.

[9]:
pbmc5k.var_names_make_unique()

We can see that adata.X contains the concatenated gene and protein expression data.

[10]:
pbmc5k.var.feature_types.astype("category").cat.categories
[10]:
Index(['Antibody Capture', 'Gene Expression'], dtype='object')

We can use scvi-tools to organize this object, which places the protein expression in adata.obms["protein_expression].

[11]:
scvi.data.organize_cite_seq_10x(pbmc5k)
/usr/local/lib/python3.6/dist-packages/anndata/_core/anndata.py:1094: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead
  if not is_categorical(df_full[k]):
[12]:
pbmc5k
[12]:
AnnData object with n_obs × n_vars = 5247 × 33538
    var: 'gene_ids', 'feature_types', 'genome'
    obsm: 'protein_expression'

Concatenate the datasets

[13]:
adata = pbmc5k.concatenate(pbmc3k)

Notice that the resulting AnnData has a batch key in .obs.

[14]:
adata.obs.head()
[14]:
batch
AAACCCAAGAGACAAG-1-0 0
AAACCCAAGGCCTAGA-1-0 0
AAACCCAGTCGTGCCA-1-0 0
AAACCCATCGTGCATA-1-0 0
AAACGAAAGACAAGCC-1-0 0

Preprocessing the data

It is common to remove outliers, and even perform feature selection before model fitting. We prefer the Scanpy preprocessing module at this stage.

[15]:
sc.pp.filter_genes(adata, min_counts=3)
sc.pp.filter_cells(adata, min_counts=3)
/usr/local/lib/python3.6/dist-packages/anndata/_core/anndata.py:1094: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead
  if not is_categorical(df_full[k]):
/usr/local/lib/python3.6/dist-packages/anndata/_core/anndata.py:1094: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead
  if not is_categorical(df_full[k]):

As it is popular to use normalize the data for many methods, we can use Scanpy for this; however, it’s important to keep the count information intact for scvi-tools models.

[16]:
adata.layers["counts"] = adata.X.copy()

Now we can proceed with common normalization methods.

[17]:
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

We can store the normalized values in .raw to keep them safe in the event the anndata gets subsetted feature-wise.

[18]:
adata.raw = adata

Register the data with scvi-tools

Now that we have an AnnData object, we need to alert scvi-tools of all the interesting data in our object. For example, now that we have batches in our AnnData, we can alert the models that we’d like to perform batch correction. Also, because we have the count data in a layer, we can use the layer argument.

Basic case

[19]:
scvi.data.setup_anndata(adata, layer="counts", batch_key="batch")
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      Successfully registered anndata object containing 7947 cells, 14309
          vars, 2 batches, 1 labels, and 0 proteins. Also registered 0 extra
          categorical covariates and 0 extra continuous covariates.
INFO      Please do not further modify adata until model is trained.
/usr/local/lib/python3.6/dist-packages/anndata/_core/anndata.py:1094: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead
  if not is_categorical(df_full[k]):
/usr/local/lib/python3.6/dist-packages/anndata/_core/anndata.py:1094: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead
  if not is_categorical(df_full[k]):

Notice the info messages notify us that batches were detected in the data. Just to demonstrate what happens if we don’t include this option:

[20]:
scvi.data.setup_anndata(adata, layer="counts")
INFO      No batch_key inputted, assuming all cells are same 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      Successfully registered anndata object containing 7947 cells, 14309
          vars, 1 batches, 1 labels, and 0 proteins. Also registered 0 extra
          categorical covariates and 0 extra continuous covariates.
/usr/local/lib/python3.6/dist-packages/anndata/_core/anndata.py:1094: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead
  if not is_categorical(df_full[k]):
INFO      Please do not further modify adata until model is trained.

Now integration-based tasks can no longer be performed in subsequent models because there is no knowledge of such information.

CITE-seq case

As PBMC5k is a CITE-seq dataset, we can use scvi-tools to register the protein expression. Note that totalVI is the only current model that uses the protein expression. The usage of registered items is model specific. As another example, registering the labels in the AnnData object will not affect totalVI or scVI, but is necessary to run scANVI.

We have not preprocessed the pbmc5k object, which we do recommend. We show how to run setup_anndata in this case for illustrative purposes.

[21]:
scvi.data.setup_anndata(pbmc5k, protein_expression_obsm_key="protein_expression")
INFO      No batch_key inputted, assuming all cells are same batch
INFO      No label_key inputted, assuming all cells have same label
INFO      Using data from adata.X
INFO      Computing library size prior per batch
INFO      Using protein expression from adata.obsm['protein_expression']
/usr/local/lib/python3.6/dist-packages/anndata/_core/anndata.py:1094: FutureWarning: is_categorical is deprecated and will be removed in a future version.  Use is_categorical_dtype instead
  if not is_categorical(df_full[k]):
INFO      Using protein names from columns of adata.obsm['protein_expression']
INFO      Successfully registered anndata object containing 5247 cells, 33538
          vars, 1 batches, 1 labels, and 32 proteins. Also registered 0 extra
          categorical covariates and 0 extra continuous covariates.
INFO      Please do not further modify adata until model is trained.

Warning

After setup_anndata has been run, the adata object should not be modified. In other words, the very next step in the workflow is to initialize and train the model of interest (e.g., scVI, totalVI). If you do modify the adata, it’s ok, just run setup_anndata again – and then reinitialize the model.

Viewing the scvi-tools data setup

[22]:
scvi.data.view_anndata_setup(pbmc5k)
Anndata setup with scvi-tools version 0.7.0b0.
              Data Summary              
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━┓
┃             Data              Count ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━┩
│            Cells              5247  │
│             Vars              33538 │
│            Labels               1   │
│           Batches               1   │
│           Proteins             32   │
│ Extra Categorical Covariates    0   │
│ Extra Continuous Covariates     0   │
└──────────────────────────────┴───────┘
                   SCVI Data Registry                    
┏━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃        Data               scvi-tools Location        ┃
┡━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│         X                       adata.X              │
│   batch_indices         adata.obs['_scvi_batch']     │
│    local_l_mean     adata.obs['_scvi_local_l_mean']  │
│    local_l_var       adata.obs['_scvi_local_l_var']  │
│       labels           adata.obs['_scvi_labels']     │
│ protein_expression  adata.obsm['protein_expression'] │
└────────────────────┴──────────────────────────────────┘
                        Label Categories                        
┏━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃      Source Location       Categories  scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['_scvi_labels']      0                0          │
└───────────────────────────┴────────────┴─────────────────────┘
                       Batch Categories                        
┏━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━┓
┃     Source Location       Categories  scvi-tools Encoding ┃
┡━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━┩
│ adata.obs['_scvi_batch']      0                0          │
└──────────────────────────┴────────────┴─────────────────────┘