Preprocessing and clustering 3k PBMCs ¶

This guide demonstrates how to use lamindb to track data lineage and manage artifacts during a standard single-cell RNA-seq analysis workflow. Adapted from the Scanpy PBMC3k tutorial, it has additional code showing how to track the notebook’s execution and how to save the raw and the the final annotated dataset in your database.

We will process a dataset of 3k PBMCs from a Healthy Donor (freely available from 10x Genomics).

Note: To run examples, if you don’t have a LaminDB instance, create one:

lamin init --storage ./test-pbmc3k

Import modules:

import matplotlib.pyplot as plt
import pandas as pd
import scanpy as sc
import lamindb as ln

Track the notebook’s execution:

ln.track()

Download and unpack the data:

from pathlib import Path
import tarfile
from urllib.request import Request, urlopen

data_dir = Path("data")
data_dir.mkdir(exist_ok=True)
archive_path = data_dir / "pbmc3k_filtered_gene_bc_matrices.tar.gz"
if not archive_path.exists():
    urls = [
        "https://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz",
        "https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz",
    ]
    for i, url in enumerate(urls):
        try:
            req = Request(url, headers={"User-Agent": "Mozilla/5.0"})
            with urlopen(req) as response, archive_path.open("wb") as f:
                f.write(response.read())
            break
        except Exception:
            if i == len(urls) - 1:
                raise
with tarfile.open(archive_path, "r:gz") as tar:
    tar.extractall(data_dir)

Read in the count matrix into an AnnData object:

adata = sc.read_10x_mtx(
    "data/filtered_gene_bc_matrices/hg19/",  # the directory with the `.mtx` file
    var_names="gene_symbols",  # use gene symbols for the variable names (variables-axis index)
    cache=True,  # write a cache file for faster subsequent reading
)
adata

Save the AnnData object with the raw data:

ln.Artifact.from_anndata(adata, key="raw_pbmc3k.h5ad").save()

Preprocessing¶

Load the raw dataset:

adata = ln.Artifact.get(key="raw_pbmc3k.h5ad").load()

Show those genes that yield the highest fraction of counts in each single cell, across all cells:

sc.pl.highest_expr_genes(adata, n_top=20)

Basic filtering:

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata

Let’s assemble some information about mitochondrial genes, which are important for quality control:

# annotate the group of mitochondrial genes as "mt"
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True)

A violin plot of some of the computed quality measures:

the number of genes expressed in the count matrix
the total counts per cell
the percentage of counts in mitochondrial genes

sc.pl.violin(
    adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    multi_panel=True,
)

Remove cells that have too many mitochondrial genes expressed or too many total counts. Do the filtering by slicing the AnnData object:

adata = adata[
    (adata.obs.n_genes_by_counts < 2500) & (adata.obs.n_genes_by_counts > 200) & (adata.obs.pct_counts_mt < 5),
    :,
].copy()
adata.layers["counts"] = adata.X.copy()
adata

Total-count normalize (library-size correct) the data matrix $\mathbf{X}$ to 10,000 reads per cell, so that counts become comparable among cells:

sc.pp.normalize_total(adata, target_sum=1e4)

Logarithmize the data:

sc.pp.log1p(adata)

Identify highly-variable genes:

sc.pp.highly_variable_genes(
    adata,
    layer="counts",
    n_top_genes=2000,
    min_mean=0.0125,
    max_mean=3,
    min_disp=0.5,
    flavor="seurat_v3",
)

Scale each gene to unit variance. Clip values exceeding standard deviation 10.

adata.layers["scaled"] = adata.X.toarray(order='C')
sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"], layer="scaled")
sc.pp.scale(adata, max_value=10, layer="scaled")

Principal component analysis¶

Reduce the dimensionality of the data by running principal component analysis (PCA), which reveals the main axes of variation and denoises the data.

sc.pp.pca(adata, layer="scaled", svd_solver="arpack")

Let us inspect the contribution of single PCs to the total variance in the data. This gives us information about how many PCs we should consider in order to compute the neighborhood relations of cells, e.g. used in the clustering function sc.tl.louvain() or tSNE sc.tl.tsne().

sc.pl.pca_variance_ratio(adata, n_pcs=20)

Computing the neighborhood graph¶

Compute the neighborhood graph of cells using the PCA representation of the data matrix.

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

Clustering the neighborhood graph¶

sc.tl.leiden(
    adata,
    resolution=0.7,
    random_state=0,
    flavor="igraph",
    n_iterations=2,
    directed=False,
)
adata.obs["leiden"] = adata.obs["leiden"].copy()
adata.uns["leiden"] = adata.uns["leiden"].copy()

Embedding the neighborhood graph in a UMAP¶

sc.tl.paga(adata)
sc.pl.paga(adata, plot=False)  # remove `plot=False` if you want to see the coarse-grained graph
sc.tl.umap(adata, init_pos='paga')
sc.tl.umap(adata)

adata.obsm["X_umap"] = adata.obsm["X_umap"].copy()
sc.pl.umap(adata, color=["leiden", "CD14", "NKG7"])

Finding marker genes¶

Let us compute a ranking for the highly differential genes in each cluster.

sc.tl.rank_genes_groups(adata, "leiden", mask_var="highly_variable", method="wilcoxon")
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

Most of these genes are found in the highly expressed genes with notable exceptions like CD8A (see the below dotplot).

Leiden Group	Markers	Cell Type
0	IL7R	CD4 T cells
1	CD14, LYZ	CD14+ Monocytes
2	MS4A1	B cells
3	CD8A	CD8 T cells
4	GNLY, NKG7	NK cells
5	FCGR3A, MS4A7	FCGR3A+ Monocytes
6	FCER1A, CST3	Dendritic Cells
7	PPBP	Megakaryocytes

Let us also define a list of marker genes for later reference.

marker_genes = [
    *["IL7R", "CD79A", "MS4A1", "CD8A", "CD8B", "LYZ", "CD14"],
    *["LGALS3", "S100A8", "GNLY", "NKG7", "KLRB1"],
    *["FCGR3A", "MS4A7", "FCER1A", "CST3", "PPBP"],
]

Actually mark the cell types.

new_cluster_names = [
    "CD4 T",
    "B",
    "CD14+ Monocytes",
    "NK",
    "CD8 T",
    "FCGR3A+ Monocytes",
    "Dendritic",
    "Megakaryocytes",
]
adata.rename_categories("leiden", new_cluster_names)
adata.obs['cell_type']=adata.obs['leiden']
adata

Visualize the cell types in a UMAP

sc.pl.umap(adata, color="cell_type", legend_loc="on data", title="", frameon=False)

Now that we annotated the cell types, let us visualize the marker genes.

sc.pl.dotplot(adata, marker_genes, groupby="cell_type")

Save the resulting AnnData object in LaminDB:

ln.Artifact.from_anndata(adata, key="result_pbmc3k.h5ad").save()

ln.finish() marks the successful completion of the current script or notebook run in LaminDB, finalizing its execution state and saving the source code to the database for reproducibility.

ln.finish()