单细胞 RNA-seq 质量控制：功能详解

深入掌握单细胞 QC 的各项技术

本章系统讲解单细胞质控的各个环节，帮助你建立完整的知识体系，应对不同场景的 QC 需求。

3.1 QC 指标详解

3.1.1 测序深度指标

n_counts（总 UMI 数）

定义：每个细胞检测到的唯一分子标识符（UMI）总数

生物学意义：

反映细胞的 RNA 含量
间接反映测序深度
与细胞大小和代谢活性相关

分布特征：

理想分布：
  - 中位数：5,000-20,000（取决于平台和组织）
  - 变异系数（CV）：0.3-0.6
  - 略右偏（存在高转录量细胞）

异常模式：
  - 双峰：可能混入不同细胞类型或技术批次
  - 极左偏：大量低质量细胞
  - 极右偏：存在双细胞或粘连体

过滤策略：

# MAD 方法（推荐）
median = np.median(adata.obs['n_counts'])
mad = np.median(np.abs(adata.obs['n_counts'] - median))
lower = median - 3 * mad
upper = median + 3 * mad

# 百分位数方法（备选）
lower = np.percentile(adata.obs['n_counts'], 2)
upper = np.percentile(adata.obs['n_counts'], 98)

组织差异：

组织类型	典型中位数	过滤阈值建议
PBMC	6,000-10,000	下限：3MAD
肝脏	8,000-15,000	下限：3MAD
肿瘤	5,000-20,000	宽松：5MAD
胚胎干细胞	3,000-7,000	下限：3MAD

n_genes_by_counts（检测到的基因数）

定义：每个细胞中表达量 > 0 的基因数量

生物学意义：

反映细胞的转录复杂性
与细胞分化和状态相关
是细胞质量的核心指标

与 n_counts 的关系：

# 理想情况：强正相关
correlation = np.corrcoef(
    adata.obs['n_counts'],
    adata.obs['n_genes_by_counts']
)[0, 1]

print(f"n_counts vs n_genes 相关性：{correlation:.3f}")
# 期望值 > 0.85

异常模式识别：

import matplotlib.pyplot as plt

plt.scatter(
    adata.obs['n_counts'],
    adata.obs['n_genes_by_counts'],
    alpha=0.3,
    s=10
)

# 标记异常区域
plt.axvline(median_counts - 3*mad_counts, color='red', linestyle='--')
plt.axhline(median_genes - 3*mad_genes, color='red', linestyle='--')

# 异常区域：
# - 左下角：低质量细胞
# - 右上角：双细胞
# - 右下角：可能的技术伪影

3.1.2 细胞活性指标

pct_counts_mt（线粒体基因比例）

生物学背景：

线粒体是细胞的"能量工厂"
线粒体基因的高表达通常反映：
- 细胞凋亡/坏死
- 细胞应激
- 膜破损（细胞裂解）
- 低温休克（某些实验条件）

计算方法：

# Scanpy 自动计算
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)

# 手动计算（验证）
mt_genes = adata.var_names.str.startswith('MT-')
mt_counts = adata[:, mt_genes].X.sum(axis=1).A1
total_counts = adata.X.sum(axis=1).A1
pct_mt_manual = (mt_counts / total_counts) * 100

# 验证
print(f"Scanpy 计算：{adata.obs['pct_counts_mt'].mean():.2f}%")
print(f"手动计算：{pct_mt_manual.mean():.2f}%")

物种差异：

# 人类
adata.var['mt'] = adata.var_names.str.startswith('MT-')

# 小鼠
adata.var['mt'] = adata.var_names.str.startswith('mt-')

# 斑马鱼
adata.var['mt'] = adata.var_names.str.startswith('mt-') | \
                 adata.var_names.str.startswith('MT-')

# 果蝇
adata.var['mt'] = adata.var_names.str.startswith('mt:') | \
                 adata.var_names.isin([
                     'COX1', 'COX2', 'COX3', 'CYTB', 'ATP6', 'ATP8',
                     'ND1', 'ND2', 'ND3', 'ND4', 'ND4L', 'ND5'
                 ])

组织特异性阈值：

组织类型	基线线粒体比例	过滤建议
免疫细胞（PBMC）	3-8%	> 10-12% 过滤
心肌细胞	15-25%	> 30% 过滤
肝细胞	8-15%	> 20% 过滤
肾细胞	5-12%	> 15% 过滤
肿瘤细胞	5-20%	视具体类型
神经元	5-10%	> 12% 过滤

