单细胞RNA-seq生信分析全流程——第四篇：质控

4. 质控Quality Control

4.1 前言

选择适合数据且不会过度校正或消除生物效应的预处理方法至关重要。
由于新的测序技术以及捕获的细胞、测量的基因和识别的细胞群数量的不断增加，用于分析单细胞RNA测序数据的工具集正在快速发展。其中许多工具专用于预处理，旨在解决以下分析步骤：双联体doublet检测、质量控制、归一化、特征feature选择和降维。在本部分中选择的工具，可能严重影响数据的下游分析和解释。例如，如果在质量控制过程中过滤掉过多的细胞，您可能会丢失稀有的细胞亚群，并错过对有趣的细胞生物学的深入了解。然而，如果指定的标准过于宽松，并且没有在预处理步骤中排除质量较差的细胞，则可能会很难注释细胞。因此，选择能够提供最佳并被证明在下游任务方面优于其他工具的方法非常重要。
本教程的起点是单细胞数据，已按照前面所述的部分进行处理。将数据比对以获得分子计数矩阵，即所谓的计数矩阵或读数矩阵。计数矩阵和读数矩阵之间的差异取决于单细胞文库构建方案中是否包含唯一分子标识符(UMI)。读数和计数矩阵的维数为条形码数 x 转录本数。值得注意的是，这里使用术语“条形码”而不是“细胞”，因为条形码可能错误地标记了多个细胞（双联体）或可能没有标记任何细胞（空滴/孔）。我们将在“双联体检测”部分对此进行详细说明。

4.2 环境配置和数据

我们使用为2021年NeurIPS会议上的单细胞数据集成生成10x Multiome数据集[Luecken et al., 2021]。该数据集捕获了在四个不同地点测量的12名健康人类捐赠者的骨髓单核细胞的单细胞多组学数据，以获得嵌套批次效应。在本教程中，我们将使用一批上述数据集（供体8的样本4）来展示scRNA-seq数据预处理的实战流程。
第一步，我们首先使用scanpy加载数据集。

import numpy as np
import scanpy as sc
import seaborn as sns
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

输出结果：

Variable names are not unique. To make them unique, call `.var_names_make_unique`.
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'

读取数据后，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'

数据集的形状为n_obs 16,934 x n_vars 36,601。这转化为barcodes x number of transcripts。我们还检查.var中有关gene_ids (Ensembl Id)、feature_types和基因组的更多信息。
大多数后续分析任务假设数据集中的每个观测值代表来自一个完整单细胞的测量值。在某些情况下，低质量细胞、无细胞RNA或双联体的污染可能会违反这一假设。本教程将指导您如何纠正和消除这种违规行为并获得高质量的数据集。

4.3 过滤低质量的读数reads

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

1.每个条形码的计数数量（计数深度）
2.每个条形码的基因数量
3.每个条形码的线粒体基因计数分数

在细胞QC中，这些协变量通过阈值过滤，因为它们可能对应于死亡细胞。正如所指出的，它们可能反映了膜破裂的细胞，其细胞质mRNA已泄漏，因此只有线粒体中的mRNA仍然存在。这些细胞可能会显示出较低的计数深度、很少检测到的基因以及较高比例的线粒体读数。然而，共同考虑三个QC协变量至关重要，否则可能会导致细胞信号的误解。例如，线粒体计数相对较高的细胞可能参与呼吸过程，不应被过滤掉。计数低或高的细胞可能对应于静止细胞群或尺寸较大的细胞。因此，优先考虑多个协变量再做出阈值决策。一般来说，建议排除较少的细胞，并尽可能避免过滤掉活细胞群或小亚群。
仅对少数或小型数据集进行QC通常以手动方式执行，方法是查看不同QC协变量的分布并识别随后将被过滤的异常值。然而，随着数据集规模的增长，这项任务变得越来越耗时，可能值得考虑通过MAD（median absolute deviations中值绝对偏差）进行自动阈值处理。如果细胞相差5 MADs，我们将其标记为异常值，这是一种相对宽松的过滤策略。我们想强调的是，在细胞注释后重新评估过滤可能是合理的。
在QC中，第一步是计算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)]"))

