单细胞 RNA-seq 质量控制：快速入门

20 分钟完成你的第一次单细胞质控

本章通过一个完整的实战案例，带你快速掌握单细胞 QC 的核心流程。我们将分析一个 PBMC（外周血单核细胞）数据集，从原始数据到高质量过滤结果。

实战案例：PBMC 数据质控

数据背景

数据类型：PBMC 单细胞 RNA-seq（10X Genomics v3）
细胞数量：5,000 个细胞
物种：人类
格式：.h5ad 文件

分析目标

计算各项 QC 指标
识别低质量细胞和异常值
应用基于 MAD 的过滤
验证过滤结果

步骤 1：环境准备与数据加载 [3 分钟]

1.1 确保环境就绪

# 检查依赖包
import scanpy as sc
import anndata as ad
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

print(f"scanpy version: {sc.__version__}")
print(f"anndata version: {ad.__version__}")

# 设置 Scanpy 参数
sc.settings.verbosity = 3  # 显示详细信息
sc.settings.set_figure_params(dpi=80, facecolor='white')

1.2 加载数据

# 加载 .h5ad 文件
adata = sc.read_h5ad('pbmc_raw.h5ad')

# 查看数据基本信息
print(f"数据形状：{adata.n_obs} 细胞 × {adata.n_vars} 基因")
print(f"数据类型：{adata.X.dtype}")
print(f"稀疏度：{(adata.X == 0).sum() / adata.X.size * 100:.1f}%")

# 查看 AnnData 结构
print("\nAnnData 对象结构：")
print(adata)

预期输出示例：

数据形状: 5000 细胞 × 32738 基因
数据类型: float64
稀疏度: 95.2%

AnnData 对象结构:
AnnData object with n_obs × n_vars = 5000 × 32738

⚠️ 重要检查：

确保 adata.X 是原始 UMI 计数（未归一化）
如果数据已经 log-normalized，需要重新加载原始数据

步骤 2：计算 QC 指标 [5 分钟]

2.1 识别线粒体、核糖体基因

# 人类线粒体基因以 'MT-' 开头
adata.var['mt'] = adata.var_names.str.startswith('MT-')

# 人类核糖体基因（常用方法）
ribosome_genes = [
    'RPS', 'RPL',  # 核糖体蛋白
    'MRPS', 'MRPL'  # 线粒体核糖体蛋白
]
adata.var['ribo'] = adata.var_names.str.startswith(tuple(ribosome_genes))

# 血红蛋白基因（针对血液样本）
adata.var['hb'] = adata.var_names.str.startswith(('HB', 'Hb'))

# 打印统计
print(f"线粒体基因数量：{adata.var['mt'].sum()}")
print(f"核糖体基因数量：{adata.var['ribo'].sum()}")
print(f"血红蛋白基因数量：{adata.var['hb'].sum()}")

物种差异：

# 小鼠：线粒体基因以 'mt-' 开头（小写）
adata.var['mt'] = adata.var_names.str.startswith('mt-')

# 其他物种：查阅对应的前缀规则
# 例如，斑马鱼：'mt-' 或 'MT-'

2.2 计算每个细胞的 QC 指标

# 计算 QC 指标
sc.pp.calculate_qc_metrics(
    adata,
    qc_vars=['mt', 'ribo', 'hb'],
    percent_top=None,  # 不计算 top 基因比例
    log1p=False,       # 不进行对数变换
    inplace=True
)

# 查看新增的 obs 列
print("\n新增 QC 指标：")
qc_cols = [col for col in adata.obs.columns if col.startswith('pct_') or col.startswith('n_')]
print(qc_cols)

新增指标说明：

n_genes_by_counts: 检测到的基因数
n_counts: 总 UMI 数
pct_counts_mt: 线粒体基因 UMI 百分比
pct_counts_ribo: 核糖体基因 UMI 百分比
pct_counts_hb: 血红蛋白基因 UMI 百分比

2.3 查看 QC 指标统计

# 查看基本统计
print("=== QC 指标统计 ===")
print(adata.obs[['n_genes_by_counts', 'n_counts', 'pct_counts_mt']].describe())

# 查看前几个细胞
print("\n前 5 个细胞的 QC 指标：")
print(adata.obs[['n_genes_by_counts', 'n_counts', 'pct_counts_mt']].head())

预期输出（示例）：

=== QC 指标统计 ===
       n_genes_by_counts      n_counts  pct_counts_mt
count        5000.000000  5.000000e+03    5000.000000
mean         1823.456000  6.234567e+03       5.234567
std           567.890123  2.345678e+03       3.456789
min           200.000000  5.000000e+02       0.500000
25%          1456.000000  4.567890e+03       3.000000
50%          1789.000000  5.987654e+03       4.800000
75%          2156.000000  7.654321e+03       6.800000
max          4567.000000  2.345678e+04      25.400000

