Marker Genes

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

Find markers for each cluster

sc.tl.rank_genes_groups(adata, groupby="leiden", method="wilcoxon")

result_df = sc.get.rank_genes_groups_df(adata, group="0")
print(result_df.head(8).to_string(index=False))

names    scores  logfoldchanges         pvals     pvals_adj
 LDHB 31.944715        2.778137 6.397028e-224 8.619356e-220
RPS12 30.740160        0.983159 1.655289e-207 7.434454e-204
RPS25 29.029331        1.140002 2.806448e-185 7.562817e-182
RPS27 28.159712        0.985933 1.822450e-174 3.069462e-171
 RPS6 27.788717        0.823007 5.937959e-170 8.889784e-167
 RPS3 27.520372        0.834473 1.001657e-166 1.349632e-163
 TPT1 27.140507        0.925139 3.277350e-162 3.396847e-159
 CD3D 26.998106        3.200113 1.555552e-160 1.497108e-157

result = adata.uns["rank_genes_groups"]
groups = result["names"].dtype.names

for group in groups:
    top = result["names"][group][:3]
    print(f"{group}: {', '.join(top)}")

0: LDHB, RPS12, RPS25
1: LYZ, S100A9, CST3
2: NKG7, CST7, CCL5
3: CD74, CD79A, HLA-DRA
4: LST1, FCER1G, FCGR3A

sc.pl.rank_genes_groups_dotplot(adata, n_genes=4, groupby="leiden",
                                show=False)
plt.savefig("outputs/marker-genes_dotplot.png", dpi=150,
            bbox_inches="tight")
plt.close()

Dot plot of the top four marker genes per cluster. Every cluster carries its own strong markers and near-zero expression elsewhere, which is what makes the labels stick.

sc.pl.rank_genes_groups_heatmap(adata, n_genes=5, show_gene_labels=False,
                                show=False)
plt.savefig("outputs/marker-genes_heatmap.png", dpi=150,
            bbox_inches="tight")
plt.close()

Heatmap of the top five marker genes per cluster. Blocks along the diagonal show each cluster expressing its own marker set and little else.