Chapter 3: Quality Control, Mitochondrial Outliers & Doublet Scrubbing

Chapter 3: Quality Control, Mitochondrial Outliers & Doublet Scrubbing

Single-cell datasets inevitably contain technical noise: apoptotic/damaged cells with ruptured membranes, empty droplets containing ambient background RNA, and droplets containing two or more encapsulated cells (multiplets). Stringent quality control is the single most critical determinant of valid downstream biological discoveries.

1. Quality Control with Seurat v5 in R

library(Seurat)
library(ggplot2)

# Load 10x matrix data
counts <- Read10X(data.dir = "filtered_feature_bc_matrix/")
pbmc <- CreateSeuratObject(counts = counts, project = "PBMC10K", min.cells = 3, min.features = 200)

# Calculate percentage of mitochondrial and ribosomal reads
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc[["percent.rb"]] <- PercentageFeatureSet(pbmc, pattern = "^RP[SL]")

# Visualize QC metrics with violin plots
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3, pt.size = 0.1)

# Apply adaptive or strict biological thresholds
# Filtering: Keep cells with >500 & <6000 genes, and <15% mitochondrial reads
pbmc_clean <- subset(pbmc, subset = nFeature_RNA > 500 & nFeature_RNA < 6000 & percent.mt < 15)
message("Retained ", ncol(pbmc_clean), " cells out of ", ncol(pbmc))

2. Quality Control with Scanpy in Python

import scanpy as sc
import numpy as np

# Read 10x HDF5 or matrix directory into an AnnData object
adata = sc.read_10x_h5("sample_filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()

# Annotate mitochondrial, ribosomal, and hemoglobin genes
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['ribo'] = adata.var_names.str.startswith(('RPS', 'RPL'))
adata.var['hb'] = adata.var_names.str.startswith(('HBA', 'HBB'))

# Calculate comprehensive metrics
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt', 'ribo', 'hb'], percent_top=None, log1p=False, inplace=True)

# Adaptive filtering using Median Absolute Deviations (MAD)
def is_outlier(adata, metric: str, nmads: int):
    M = adata.obs[metric]
    return (M < np.median(M) - nmads * np.median(np.abs(M - np.median(M)))) | \
           (M > np.median(M) + nmads * np.median(np.abs(M - np.median(M))))

adata = adata[~is_outlier(adata, 'pct_counts_mt', 3) & (adata.obs.pct_counts_mt < 15)]
adata = adata[(adata.obs.n_genes_by_counts > 400) & (adata.obs.n_genes_by_counts < 6500)]

3. Computational Doublet Scrubbing with Scrublet & DoubletFinder

Homotypic doublets (two cells of the same type) alter gene expression magnitude, but heterotypic doublets (e.g. a B-cell and a T-cell in the same droplet) create artificial "hybrid" intermediate clusters that can be mistaken for novel transitional cell types.

# Doublet detection using Scrublet in Python
import scrublet as scr

# Initialize Scrublet with raw unnormalized count matrix
scrub = scr.Scrublet(adata.X, expected_doublet_rate=0.06)
doublet_scores, predicted_doublets = scrub.scrub_doublets(min_counts=2, min_cells=3, min_gene_variability_pctl=85, n_prin_comps=30)

adata.obs['doublet_score'] = doublet_scores
adata.obs['is_doublet'] = predicted_doublets

# Filter out predicted doublets
adata_singlets = adata[~adata.obs['is_doublet']].copy()
print(f"Removed {sum(predicted_doublets)} heterotypic doublets.")