步骤 3：可视化 QC 分布 [5 分钟]

3.1 小提琴图（Violin Plot）

# 创建小提琴图
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# n_genes_by_counts
sc.pl.violin(
    adata,
    ['n_genes_by_counts'],
    groupby=None,
    jitter=0.4,
    multi_panel=False,
    ax=axes[0],
    show=False
)
axes[0].set_title('检测到的基因数')

# n_counts
sc.pl.violin(
    adata,
    ['n_counts'],
    groupby=None,
    jitter=0.4,
    multi_panel=False,
    ax=axes[1],
    show=False
)
axes[1].set_title('总 UMI 数')

# pct_counts_mt
sc.pl.violin(
    adata,
    ['pct_counts_mt'],
    groupby=None,
    jitter=0.4,
    multi_panel=False,
    ax=axes[2],
    show=False
)
axes[2].set_title('线粒体基因比例')

plt.tight_layout()
plt.show()

解读要点：

理想分布：单峰，略右偏（多数细胞聚集在中位数附近）
警惕信号：双峰分布（可能混入低质量细胞群）
关注离群值：右侧极端值可能是双细胞

3.2 散点图（Scatter Plot）

# n_genes vs n_counts
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 左图：按线粒体比例着色
sc.pl.scatter(
    adata,
    x='n_counts',
    y='n_genes_by_counts',
    color='pct_counts_mt',
    ax=axes[0],
    show=False
)

# 右图：简单的 matplotlib 散点图（更灵活）
axes[1].scatter(
    adata.obs['n_counts'],
    adata.obs['n_genes_by_counts'],
    c=adata.obs['pct_counts_mt'],
    cmap='viridis',
    alpha=0.5,
    s=10
)
axes[1].set_xlabel('Total UMI counts')
axes[1].set_ylabel('Number of genes')
axes[1].set_title('Gene counts vs UMI counts')
plt.colorbar(axes[1].collections[0], ax=axes[1], label='MT %')

plt.tight_layout()
plt.show()

解读要点：

正相关：n_genes 和 n_counts 应该呈正相关
异常点：
- 右下角（高 UMI，低基因数）：可能是双细胞或粘连体
- 左下角（低 UMI，低基因数）：低质量/破碎细胞
- 高 MT%（黄色点）：死亡/应激细胞

步骤 4：基于 MAD 的异常值检测 [5 分钟]

4.1 定义 MAD 过滤函数

def mad_outliers(adata, metric, n_mads=3, direction='both'):
    """
    基于 MAD 的异常值检测

    参数：
        adata: AnnData 对象
        metric: 要检查的指标名（如 'n_genes_by_counts'）
        n_mads: MAD 倍数（推荐 3-5）
        direction: 'both', 'lower', 或 'upper'

    返回：
        布尔 Series（True 表示是异常值）
    """
    data = adata.obs[metric]
    median = np.median(data)
    mad = np.median(np.abs(data - median))

    # 计算阈值
    if direction in ['both', 'lower']:
        lower_threshold = median - n_mads * mad
    else:
        lower_threshold = -np.inf

    if direction in ['both', 'upper']:
        upper_threshold = median + n_mads * mad
    else:
        upper_threshold = np.inf

    # 识别异常值
    outliers = (data < lower_threshold) | (data > upper_threshold)

    # 打印信息
    n_outliers = outliers.sum()
    print(f"{metric}:")
    print(f"  中位数 = {median:.1f}")
    print(f"  MAD = {mad:.1f}")
    print(f"  阈值 = [{lower_threshold:.1f}, {upper_threshold:.1f}]")
    print(f"  异常值数量 = {n_outliers} ({n_outliers/len(data)*100:.1f}%)\n")

    return outliers

4.2 应用 MAD 过滤

# 设置 MAD 倍数（推荐 3）
N_MADS = 3

# 检测 n_genes_by_counts 的异常值（双向）
outliers_genes_low = mad_outliers(adata, 'n_genes_by_counts', n_mads=N_MADS, direction='lower')
outliers_genes_high = mad_outliers(adata, 'n_genes_by_counts', n_mads=N_MADS, direction='upper')

# 检测 n_counts 的异常值（双向）
outliers_counts_low = mad_outliers(adata, 'n_counts', n_mads=N_MADS, direction='lower')
outliers_counts_high = mad_outliers(adata, 'n_counts', n_mads=N_MADS, direction='upper')

