Annotation

Load, filter, embed, cluster

import numpy as np
import scanpy as sc
import matplotlib.pyplot as plt

sc.settings.verbosity = 3
np.random.seed(42)

adata = sc.read_10x_mtx("/opt/data/pbmc3k", var_names="gene_symbols",
                        make_unique=True)
adata.layers["counts"] = adata.X.copy()
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)


def is_outlier(adata, metric, nmads):
    """Flag cells whose metric sits more than nmads MADs from the median.

    Args:
        adata: AnnData object holding the metric in .obs.
        metric: Column of adata.obs to test.
        nmads: Number of median absolute deviations for the cutoff.

    Returns:
        Boolean Series, True where the cell is an outlier.
    """
    values = adata.obs[metric]
    median = values.median()
    mad = (values - median).abs().median()
    return (values < median - nmads * mad) | (values > median + nmads * mad)


adata.obs["log1p_total_counts"] = np.log1p(adata.obs["total_counts"])
adata.obs["log1p_n_genes_by_counts"] = np.log1p(
    adata.obs["n_genes_by_counts"]
)
adata.obs["outlier"] = (
    is_outlier(adata, "log1p_total_counts", 5)
    | is_outlier(adata, "log1p_n_genes_by_counts", 5)
    | is_outlier(adata, "pct_counts_mt", 3)
)
adata = adata[~adata.obs["outlier"]].copy()
sc.pp.filter_genes(adata, min_cells=3)

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, layer="counts", n_top_genes=2000,
                            flavor="seurat_v3")
adata.raw = adata
adata = adata[:, adata.var["highly_variable"]]
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50, random_state=42)
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30, random_state=42)
sc.tl.umap(adata, random_state=42)
sc.tl.leiden(adata, resolution=0.5, key_added="leiden", random_state=42,
             flavor="igraph", n_iterations=2, directed=False)
print(adata.obs["leiden"].value_counts().sort_index())

leiden
0    1099
1     456
2     412
3     315
4     146
Name: count, dtype: int64

Marker genes

markers = {
    "T cells": ["CD3D", "CD3E"],
    "B cells": ["MS4A1", "CD79A"],
    "Monocytes": ["CD14", "LYZ"],
    "NK cells": ["GNLY", "NKG7"],
    "DC": ["FCER1A", "CST3"],
}

sc.pl.umap(adata, color=["CD3D", "MS4A1", "CD14", "GNLY"], use_raw=True,
           show=False)
plt.savefig("outputs/annotation_marker-umap.png", dpi=150,
            bbox_inches="tight")
plt.close()

UMAP colored by four marker genes. CD3D marks the T cell arm, MS4A1 a tight B cell cluster, CD14 the monocyte region and GNLY the NK cells.

sc.pl.dotplot(adata, markers, groupby="leiden", use_raw=True, show=False)
plt.savefig("outputs/annotation_dotplot.png", dpi=150,
            bbox_inches="tight")
plt.close()

Dot plot of marker genes per Leiden cluster. Each cluster shows one strong marker row, and the DC markers stay weak everywhere except the smallest cluster.

Marker scores

# score_genes averages a marker set and subtracts a matched control
# set, so each cluster gets one number per candidate identity.
marker_sets = {
    "T cells": ["CD3D", "CD3E", "IL7R"],
    "B cells": ["MS4A1", "CD79A", "CD79B"],
    "Monocytes": ["CD14", "LYZ", "S100A9"],
    "NK cells": ["GNLY", "NKG7", "KLRD1"],
    "DC": ["FCER1A", "CST3", "CLEC10A"],
}

for name, genes in marker_sets.items():
    sc.tl.score_genes(adata, gene_list=genes, score_name=f"score_{name}",
                      use_raw=True)
    scores = adata.obs.groupby("leiden", observed=True)[
        f"score_{name}"
    ].mean()
    print(f"{name:10s} top cluster {scores.idxmax()} ({scores.max():.3f})")

T cells    top cluster 0 (0.555)
B cells    top cluster 3 (1.891)
Monocytes  top cluster 1 (2.770)
NK cells   top cluster 2 (1.572)
DC         top cluster 1 (0.808)

From clusters to cell types

# Cluster 2 carries CD3D together with CD8A and GNLY, so it holds
# both CD8 T cells and NK cells. The resolution sweep on the
# clustering page splits it once the resolution rises above 0.5.
cluster_labels = {
    "0": "CD4 T cells",
    "1": "CD14+ Monocytes",
    "2": "CD8 T and NK cells",
    "3": "B cells",
    "4": "FCGR3A+ Monocytes",
}

adata.obs["cell_type"] = adata.obs["leiden"].map(cluster_labels)
print(adata.obs["cell_type"].value_counts())

cell_type
CD4 T cells           1099
CD14+ Monocytes        456
CD8 T and NK cells     412
B cells                315
FCGR3A+ Monocytes      146
Name: count, dtype: int64

sc.pl.umap(adata, color="cell_type", show=False)
plt.savefig("outputs/annotation_cell-types.png", dpi=150,
            bbox_inches="tight")
plt.close()

UMAP colored by the final cell type labels. The CD4 T cells hold the largest region, the monocytes split into two neighboring groups, and the merged CD8 T and NK population holds its own arm.