现在我们可以用scanpy计算对应的QC指标。

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', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts'

正如我们所看到的，该函数向.var和.obs添加了几个附加列。我们想在这里重点介绍其中的一些，有关不同指标的更多信息可以在scanpy文档中找到：

.obs中的n_genes_by_counts是细胞中计数为正的基因数量。
total_counts是细胞的计数总数，这也可能称为库大小。
pct_counts_mt是线粒体细胞总数的比例。
现在，我们绘制每个样本的三个QC协变量n_genes_by_counts、total_counts和pct_counts_mt，以评估各个细胞的捕获情况。

p1 = sns.displot(adata.obs["total_counts"], bins=100, kde=False)
# sc.pl.violin(adata, 'total_counts')
p2 = sc.pl.violin(adata, "pct_counts_mt")
p3 = sc.pl.scatter(adata, "total_counts", "n_genes_by_counts", color="pct_counts_mt")

输出结果：

... storing 'feature_types' as categorical
... storing 'genome' as categorical

这些图表明，一些读数具有相对较高百分比的线粒体计数，这通常与细胞裂解相关。但由于每个细胞的计数数量足够高，并且大多数细胞的线粒体读数百分比低于 20%，我们仍然可以继续下一步分析数据。基于这些图，现在还可以定义用于过滤细胞的手动阈值。在这里，我们将展示基于MAD的自动阈值和过滤的QC。
首先，我们定义一个metric函数，即.obs中的一列和过滤策略中仍然允许的MAD(nmad) 数量。

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_genesQC协变量，每个协变量的阈值为5 MADs。

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()

输出结果：

False    16065
True       869
Name: outlier, dtype: int64

pct_counts_Mt使用3 MADs进行过滤。此外，线粒体计数百分比超过8%的细胞也会被滤除。

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

输出结果：

False    15240
True      1694
Name: mt_outlier, dtype: int64

现在，我们根据这两个附加列来过滤AnnData对象。

print(f"Total number of cells: {adata.n_obs}")
adata = adata[(~adata.obs.outlier) & (~adata.obs.mt_outlier)].copy()

print(f"Number of cells after filtering of low quality cells: {adata.n_obs}")

输出结果：

Total number of cells: 16934
Number of cells after filtering of low quality cells: 14814

p1 = sc.pl.scatter(adata, "total_counts", "n_genes_by_counts", color="pct_counts_mt")

4.4 环境RNA校正

对于基于液滴的单细胞RNA-seq实验，稀释液中存在一定量的背景mRNA，这些mRNA会与细胞一起分布到液滴中并与其一起测序。其最终效果是产生背景污染，该背景污染代表的表达不是来自液滴内包含的细胞，而是来自包含细胞的溶液。
基于液滴的scRNA-seq生成跨多个细胞的基因的唯一分子标识符 (UMI) 计数，旨在识别每个基因和每个细胞的分子数量。它假设每个液滴都含有来自单个细胞的mRNA。双联体、空滴和无细胞RNA可能违反这一假设。无细胞mRNA分子代表稀释液中存在的背景mRNA。这些分子沿着液滴分布，并与它们一起测序。输入溶液中无细胞mRNA的这种污染通常称为“汤”，是由细胞裂解产生的。

无细胞mRNA分子（也称为环境RNA）可能会混淆观察到的计数数量，并可被视为背景污染。纠正基于液滴的scRNA-seq数据集以获得无细胞mRNA非常重要，因为它可能会扭曲我们下游分析中数据的解释。一般来说，每个输入溶液的汤都不同，并且取决于数据集中各个细胞的表达模式。去除环境mRNA的方法（例如SoupX和 DecontX）旨在估计汤的成分，并根据汤的表达校正计数表。
第一步，SoupX计算汤的profile。它根据未过滤的Cellranger矩阵给出的空液滴估计环境mRNA表达谱。接下来，SoupX估计细胞特定的污染分数。最后，它根据环境mRNA表达谱和估计的污染来校正表达矩阵。
SoupX的输出是修改后的计数矩阵，可用于任何下游分析工具。
我们现在加载运行SoupX所需的python和R包。