# 检测 pct_counts_mt 的异常值（仅高值）
outliers_mt = mad_outliers(adata, 'pct_counts_mt', n_mads=N_MADS, direction='upper')

预期输出（示例）：

n_genes_by_counts:
  中位数 = 1789.0
  MAD = 285.6
  阈值 = [932.2, 2645.8]
  异常值数量 = 145 (2.9%)

n_counts:
  中位数 = 5987.5
  MAD = 1234.2
  阈值 = [2284.9, 9690.1]
  异常值数量 = 156 (3.1%)

pct_counts_mt:
  中位数 = 4.8
  MAD = 1.9
  阈值 = [-inf, 10.5]
  异常值数量 = 87 (1.7%)

4.3 可视化过滤阈值

# 创建直方图并标记阈值
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# n_genes_by_counts
axes[0].hist(adata.obs['n_genes_by_counts'], bins=50, alpha=0.7)
axes[0].axvline(
    np.median(adata.obs['n_genes_by_counts']) - N_MADS * np.median(np.abs(adata.obs['n_genes_by_counts'] - np.median(adata.obs['n_genes_by_counts']))),
    color='red', linestyle='--', label='Lower threshold'
)
axes[0].axvline(
    np.median(adata.obs['n_genes_by_counts']) + N_MADS * np.median(np.abs(adata.obs['n_genes_by_counts'] - np.median(adata.obs['n_genes_by_counts']))),
    color='red', linestyle='--', label='Upper threshold'
)
axes[0].set_xlabel('Number of genes')
axes[0].set_title('Gene counts distribution')
axes[0].legend()

# n_counts
axes[1].hist(adata.obs['n_counts'], bins=50, alpha=0.7, color='green')
axes[1].axvline(
    np.median(adata.obs['n_counts']) - N_MADS * np.median(np.abs(adata.obs['n_counts'] - np.median(adata.obs['n_counts']))),
    color='red', linestyle='--', label='Lower threshold'
)
axes[1].axvline(
    np.median(adata.obs['n_counts']) + N_MADS * np.median(np.abs(adata.obs['n_counts'] - np.median(adata.obs['n_counts']))),
    color='red', linestyle='--', label='Upper threshold'
)
axes[1].set_xlabel('Total UMI counts')
axes[1].set_title('UMI counts distribution')
axes[1].legend()

# pct_counts_mt
axes[2].hist(adata.obs['pct_counts_mt'], bins=50, alpha=0.7, color='orange')
threshold_mt = np.median(adata.obs['pct_counts_mt']) + N_MADS * np.median(np.abs(adata.obs['pct_counts_mt'] - np.median(adata.obs['pct_counts_mt'])))
axes[2].axvline(threshold_mt, color='red', linestyle='--', label='Threshold')
axes[2].set_xlabel('Mitochondrial %')
axes[2].set_title('Mitochondrial percentage distribution')
axes[2].legend()

plt.tight_layout()
plt.show()

步骤 5：应用过滤并验证 [5 分钟]

5.1 创建过滤 mask

# 综合所有过滤条件
filter_mask = (
    ~outliers_genes_low &      # 基因数不能太少
    ~outliers_genes_high &     # 基因数不能太多（双细胞）
    ~outliers_counts_low &     # UMI 数不能太少
    ~outliers_counts_high &    # UMI 数不能太多（双细胞）
    ~outliers_mt               # 线粒体比例不能太高
)

# 打印统计
print(f"过滤前细胞数：{adata.n_obs}")
print(f"过滤后细胞数：{filter_mask.sum()}")
print(f"保留率：{filter_mask.sum() / adata.n_obs * 100:.1f}%")

5.2 应用过滤

# 创建过滤后的数据副本
adata_filtered = adata[filter_mask].copy()

print(f"\n过滤后数据形状：{adata_filtered.n_obs} 细胞 × {adata_filtered.n_vars} 基因")

5.3 过滤低表达基因

# 过滤在少于 10 个细胞中表达的基因
print(f"过滤前基因数：{adata_filtered.n_vars}")

sc.pp.filter_genes(adata_filtered, min_cells=10)

print(f"过滤后基因数：{adata_filtered.n_vars}")
print(f"过滤掉的基因数：{adata.n_vars - adata_filtered.n_vars}")

5.4 验证过滤结果

# 重新计算 QC 指标并可视化
sc.pp.calculate_qc_metrics(
    adata_filtered,
    qc_vars=['mt'],
    percent_top=None,
    log1p=False,
    inplace=True
)

# 可视化过滤后的分布
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

sc.pl.violin(
    adata_filtered,
    ['n_genes_by_counts'],
    jitter=0.4,
    multi_panel=False,
    ax=axes[0],
    show=False
)
axes[0].set_title('过滤后：检测到的基因数')

