《单细胞最佳实践》—重新回顾质控部分的新理解

生信菜鸟团

发布于 2025-08-04 09:57:06

24600

代码可运行

运行总次数：0

代码可运行

在我接触单细胞数据分析的时候，早就已经有一套非常完整的 seurat 的流程了，直接照着官网给的代码示例分析就好了，为什么要这么分析呢我也不清楚，只知道照着做，而且目前单细胞数据有那么多，每篇文章里各种参数都有一些细节上的差别，很多参数都是根据经验来确定的，没有一个明确的标准。为了丰富这方面的知识以及加深理解，《Single Cell Best Practices》这本在线书籍可能有一定的帮助，链接是 https://www.sc-best-practices.org/preamble.html

这本书全面概述了从预处理到可视化和统计评估等基本分析的步骤，而且还在更新中，示例代码以 python 为主（如果想要以 R 为主的可以看 Orchestrating Single-Cell Analysis with Bioconductor 这本书），这本书前面还介绍了一下 scRNA-seq 的原理和原始数据的处理，我就直接跳到质量控制部分。

1.1 Motivation（为什么要进行质控）

scRNA-seq 数据存在 drop-out 现象，这意味着数据是含有过多 0 的稀疏矩阵。如果在质量控制期间过滤掉过多细胞，可能会丢失稀有的细胞亚群，消除掉一些生物学效应；而如果筛选过于宽松，在预处理流程中没有排除低质量细胞，那么在后续注释细胞时可能会遇到困难。因此，选择合适的预处理方法至关重要。

1.2 环境设置与数据

首先是加载数据集，这个数据是托管在 figshare 上的一个健康样本的人类骨髓单核细胞数据。

import numpy as np
import scanpy as sc
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import median_abs_deviation

sc.settings.verbosity = 0
sc.settings.set_figure_params(
    dpi=80,
    facecolor="white",
    frameon=False,
)

adata = sc.read_10x_h5(
    filename="filtered_feature_bc_matrix.h5",
    backup_url="https://figshare.com/ndownloader/files/39546196",
)
adata

UserWarning: Variable names are not unique. To make them unique, call `.var_names_make_unique`

AnnData object with n_obs × n_vars = 16934 × 36601
    var: 'gene_ids', 'feature_types', 'genome', 'interval'

读取数据后，scanpy 显示警告，并非所有变量名都是唯一的。这表明一些变量（基因）出现了不止一次，这可能导致下游分析任务出现错误或意外行为。执行 var_names_make_unique()，通过将数字字符串附加到每个重复的索引元素（“1”、“2” 等）上来使变量名唯一。

adata.var_names_make_unique()
adata

AnnData object with n_obs × n_vars = 16934 × 36601
    var: 'gene_ids', 'feature_types', 'genome', 'interval'

这个数据集是 n_obs × n_vars = 16934 × 36601 ，其中 n_obs 表示 barcodes，n_vars 为 number of transcripts 转录本的数目。这里的 barcodes/observation 在大多数的下游分析中假设数据集中的每个观测值代表来自一个完整单细胞的测量结果，在某些情况下这种假设可能被低质量细胞、细胞外游离RNA或者双细胞所破坏。另外还有 .var 中包括 Ensemble ID、特征类型、基因组和所在位置在内的更多信息，可以简单看一下：

adata.var.head()

1.3 过滤低质量细胞

质量控制的第一步是从数据集中移除低质量细胞。当一个细胞检测到的基因数量少、计数深度低且线粒体计数比例高时，它可能具有破裂的细胞膜，这可能表明细胞正在死亡。由于这些细胞通常不是分析的主要目标，并且可能扭曲下游分析，我们会在质量控制期间将其移除。为了识别它们，我们定义细胞质量控制（QC）阈值。细胞QC通常基于以下三个 QC 协变量进行：

每个条形码的计数数量（计数深度）
每个条形码检测到的基因数量
每个条形码中线粒体基因计数的比例