⚠️ 注意：心肌细胞、褐色脂肪细胞等高能量需求组织，线粒体比例天生较高，需要调整阈值。

pct_counts_ribo（核糖体基因比例）

生物学意义：

反映蛋白质合成活性
高值通常表示：
- 细胞增殖（细胞周期活跃）
- 高代谢状态
- 某些分化阶段的特征

核糖体基因识别：

# 哺乳动物
ribosomal_prefixes = ['RPS', 'RPL', 'MRPS', 'MRPL']
adata.var['ribo'] = adata.var_names.str.startswith(
    tuple(ribosomal_prefixes)
)

# 或使用基因列表（更准确）
ribosomal_genes = [
    'RPL1', 'RPL2', 'RPL3', ...  # 核糖体蛋白 L 家族
    'RPS2', 'RPS3', 'RPS4', ...  # 核糖体蛋白 S 家族
]
adata.var['ribo'] = adata.var_names.isin(ribosomal_genes)

应用场景：

# 区分增殖和静息细胞
proliferating = adata[
    adata.obs['pct_counts_ribo'] >
    np.percentile(adata.obs['pct_counts_ribo'], 75)
]

quiescent = adata[
    adata.obs['pct_counts_ribo'] <
    np.percentile(adata.obs['pct_counts_ribo'], 25)
]

print(f"增殖细胞：{len(proliferating)}")
print(f"静息细胞：{len(quiescent)}")

3.1.3 组织特异性指标

血红蛋白基因（血液样本）

# 人类血红蛋白基因
adata.var['hb'] = adata.var_names.str.startswith(('HBA', 'HBB', 'HBD', 'HBE', 'HBG', 'HBM', 'HBQ', 'HBZ'))

# 应用场景
# 1. 红细胞污染检测
# 2. 贫血研究
# 3. 血液发育研究

# 过滤高血红蛋白细胞（如果研究非红细胞类型）
high_hb_mask = adata.obs['pct_counts_hb'] > np.percentile(
    adata.obs['pct_counts_hb'],
    95
)
adata_no_rbc = adata[~high_hb_mask]

其他组织特异性基因

# 肝脏：白蛋白、AFP
adata.var['liver_specific'] = adata.var_names.isin(['ALB', 'AFP'])

# 胰腺：胰岛素、胰高血糖素
adata.var['pancreas_specific'] = adata.var_names.isin(['INS', 'GCG'])

# 神经：MAP2, RBFOX3 (NeuN), GFAP
adata.var['neuronal_markers'] = adata.var_names.isin(['MAP2', 'RBFOX3', 'GFAP'])

3.2 高级过滤策略

3.2.1 双细胞检测

双细胞（Doublets）：两个细胞被封装在同一个液滴中

特征：

异常高的 n_counts 和 n_genes
可能表达两个不同细胞类型的标志基因

检测方法 1：基于统计阈值（简单）

# 使用 MAD 方法
n_genes_median = np.median(adata.obs['n_genes_by_counts'])
n_genes_mad = np.median(np.abs(adata.obs['n_genes_by_counts'] - n_genes_median))

doublet_threshold = n_genes_median + 5 * n_genes_mad  # 使用 5 MAD

doublets = adata.obs['n_genes_by_counts'] > doublet_threshold

print(f"检测到 {doublets.sum()} 个潜在双细胞 ({doublets.sum()/len(doublets)*100:.1f}%)")

检测方法 2：Scrublet（推荐）

# 安装 Scrublet
# pip install scrublet

import scrublet as scr

# 转换为 Scrublet 格式
counts_matrix = adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X

# 创建 Scrublet 对象
scrub = scr.Scrublet(counts_matrix)

# 检测双细胞
doublet_scores, predicted_doublets = scrub.scrub_doublets()

# 添加到 AnnData
adata.obs['doublet_score'] = doublet_scores
adata.obs['predicted_doublet'] = predicted_doublets

# 可视化
scrub.plot_doublet_score_histogram()

# 过滤双细胞
adata_no_doublets = adata[~adata.obs['predicted_doublet']].copy()

3.2.2 空液滴检测（适用于 droplet-based 数据）

问题：液滴可能不含细胞，仅捕获环境 RNA