sc.pl.violin(
    adata_filtered,
    ['n_counts'],
    jitter=0.4,
    multi_panel=False,
    ax=axes[1],
    show=False
)
axes[1].set_title('过滤后：总 UMI 数')

sc.pl.violin(
    adata_filtered,
    ['pct_counts_mt'],
    jitter=0.4,
    multi_panel=False,
    ax=axes[2],
    show=False
)
axes[2].set_title('过滤后：线粒体基因比例')

plt.tight_layout()
plt.show()

期望结果：

分布更加集中
极端值明显减少
线粒体比例降低

步骤 6：保存结果 [2 分钟]

# 保存过滤后的数据
adata_filtered.write('pbmc_filtered.h5ad')

# 保存 QC 报告（可选）
qc_report = {
    '原始细胞数': adata.n_obs,
    '过滤后细胞数': adata_filtered.n_obs,
    '保留率': f"{adata_filtered.n_obs / adata.n_obs * 100:.1f}%",
    '过滤前基因数': adata.n_vars,
    '过滤后基因数': adata_filtered.n_vars,
    '使用的阈值': {
        'n_genes_MAD': N_MADS,
        'n_counts_MAD': N_MADS,
        'pct_mt_MAD': N_MADS
    }
}

import json
with open('qc_report.json', 'w') as f:
    json.dump(qc_report, f, indent=2)

print("QC 完成！结果已保存。")
print(f"- 过滤后数据：pbmc_filtered.h5ad")
print(f"- QC 报告：qc_report.json")

完整代码清单

# ============ 单细胞 QC 完整流程 ============
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt

# 1. 加载数据
adata = sc.read_h5ad('pbmc_raw.h5ad')

# 2. 识别特殊基因
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['ribo'] = adata.var_names.str.startswith(('RPS', 'RPL'))

# 3. 计算 QC 指标
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt', 'ribo'], inplace=True)

# 4. MAD 过滤
N_MADS = 3
def mad_outliers(adata, metric, n_mads=3, direction='both'):
    data = adata.obs[metric]
    median = np.median(data)
    mad = np.median(np.abs(data - median))
    if direction in ['both', 'lower']:
        lower = median - n_mads * mad
    else:
        lower = -np.inf
    if direction in ['both', 'upper']:
        upper = median + n_mads * mad
    else:
        upper = np.inf
    return (data < lower) | (data > upper)

filter_mask = (
    ~mad_outliers(adata, 'n_genes_by_counts', N_MADS, 'lower') &
    ~mad_outliers(adata, 'n_genes_by_counts', N_MADS, 'upper') &
    ~mad_outliers(adata, 'n_counts', N_MADS, 'lower') &
    ~mad_outliers(adata, 'n_counts', N_MADS, 'upper') &
    ~mad_outliers(adata, 'pct_counts_mt', N_MADS, 'upper')
)

# 5. 应用过滤
adata_filtered = adata[filter_mask].copy()

# 6. 过滤低表达基因
sc.pp.filter_genes(adata_filtered, min_cells=10)

# 7. 保存
adata_filtered.write('pbmc_filtered.h5ad')
print(f"QC 完成！保留 {adata_filtered.n_obs} 个细胞")

常见问题

❌ 问题 1：找不到线粒体基因

症状：adata.var['mt'].sum() 返回 0

解决方案：

# 检查基因名格式
print(adata.var_names[:20])

# 可能的问题：
# 1. 基因名是 Ensembl ID（如 ENSG00000244734）
#    → 需要转换为基因符号
# 2. 基因名全小写（如 mt-co1）
#    → 使用小写前缀
adata.var['mt'] = adata.var_names.str.lower().str.startswith('mt-')

❌ 问题 2：过滤后保留率过低（< 50%）

可能原因：

数据质量差
N_MADS 参数过严
组织特异性（如某些组织本身线粒体含量高）

解决方案：

# 放宽阈值（使用 5 个 MAD）
N_MADS = 5

# 或仅过滤极端异常值（使用 percentile）
lower = np.percentile(adata.obs['n_genes_by_counts'], 1)
upper = np.percentile(adata.obs['n_genes_by_counts'], 99)

下一步

恭喜你完成第一次单细胞 QC！现在你应该能够：

✅ 计算和解读 QC 指标
✅ 使用 MAD 方法进行自适应过滤
✅ 可视化和验证过滤结果

在第三章中，我们将深入学习：

不同组织类型的 QC 策略
高级过滤方法（如双细胞检测）
自动化 QC 工作流

→ 继续阅读：第三章 - 功能详解