在细胞 QC 中，这些协变量通过阈值设定进行过滤。低计数深度、检测到的基因少以及高比例的线粒体 reads 可能反映了细胞膜破裂的细胞，其胞质mRNA已泄漏出去，因此只有线粒体中的mRNA仍然存在。然而，必须联合考虑这三个 QC 协变量，否则可能导致对细胞信号的错误解释。例如，具有相对较高线粒体计数比例的细胞可能参与呼吸过程，不应被过滤掉。而具有低或高计数的细胞可能对应于静止细胞群或体积较大的细胞。因此，设定阈值时最好考虑多个协变量。一般来说，建议排除更少的细胞，并尽可能宽松，以避免过滤掉有活力的细胞群或小的亚群。

对于少量或小型数据集的 QC 通常通过手动方式进行，即查看不同 QC 协变量的分布并识别离群值，随后进行过滤。然而，随着数据集规模增大，这项任务变得越来越耗时，可能值得考虑使用MAD（中位数绝对偏差）进行自动阈值设定。

其中 Xi 是观测值 i 对应的 QC 指标。如果细胞偏离中位数超过 5 个 MAD ，我们将其标记为离群值，这是一种相对宽松的过滤策略。需要强调的是，在细胞注释后重新评估过滤策略可能是合理的。

计算 QC 协变量使用的是 scanpy 函数 sc.pp.calculate_qc_metrics，该函数还可以计算特定基因群计数的比例。因此，我们定义线粒体基因、核糖体基因和血红蛋白基因。需要注意的是，根据数据集中考虑的物种，线粒体计数要么以 "mt-" 为前缀，要么以 "MT-" 为前缀进行注释。这里使用的数据集是人类骨髓，因此以 "MT-" 为前缀，小鼠数据集前缀通常是小写 "mt-"。

# mitochondrial genes 线粒体基因
adata.var["mt"] = adata.var_names.str.startswith("MT-")
# ribosomal genes 核糖体基因
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
# hemoglobin genes 血红蛋白基因
adata.var["hb"] = adata.var_names.str.contains("^HB[^(P)]")

sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt", "ribo", "hb"], inplace=True, percent_top=[20], log1p=True
)
adata

AnnData object with n_obs × n_vars = 16934 × 36601
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb'
    var: 'gene_ids', 'feature_types', 'genome', 'interval', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'

可以看到这个函数向 .obs 和 .var 添加了一些额外的列。重点介绍其中的几个（更多关于不同指标的信息可以在 scanpy 文档中找到）：

n_genes_by_counts 是一个细胞中有表达的基因数量。
total_counts 是一个细胞的总计数数量，这也可能被称为文库大小（library size）。
pct_counts_mt 是一个细胞中线粒体计数占总计数的比例。

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

sns.histplot(adata.obs["total_counts"], bins=100, ax=axes[0])
sc.pl.violin(adata, 'pct_counts_mt', ax=axes[1], show=False)
sc.pl.scatter(adata, 'total_counts', 'n_genes_by_counts', ax=axes[2], show=False)

plt.tight_layout()
plt.show()

这些图表明一些细胞具有相对较高的线粒体计数比例，这通常与细胞退化有关。但由于每个细胞的计数数量足够高，且大多数细胞的线粒体 reads 比例低于 20%，我们仍然可以处理这些数据。基于这些图，可以定义用于过滤细胞的手动阈值。接下来，我们将展示基于 MAD 的自动阈值设定和过滤。

首先，我们定义一个函数 is_outlier() ，它有两个参数 metric（.obs 中的一列）和 nmads（过滤允许的 MAD 数量）：

def is_outlier(adata, metric: str, nmads: int):
    M = adata.obs[metric]
    outlier = (M < np.median(M) - nmads * median_abs_deviation(M)) | (
        np.median(M) + nmads * median_abs_deviation(M) < M
    )
    return outlier

我们应用于 log1p_total_counts、log1p_n_genes_by_counts 和 pct_counts_in_top_20_genes 这三个 QC 协变量，每个都使用 5 个 MAD 的阈值。

adata.obs["outlier"] = (
    is_outlier(adata, "log1p_total_counts", 5)
    | is_outlier(adata, "log1p_n_genes_by_counts", 5)
    | is_outlier(adata, "pct_counts_in_top_20_genes", 5)
)
adata.obs.outlier.value_counts()

