Feature Selection

Load and normalize

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)

# The seurat_v3 and pearson_residuals flavours need raw counts, the
# seurat and cell_ranger flavours need log-normalized data.
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
print(adata.shape)

(2428, 13474)

Highly variable genes

# Find 2,000 highly variable genes on the raw counts.
sc.pp.highly_variable_genes(adata, layer="counts", n_top_genes=2000,
                            flavor="seurat_v3")
print(adata.var["highly_variable"].sum())

sc.pl.highly_variable_genes(adata, show=False)
plt.savefig("outputs/feature-selection_hvg.png", dpi=150,
            bbox_inches="tight")
plt.close()

Highly variable genes on the mean-variance plane. The 2,000 selected genes sit on the upper edge of the cloud, where variance exceeds what the mean alone predicts.

Comparing the flavours

# sc-best-practices recommends testing more than one flavour. The
# seurat and cell_ranger flavours run on the log-normalized data, the
# pearson_residuals flavour on the raw counts.
flavors = {}
for flavor in ("seurat", "cell_ranger"):
    sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor=flavor)
    flavors[flavor] = set(adata.var_names[adata.var["highly_variable"]])

sc.experimental.pp.highly_variable_genes(adata, layer="counts",
                                          n_top_genes=2000,
                                          flavor="pearson_residuals")
flavors["pearson_residuals"] = set(
    adata.var_names[adata.var["highly_variable"]]
)

sc.pp.highly_variable_genes(adata, layer="counts", n_top_genes=2000,
                            flavor="seurat_v3")
flavors["seurat_v3"] = set(adata.var_names[adata.var["highly_variable"]])

for name, genes in flavors.items():
    print(f"{name:20s} {len(genes)} genes")
print("seurat and seurat_v3 share:",
      len(flavors["seurat"] & flavors["seurat_v3"]))
print("seurat_v3 and pearson_residuals share:",
      len(flavors["seurat_v3"] & flavors["pearson_residuals"]))

seurat               2000 genes
cell_ranger          2000 genes
pearson_residuals    2000 genes
seurat_v3            2000 genes
seurat and seurat_v3 share: 1576
seurat_v3 and pearson_residuals share: 1460

Keeping the full data

# Keep the full dataset in .raw, then restrict the object to the
# highly variable genes for the downstream steps.
adata.raw = adata
adata = adata[:, adata.var["highly_variable"]]
print(adata.shape)

(2428, 2000)