import anndata2ri
import logging

import rpy2.rinterface_lib.callbacks as rcb
import rpy2.robjects as ro

rcb.logger.setLevel(logging.ERROR)
ro.pandas2ri.activate()
anndata2ri.activate()

%load_ext rpy2.ipython

%%R
library(SoupX)

SoupX可以在没有聚类信息的情况下运行，但是如果提供基本聚类，结果会更好。SoupX可以与cellranger生成的默认cluster一起使用，也可以通过手动定义cluster来使用。我们将在本篇中展示后者，因为SoupX的结果对所使用的聚类并不强烈敏感。
现在，我们创建AnnData对象的副本，对其进行标准化，降低其维度，并在处理后的副本上计算默认leiden簇。后续章节将更详细地介绍聚类。现在我们只需要知道leiden聚类为我们提供了数据集中细胞的分区（社区）。我们将获得的簇保存为soupx_groups并删除AnnData对象的副本以节省内存。
首先，我们生成AnnData对象的副本，对其进行标准化和log1p转换。此时我们使用简单的移位对数归一化。

adata_pp = adata.copy()
sc.pp.normalize_per_cell(adata_pp)
sc.pp.log1p(adata_pp)

接下来，我们计算数据的主成分以获得较低维的形式。然后使用该形式来生成数据的邻域图并在KNN 图上运行leiden聚类。我们将簇作为soupx_groups添加到.obs并将它们保存为向量。

sc.pp.pca(adata_pp)
sc.pp.neighbors(adata_pp)
sc.tl.leiden(adata_pp, key_added="soupx_groups")

# Preprocess variables for SoupX
soupx_groups = adata_pp.obs["soupx_groups"]

现在，我们可以删除AnnData对象的副本，因为我们生成了可在soupX中使用的簇向量。

del adata_pp

接下来，我们保存细胞名称、基因名称和过滤后的cellranger输出的数据矩阵。SoupX需要features x barcodes的矩阵，因此我们必须转置.X。

cells = adata.obs_names
genes = adata.var_names
data = adata.X.T

SoupX还需要细胞矩阵的原始基因，该矩阵在cellranger输出中通常称为raw_feature_bc_matrix.h5。我们像之前一样使用scanpy加载filtered_feature_bc_matrix.h5并在对象上运行.var_names_make_unique()并转置相应的.X。

adata_raw = sc.read_10x_h5(
    filename="raw_feature_bc_matrix.h5",
    backup_url="https://figshare.com/ndownloader/files/39546217",
)
adata_raw.var_names_make_unique()
data_tod = adata_raw.X.T

输出结果：

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

del adata_raw

现在，我们已准备好运行SoupX的一切。输入是经过过滤的条形码 x 细胞的cellranger矩阵、来自cellranger的barcodes x droplets液滴原始表、基因和细胞名称以及通过简单的leiden聚类获得的簇。输出将是校正后的计数矩阵。
我们首先从液滴表和细胞表构建一个SoupChannel。接下来，我们将元数据添加到SoupChannel对象中，该对象可以是data.frame形式的任何元数据。我们在这里添加:

%%R -i data -i data_tod -i genes -i cells -i soupx_groups -o out 

# specify row and column names of data
rownames(data) = genes
colnames(data) = cells
# ensure correct sparse format for table of counts and table of droplets
data <- as(data, "sparseMatrix")
data_tod <- as(data_tod, "sparseMatrix")

# Generate SoupChannel Object for SoupX 
sc = SoupChannel(data_tod, data, calcSoupProfile = FALSE)

# Add extra meta data to the SoupChannel object
soupProf = data.frame(row.names = rownames(data), est = rowSums(data)/sum(data), counts = rowSums(data))
sc = setSoupProfile(sc, soupProf)
# Set cluster information in SoupChannel
sc = setClusters(sc, soupx_groups)