解决方案 1：EmptyDrops（R）

# 在 R 中运行
library(DropletUtils)

# 读取 10X 数据
counts <- read10xCounts('path/to/10X/output')

# 运行 EmptyDrops
empty <- emptyDrops(counts)

# 过滤
cells_keep <- which(empty$FDR <= 0.01)
filtered_counts <- counts[, cells_keep]

解决方案 2：简单阈值过滤

# 基于总 UMI 数
min_counts = 500  # 根据数据调整
adata = adata[adata.obs['n_counts'] > min_counts].copy()

# 基于基因数
min_genes = 200
adata = adata[adata.obs['n_genes_by_counts'] > min_genes].copy()

3.2.3 细胞周期评分

目的：评估细胞是否受细胞周期效应主导

# 1. 获取细胞周期基因
# Scanpy 提供了小鼠和人类的细胞周期基因列表
import scanpy as sc

# 人类
s_genes = sc.datasets.cell_cycle_genes('human')['s']
g2m_genes = sc.datasets.cell_cycle_genes('human')['g2m']

# 2. 细胞周期评分
sc.tl.score_genes_cell_cycle(
    adata,
    s_genes=s_genes,
    g2m_genes=g2m_genes
)

# 3. 可视化
sc.pl.umap(adata, color=['phase', 'S_score', 'G2M_score'])

# 4. 回归（如果需要）
# sc.pp.regress_out(adata, ['S_score', 'G2M_score'])

3.3 自动化 QC 工作流

3.3.1 完整 QC 管道

def auto_qc(
    adata,
    species='human',
    n_mads=3,
    min_cells=10,
    min_counts=None,
    min_genes=None,
    max_mt_pct=None,
    verbose=True
):
    """
    自动化 QC 流程

    参数：
        adata: AnnData 对象
        species: 'human' 或 'mouse'
        n_mads: MAD 倍数
        min_cells: 过滤基因的最小细胞数
        min_counts: 可选，覆盖 MAD 的下限
        min_genes: 可选，覆盖 MAD 的下限
        max_mt_pct: 可选，覆盖 MAD 的上限
        verbose: 是否打印信息

    返回：
        过滤后的 AnnData 对象
    """
    # 1. 识别特殊基因
    if species == 'human':
        adata.var['mt'] = adata.var_names.str.startswith('MT-')
        ribo_prefixes = ['RPS', 'RPL', 'MRPS', 'MRPL']
    elif species == 'mouse':
        adata.var['mt'] = adata.var_names.str.startswith('mt-')
        ribo_prefixes = ['Rps', 'Rpl', 'Mrps', 'Mrpl']
    else:
        raise ValueError(f"Unknown species: {species}")

    adata.var['ribo'] = adata.var_names.str.startswith(tuple(ribo_prefixes))

    # 2. 计算 QC 指标
    sc.pp.calculate_qc_metrics(
        adata,
        qc_vars=['mt', 'ribo'],
        percent_top=None,
        log1p=False,
        inplace=True
    )

    # 3. 基于 MAD 的过滤
    def mad_filter(data, n_mads, direction='both'):
        median = np.median(data)
        mad = np.median(np.abs(data - median))
        if direction == 'lower':
            threshold = median - n_mads * mad
            return data > threshold
        elif direction == 'upper':
            threshold = median + n_mads * mad
            return data < threshold
        else:
            lower = median - n_mads * mad
            upper = median + n_mads * mad
            return (data > lower) & (data < upper)

    # 应用过滤
    filter_mask = np.ones(adata.n_obs, dtype=bool)

    # n_counts（双向）
    counts_filter = mad_filter(adata.obs['n_counts'], n_mads, 'both')
    filter_mask &= counts_filter

    # n_genes（双向）
    genes_filter = mad_filter(adata.obs['n_genes_by_counts'], n_mads, 'both')
    filter_mask &= genes_filter

    # mt%（仅上限）
    mt_filter = mad_filter(adata.obs['pct_counts_mt'], n_mads, 'upper')
    filter_mask &= mt_filter

    # 可选：覆盖 MAD 阈值
    if min_counts:
        filter_mask &= (adata.obs['n_counts'] >= min_counts)
    if min_genes:
        filter_mask &= (adata.obs['n_genes_by_counts'] >= min_genes)
    if max_mt_pct:
        filter_mask &= (adata.obs['pct_counts_mt'] <= max_mt_pct)

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

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

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

    if verbose:
        print(f"\n基因数：{adata.n_vars} → {adata_filtered.n_vars}")

    return adata_filtered