outlier
False    16065
True       869
Name: count, dtype: int64

pct_counts_mt 使用 3 个 MAD 进行过滤。此外，过滤掉线粒体计数比例超过 8% 的细胞。

adata.obs["mt_outlier"] = is_outlier(adata, "pct_counts_mt", 3) | (
    adata.obs["pct_counts_mt"] > 8
)
adata.obs.mt_outlier.value_counts()

mt_outlier
False    15240
True      1694
Name: count, dtype: int64

adata = adata[(~adata.obs.outlier) & (~adata.obs.mt_outlier)].copy()
sc.pl.scatter(adata, "total_counts", "n_genes_by_counts", color="pct_counts_mt")

过滤掉低质量细胞之后还剩 14814 个细胞。

1.4 环境RNA（cell-free RNA）校正

对于基于液滴的单细胞 RNA 测序实验，稀释液中存在一定量的背景 mRNA，这些 mRNA 与细胞一起被分配到液滴中并进行测序。这样做的最终效果是产生背景污染，该污染代表的不是来自液滴内细胞的表达，而是来自包含细胞的溶液。游离 mRNA 分子代表稀释液中存在的背景 mRNA。输入溶液中游离 mRNA 的这种污染通常称为 “soup”，它是由细胞裂解产生的。

什么时候需要去除环境 RNA 呢？各种实验因素导致样本中死细胞比例过高，观察到明显的“汤效应”（Soup-like Profile），例如在初步分析中，发现许多细胞都表达一些通常只在特定稀有细胞类型中高表达或普遍低表达的基因，或者 QC 指标中每个细胞的线粒体基因比例异常低（因为环境RNA主要来自胞浆 mRNA），这强烈提示环境RNA污染严重。

去除环境 mRNA 的方法有 SoupX 和 DecontX，都是基于 R 的，详细内容可以参考 https://www.sc-best-practices.org/preprocessing_visualization/quality_control.html。

1.5 Doublets 检测

Doublets 被定义为在同一个细胞条形码下测序的两个细胞，它们被捕获在同一个液滴中。如果 doublets 由相同细胞类型（但来自不同个体）形成，就是同型（homotypic）的，否则就是异型（heterotypic）的。同型双细胞不一定可以从计数矩阵中识别，并且通常被认为是无害的，因为它们可以通过 cell hashing 或 SNP 来识别。异性双细胞的识别至关重要，因为它们很可能被错误分类，并可能导致下游分析步骤失真。

《Benchmarking Computational Doublet-Detection Methods for Single-Cell RNA Sequencing Data》这篇文章对九种不同的双细胞检测方法进行了基准测试，并评估了它们在计算效率和双峰检测精度方面的性能，链接：https://www.sciencedirect.com/science/article/pii/S2405471220304592

这本书里展示了用 scDblFinder R 包实现双细胞检测的过程，详情可看 https://www.sc-best-practices.org/preprocessing_visualization/quality_control.html

也可以直接使用 scanpy 包中的 scrublet 方法：

adata = sc.pp.scrublet(adata, random_state=123, expected_doublet_rate=0.05)
sc.external.pl.scrublet_score_distribution(adata)

然后直接过滤掉就可以了

adata_clean = adata[adata.obs["predicted_doublet"] == False, :].copy()

1.6 关键

对低质量细胞的过滤应该基于绝对中位偏差和宽松的阈值，以避免对较小亚群产生偏差。
可以用 scDblFinder 等方法有效识别双细胞。
双细胞检测方法不应该在多个批次的整合 scRNA-seq 数据上运行。
Feature-based filtering does not show benefits for downstream tasks. 基于特征的过滤对下游任务没有显示出益处。“Feature-based filtering” 应该指的是 sc.pp.filter_genes(adata, min_cells=3) 这种根据特定标准移除表达量过低或在太少细胞中表达的基因，这里的核心意思指的应该是，这种基于基因的过滤步骤未能被证明能显著改善后续分析任务（如下游任务）的质量或结果。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2025-08-01，如有侵权请联系 cloudcommunity@tencent.com 删除

var