Scanpy 1.10 单细胞质控实战：3指标联合过滤与Scrublet双细胞识别（附Python代码）

1. 单细胞RNA测序质控的核心逻辑

单细胞RNA测序(scRNA-seq)数据的质量控制(QC)是分析流程中至关重要的第一步。与bulk RNA-seq不同，单细胞数据具有更高的技术噪音和生物学异质性，这使得质控环节成为确保下游分析可靠性的关键屏障。

在单细胞数据中，常见的低质量细胞来源包括：

• 破损细胞 ：细胞膜破裂导致RNA大量流失
• 凋亡细胞 ：正在经历程序性死亡的细胞
• 双细胞/多细胞 (doublets/multiplets)：一个液滴中包含多个细胞
• 空液滴 ：未捕获细胞的液滴

这些低质量细胞会引入系统性偏差，影响后续的细胞分群、差异表达分析等结果。因此，建立科学的质控标准对单细胞分析至关重要。

2. 核心质控指标与阈值设定

2.1 三大核心质控指标

在Scanpy中，我们主要依赖三个核心指标进行质控：

指标类型	计算方式	生物学意义	异常表现
每个细胞的基因数	n_genes_by_counts	细胞转录组复杂度	过低：破损细胞；过高：双细胞
每个细胞的UMI总数	total_counts	测序深度	过低：低质量细胞；过高：双细胞
线粒体基因占比	pct_counts_mt	细胞活性	过高：破损或凋亡细胞

线粒体基因的高表达通常意味着细胞应激或凋亡，因为当细胞膜破损时，胞质RNA会流失，而线粒体RNA由于位置固定得以保留。

2.2 基于MAD的自动阈值设定

传统质控通常使用固定阈值，但这种方法无法适应不同数据集的特异性。Scanpy 1.10推荐使用中值绝对偏差(MAD)方法自动计算阈值：

# 计算MAD自动阈值
def mad_threshold(series, n_mads=3):
    median = series.median()
    mad = (series - median).abs().median() * 1.4826
    return median + n_mads * mad

# 应用自动阈值
min_genes = mad_threshold(adata.obs['n_genes_by_counts'], n_mads=-3) 
max_genes = mad_threshold(adata.obs['n_genes_by_counts'], n_mads=3)
min_counts = mad_threshold(adata.obs['total_counts'], n_mads=-3)
max_counts = mad_threshold(adata.obs['total_counts'], n_mads=3)
max_mt = mad_threshold(adata.obs['pct_counts_mt'], n_mads=3)

这种方法能根据数据分布自动调整阈值，特别适合处理不同来源或不同实验条件的数据集。

3. 质控流程实战演示

3.1 数据加载与初步质控

首先加载数据并计算基础指标：

import scanpy as sc
import numpy as np

# 读取数据
adata = sc.read_10x_mtx('filtered_feature_bc_matrix/', var_names='gene_symbols')

# 计算基础QC指标
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)

# 可视化QC指标分布
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True)

3.2 三指标联合过滤

基于计算得到的阈值进行过滤：

# 设置过滤条件
min_genes = 500    # 最少检测到的基因数
max_genes = 6000   # 最多检测到的基因数
max_mt = 15        # 最大线粒体基因占比

# 应用过滤
sc.pp.filter_cells(adata, min_genes=min_genes)
sc.pp.filter_cells(adata, max_genes=max_genes)
adata = adata[adata.obs['pct_counts_mt'] < max_mt, :]

# 记录过滤后数据
print(f"过滤后细胞数: {adata.n_obs}")

3.3 Scrublet双细胞识别

双细胞(doublets)是单细胞测序中的常见问题，Scrublet通过模拟双细胞表达谱来识别真实数据中的双细胞：

# 安装Scrublet
!pip install scrublet

# 运行Scrublet
import scrublet as scr
scrub = scr.Scrublet(adata.X)
doublet_scores, predicted_doublets = scrub.scrub_doublets()

# 将结果添加到adata对象
adata.obs['doublet_score'] = doublet_scores
adata.obs['predicted_doublet'] = predicted_doublets

# 可视化双细胞得分
sc.pl.scatter(adata, x='total_counts', y='doublet_score', color='predicted_doublet')

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

4. 质控后数据标准化与可视化

完成质控后，我们需要对数据进行标准化处理：

# 基本预处理流程
sc.pp.normalize_total(adata, target_sum=1e4)  # CPM标准化
sc.pp.log1p(adata)  # 对数转换
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)  # 筛选高变基因
adata = adata[:, adata.var.highly_variable]  # 保留高变基因

# PCA降维
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver='arpack')

# UMAP可视化
sc.pp.neighbors(adata, n_pcs=30, n_neighbors=20)
sc.tl.umap(adata)
sc.pl.umap(adata, color=['n_genes_by_counts', 'total_counts', 'pct_counts_mt'])

5. 质控结果解读与注意事项

5.1 质控决策要点

1. 阈值设置的灵活性 ：不同组织类型、不同实验protocol可能需要调整阈值
2. 指标间的相关性 ：高UMI数常伴随高基因数，但线粒体占比应独立评估
3. 批次效应影响 ：跨批次数据应分别质控后再整合

5.2 常见问题解决方案

问题现象	可能原因	解决方案
过滤后细胞数过少	阈值设置过严	放宽MAD系数或手动调整阈值
双细胞比例异常高	细胞悬液浓度过高	检查实验记录，考虑重新分析
线粒体基因普遍偏高	样本处理不当	检查样本采集到建库的时间间隔

6. 完整代码实现

以下是整合所有步骤的完整代码：

# 单细胞质控全流程
import scanpy as sc
import numpy as np
import scrublet as scr

# 1. 数据加载
adata = sc.read_10x_mtx('data/', var_names='gene_symbols')

# 2. 计算QC指标
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)

# 3. 自动阈值计算
def mad_threshold(series, n_mads=3):
    median = series.median()
    mad = (series - median).abs().median() * 1.4826
    return median + n_mads * mad

min_genes = mad_threshold(adata.obs['n_genes_by_counts'], n_mads=-3)
max_genes = mad_threshold(adata.obs['n_genes_by_counts'], n_mads=3)
min_counts = mad_threshold(adata.obs['total_counts'], n_mads=-3)
max_counts = mad_threshold(adata.obs['total_counts'], n_mads=3)
max_mt = 15  # 线粒体阈值通常固定

# 4. 细胞过滤
sc.pp.filter_cells(adata, min_genes=min_genes)
sc.pp.filter_cells(adata, max_genes=max_genes)
adata = adata[adata.obs['pct_counts_mt'] < max_mt, :]

# 5. 双细胞检测
scrub = scr.Scrublet(adata.X)
doublet_scores, predicted_doublets = scrub.scrub_doublets()
adata.obs['doublet_score'] = doublet_scores
adata.obs['predicted_doublet'] = predicted_doublets
adata = adata[~adata.obs['predicted_doublet'], :].copy()

# 6. 数据标准化与可视化
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable]

sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver='arpack')
sc.pp.neighbors(adata, n_pcs=30, n_neighbors=20)
sc.tl.umap(adata)

# 7. 结果可视化
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'], jitter=0.4)
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt', color='n_genes_by_counts')
sc.pl.umap(adata, color=['n_genes_by_counts', 'total_counts', 'pct_counts_mt'])

在实际项目中，我们还需要考虑样本特异性、实验批次等因素，可能需要针对不同样本单独进行质控后再整合分析。Scanpy提供的这套质控流程已经过大量实际项目验证，能够为后续分析提供可靠的数据基础。