# 使用
adata_qc = auto_qc(
    adata,
    species='human',
    n_mads=3,
    min_cells=10
)

3.3.2 QC 报告生成

def generate_qc_report(adata, adata_filtered, output_path='qc_report.html'):
    """
    生成 HTML 格式的 QC 报告
    """
    import matplotlib.pyplot as plt
    from matplotlib.backends.backend_pdf import PdfPages

    # 创建 PDF
    with PdfPages(output_path) as pdf:
        # 页面 1：过滤前概览
        fig, axes = plt.subplots(2, 3, figsize=(15, 10))

        # n_genes
        axes[0, 0].violinplot([adata.obs['n_genes_by_counts']])
        axes[0, 0].set_title('n_genes (过滤前)')

        # n_counts
        axes[0, 1].violinplot([adata.obs['n_counts']])
        axes[0, 1].set_title('n_counts (过滤前)')

        # mt%
        axes[0, 2].violinplot([adata.obs['pct_counts_mt']])
        axes[0, 2].set_title('MT% (过滤前)')

        # scatter
        axes[1, 0].scatter(
            adata.obs['n_counts'],
            adata.obs['n_genes_by_counts'],
            c=adata.obs['pct_counts_mt'],
            cmap='viridis',
            alpha=0.5
        )
        axes[1, 0].set_xlabel('n_counts')
        axes[1, 0].set_ylabel('n_genes')

        # histogram
        axes[1, 1].hist(adata.obs['pct_counts_mt'], bins=50)
        axes[1, 1].set_xlabel('MT%')

        # summary stats
        stats_text = f"""
        过滤前统计：
        - 细胞数：{adata.n_obs}
        - 基因数：{adata.n_vars}
        - n_genes 中位数：{adata.obs['n_genes_by_counts'].median():.0f}
        - n_counts 中位数：{adata.obs['n_counts'].median():.0f}
        - MT% 中位数：{adata.obs['pct_counts_mt'].median():.1f}%
        """
        axes[1, 2].text(0.1, 0.5, stats_text, fontsize=12, verticalalignment='center')
        axes[1, 2].axis('off')

        plt.suptitle('QC Report - Before Filtering')
        pdf.savefig(fig)
        plt.close()

        # 页面 2：过滤后概览
        # （类似代码，使用 adata_filtered）

    print(f"QC 报告已保存：{output_path}")

# 使用
generate_qc_report(adata, adata_filtered, 'qc_report.pdf')

3.4 特殊场景的 QC 策略

3.4.1 肿瘤样本

特点：

细胞异质性极高
可能含有循环肿瘤细胞（CTC，低 RNA 含量）
坏死区域产生大量低质量细胞

策略：

# 1. 更宽松的过滤（避免丢失稀有细胞群）
n_mads = 5  # 使用 5 MAD 而非 3

# 2. 不过滤低转录量细胞（可能是静息肿瘤细胞）
# 仅基于 mt% 过滤
filter_mask = (
    (adata.obs['pct_counts_mt'] < np.percentile(adata.obs['pct_counts_mt'], 95))
)

# 3. 后续使用专门工具检测恶性细胞
# 例如 inferCNV、CopyKAT

3.4.2 胚胎/干细胞

特点：

转录活性较低
高度动态的基因表达
稀有亚群频率低

策略：

# 1. 极宽松的过滤
n_mads = 5

# 2. 不设置硬性下限
# 允许低 UMI 细胞存在

# 3. 关注双细胞
# 干细胞容易形成聚集物

3.4.3 空间转录组

特点：

每个 spot 含多个细胞
基因检出数通常高于单细胞
空间位置信息重要

策略：

# 1. 调整阈值（更高）
min_genes = 500  # 单细胞通常是 200
min_counts = 1000

# 2. 不过滤高 UMI（可能是多细胞 spot）
# 仅过滤极低质量

# 3. 保留空间坐标
adata_spatial = adata[filter_mask].copy()

通过本章的学习，你已经掌握了单细胞 QC 的各项技术细节。在第四章中，我们将通过不同物种和组织的真实案例，学习如何灵活应用这些技术。

→ 继续阅读：第四章 - 应用案例