# Estimate contamination fraction
sc  = autoEstCont(sc, doPlot=FALSE)
# Infer corrected table of counts and rount to integer
out = adjustCounts(sc, roundToInt = TRUE)

SoupX成功推断出校正后的计数，我们现在可以将其存储为附加层layer。在以下所有分析步骤中，我们希望使用SoupX校正计数矩阵，因此我们用soupX校正矩阵覆盖.X。

adata.layers["counts"] = adata.X
adata.layers["soupX_counts"] = out.T
adata.X = adata.layers["soupX_counts"]

接下来，我们另外过滤掉少于20个细胞中未检测到的基因，因为这些基因没有提供信息。

print(f"Total number of genes: {adata.n_vars}")

# Min 20 cells - filters out 0 count genes
sc.pp.filter_genes(adata, min_cells=20)
print(f"Number of genes after cell filter: {adata.n_vars}")

输出结果：

Total number of genes: 36601
Number of genes after cell filter: 20171

4.5 双联体检测

双联体被定义为在相同细胞条形码下测序的两个细胞，它们是在同一液滴中捕获的。这就是为什么我们到目前为止使用术语“条形码”而不是“细胞”。如果双联体由相同细胞类型（但来自不同个体）形成，则称为同型，否则称为异型。同型双联体不一定可以从计数矩阵中识别出来，并且通常被认为是无害的，因为它们可以通过细胞散列或SNP来识别。因此，它们的识别不是双联体检测方法的主要目标。
由不同细胞类型或状态形成的双联体称为异型。它们的识别至关重要，因为它们很可能被错误分类，并可能导致下游分析步骤的错误。因此，双联体检测和去除通常是初始预处理步骤。双联体可以通过其大量读取和检测到的特征来识别，也可以通过创建人工双联体并将其与数据集中存在的细胞进行比较的方法来识别。双联体检测方法计算效率高，并且存在多个用于此任务的软件包。
在本教程中，我们将展示scDblFinder R包。scDblFinder随机选择两个液滴，并通过平均它们的基因表达谱来创建人工双联体。然后，双联体得分定义为主成分空间中每个液滴的k最近邻图中人工双联体的比例。

我们首先载入一些python和R包。

%%R
library(Seurat)
library(scater)
library(scDblFinder)
library(BiocParallel)

data_mat = adata.X.T

现在，我们可以使用data_mat作为SingleCellExperiment中scDblFinder的输入来启动双联体检测。scBblFinder向sce的colData添加几列。其中三个主要的组成：

sce$scDblFinder.score：最终的双联体得分（越高，细胞越有可能是双联体）
sce$scDblFinder.ratio：细胞邻域中人工双联体的比例
sce$scDblFinder.class：分类（双联体或单联体）

我们只会输出参数并将其存储在.obs的AnnData对象中。其他参数可以类似地添加到AnnData对象。

%%R -i data_mat -o doublet_score -o doublet_class

set.seed(123)
sce = scDblFinder(
    SingleCellExperiment(
        list(counts=data_mat),
    ) 
)
doublet_score = sce$scDblFinder.score
doublet_class = sce$scDblFinder.class

scDblFinder输出具有分类Singlet(1) 和Doublet (2) 的类。我们将其添加到.obs中的AnnData对象中。

adata.obs["scDblFinder_score"] = doublet_score
adata.obs["scDblFinder_class"] = doublet_class
adata.obs.scDblFinder_class.value_counts()

输出结果：

singlet    11956
doublet     2858
Name: scDblFinder_class, dtype: int64

我们建议暂时将已识别的双联体留在数据集中，并在可视化过程中检查双联体。
在下游聚类过程中，重新评估质量控制和所选参数可能会很有用，以过滤掉更多或更少的细胞。我们现在可以保存数据集并继续之后的标准化过程。

adata.write("s4d8_quality_control.h5ad")