质控
scRNA-seq数据中的低质量文库可能有多种来源，例如解离过程中的细胞损伤或文库制备失败（例如无效的逆转录或PCR扩增）。这些通常表现为总计数低，表达的基因很少，线粒体或spike-in比例高的“细胞”。这些低质量的库是有问题的，因为它们可能导致下游分析中的误导性结果：
1.它们形成了自己独特的细胞群，使对结果的解释变得复杂。最明显的原因是细胞损伤后线粒体比例增加或核RNA富集。在最坏的情况下，由不同细胞类型产生的低质量文库可以基于损伤诱导表达谱中的相似性聚集在一起，从而在其他不同亚群之间形成人为的中间状态或轨迹。此外，由于转换后平均值的变化，非常小的库可以形成自己的集群。
2.他们在方差估计或主成分分析过程中扭曲了集群异质性的特征。前几个主要成分将捕获质量差异而不是生物学差异，从而降低降维效果。同样，差异最大的基因将由低质量细胞与高质量细胞之间的差异驱动。最明显的例子：计数非常低的低质量文库，其中归一化放大了那些库中恰好具有非零计数的基因的表观变异。
3.它们包含的基因似乎由于主动缩放以针对小文库进行标准化而被强烈“上调”。这对于污染的以低但恒定水平存在于所有文库中的转录本（例如，来自环境溶液）最成问题。在低质量文库中增加的缩放比例会以较大的标准化表达值将这些转录物的计数转换为少量，从而导致与其他细胞相比明显的上调。这可能会产生误导，因为受影响的基因通常具有生物学敏感性，但实际上在另一个亚群中表达。
为了避免（或至少减轻）这些问题，我们需要在分析开始时删除这些细胞。此步骤通常称为细胞层面的质量控制（QC）。 （这里我们将互换使用“库”和“细胞”，尽管在处理基于液滴的数据时他们的区别将变得很重要。）我们将演示使用A. T. L. Lun等人的小型scRNA-seq数据集。 该版本没有事先进行质量控制，因此我们可以应用。
#--- loading ---#
library(scRNAseq)
sce.416b <- LunSpikeInData(which="416b")
sce.416b$block <- factor(sce.416b$block)
sce.416b
## class: SingleCellExperiment
## dim: 46604 192
## metadata(0):
## assays(1): counts
## rownames(46604): ENSMUSG00000102693 ENSMUSG00000064842 ...
## ENSMUSG00000095742 CBFB-MYH11-mcherry
## rowData names(1): Length
## colnames(192): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
## colData names(9): Source Name cell line ... spike-in addition block
## reducedDimNames(0):
## altExpNames(2): ERCC SIRV
质控指标的选择
我们使用几种常见的QC指标，根据表达谱来鉴定低质量的细胞。这些指标将在下面针对SMART-seq2数据的reads进行质控，但相同的过程适用于由其他技术（如MARS-seq和基于液滴的方法）生成的UMI数据。
1.库大小定义为每个细胞的所有相关feature的总计数之和。在这里，我们将相关feature视为内源基因。具有小文库的细胞质量低下，因为在文库制备过程中某个时刻由于细胞裂解或无效的cDNA捕获和扩增导致RNA丢失。
2.每个细胞中表达feature的数量定义为该细胞计数为非零的内源基因的数量。由于尚未成功捕获多样化的转录物种群，任何具有很少表达基因的细胞都可能质量较差。
3.计算“映射到spike-in转录本的reads数量”相对于“每个细胞所有feature（包括spike-in）reads的总数”的比例。由于应该向每个细胞中添加相同量的spike-in RNA，因此spike-in计数的任何富集都是内源RNA丢失的征兆。因此，高比例表明劣质细胞，可能由于部分细胞裂解或解离过程中的RNA降解而丢失了内源RNA。
- 在没有spike-in转录物的情况下，可以使用定位到线粒体基因组中基因的reads的比例。高比例表明细胞质量较差，可能是由于穿孔细胞的细胞质RNA丢失。理由是，在存在适度损害的情况下，细胞膜上的孔允许单个转录物分子外排，但过小而无法使线粒体逸出，从而导致线粒体转录物相对富集。对于单核RNA-seq实验，高比例也很有用，因为它们可以标记尚未成功剥离细胞质的细胞。
对于每个细胞，我们使用来自scater包的perCellQCMetrics（）函数来计算这些QC指标。sum列包含每个细胞的总数，detected列包含检测到的基因数。subsets_Mito_percent列包含映射到线粒体转录本的reads的百分比。 （出于演示目的，我们展示了两种确定每个转录本的基因组位置的方法。）最后，altexps_ERCC_percent列包含映射到ERCC转录本的reads的百分比。
# Retrieving the mitochondrial transcripts using genomic locations included in
# the row-level annotation for the SingleCellExperiment.
location <- rowRanges(sce.416b)
is.mito <- any(seqnames(location)=="MT")
# ALTERNATIVELY: using resources in AnnotationHub to retrieve chromosomal
# locations given the Ensembl IDs; this should yield the same result.
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
chr.loc <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
keytype="GENEID", column="SEQNAME")
is.mito.alt <- which(chr.loc=="MT")
library(scater)
df <- perCellQCMetrics(sce.416b, subsets=list(Mito=is.mito))
df
## DataFrame with 192 rows and 12 columns
## sum detected subsets_Mito_sum
## <numeric> <numeric> <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 865936 7618 78790
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 1076277 7521 98613
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 1180138 8306 100341
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 1342593 8143 104882
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 1668311 7154 129559
## ... ... ... ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 776622 8174 48126
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 1299950 8956 112225
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 1800696 9530 135693
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 46731 6649 3505
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 1866692 10964 150375
## subsets_Mito_detected
## <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 20
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 20
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 19
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 20
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 22
## ... ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 20
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 25
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 23
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 16
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 29
## subsets_Mito_percent altexps_ERCC_sum
## <numeric> <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 9.09882 65278
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 9.16242 74748
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 8.50248 60878
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 7.81190 60073
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 7.76588 136810
## ... ... ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 6.19684 61575
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 8.63302 94982
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 7.53559 113707
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 7.50037 7580
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 8.05569 48664
## altexps_ERCC_detected
## <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 39
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 40
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 42
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 42
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 44
## ... ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 39
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 41
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 40
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 44
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 39
## altexps_ERCC_percent altexps_SIRV_sum
## <numeric> <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 6.80658 27828
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 6.28030 39173
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 4.78949 30058
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 4.18567 32542
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 7.28887 71850
## ... ... ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 7.17620 19848
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 6.65764 31729
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 5.81467 41116
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 13.48898 1883
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 2.51930 16289
## altexps_SIRV_detected
## <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 7
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 7
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 7
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 7
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 7
## ... ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 7
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 7
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 7
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 7
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 7
## altexps_SIRV_percent total
## <numeric> <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 2.90165 959042
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 3.29130 1190198
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1 2.36477 1271074
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1 2.26741 1435208
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1 3.82798 1876971
## ... ... ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1 2.313165 858045
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1 2.224004 1426661
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1 2.102562 1955519
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 3.350892 56194
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1 0.843271 1931645
或者，有的人可能更喜欢使用addPerCellQC（）函数。 这将计算每个细胞的QC统计信息并将其附加到SingleCellExperiment对象的colData中，从而使我们能够将所有相关信息保留在单个对象中，以供以后处理。
sce.416b <- addPerCellQC(sce.416b, subsets=list(Mito=is.mito))
colnames(colData(sce.416b))
## [1] "Source Name" "cell line"
## [3] "cell type" "single cell well quality"
## [5] "genotype" "phenotype"
## [7] "strain" "spike-in addition"
## [9] "block" "sum"
## [11] "detected" "subsets_Mito_sum"
## [13] "subsets_Mito_detected" "subsets_Mito_percent"
## [15] "altexps_ERCC_sum" "altexps_ERCC_detected"
## [17] "altexps_ERCC_percent" "altexps_SIRV_sum"
## [19] "altexps_SIRV_detected" "altexps_SIRV_percent"
## [21] "total"
此处的关键假设是QC指标与每个细胞的生物学状态无关。 差异（例如，低文库大小，高线粒体比例）被认为是由技术因素而不是生物学过程驱动的，这意味着随后去除的细胞不会在下游分析中歪曲生物学过程。 严重违反此假设可能会导致细胞类型丢失，例如该实验体系本身具有较低的RNA含量或大量的线粒体。 我们可以使用其他诊断工具来检查此类现象（后续高级分析中讲解）。
鉴定低质量细胞
识别低质量单元的最简单方法是在QC指标上应用阈值。 例如，如果库大小小于100,000reads，我们可能会认为这些单元的质量较低； 表达少于5,000个基因； spike-in率超过10％; 或线粒体比例超过10％。
c.lib <- df$sum < 1e5
qc.nexprs <- df$detected < 5e3
qc.spike <- df$altexps_ERCC_percent > 10
qc.mito <- df$subsets_Mito_percent > 10
discard <- qc.lib | qc.nexprs | qc.spike | qc.mito
# Summarize the number of cells removed for each reason.
DataFrame(LibSize=sum(qc.lib), NExprs=sum(qc.nexprs),
SpikeProp=sum(qc.spike), MitoProp=sum(qc.mito), Total=sum(discard))
## DataFrame with 1 row and 5 columns
## LibSize NExprs SpikeProp MitoProp Total
## <integer> <integer> <integer> <integer> <integer>
## 1 3 0 19 14 33
虽然简单，但此策略需要大量经验才能确定每种实验方案和生物系统的合适阈值。 基于reads计数的数据的阈值根本不适用于基于UMI的数据，反之亦然。 线粒体活性或总RNA含量的差异要求分别针对不同的生物系统不断调整线粒体阈值和spike-in阈值。 确实，即使使用相同的方法和系统，由于cDNA捕获效率和每细胞测序深度的差异，适当的阈值也可能因运行而异。
自适应阈值
识别异常值
为了获得合适的阈值，我们假设大多数数据集由高质量细胞组成。然后，我们根据所有细胞中每个指标的中位数的中位数绝对偏差（MAD），确定对于各种QC指标而言异常的细胞。具体来说，如果某个值在“问题”方向上距离中位数多于3个MAD，则认为该值是异常值。这种过滤器将保留遵循正态分布的99％的非异常值。
对于416B数据，我们确定对数转换的库大小比中位数低3个MAD的细胞。对数转换用type =“ lower”时以较小的值提高分辨率。具体而言，它保证阈值不是负值，这对于非负矩阵而言将毫无意义。此外，文库大小的分布表现出较重的右尾并不罕见。对数转换可避免MAD膨胀，从而可能损害左尾的异常检测。 （更一般而言，它证明上述99％的理由是合理的。）
qc.lib2 <- isOutlier(df$sum, log=TRUE, type="lower")
我们对经对数转化的表达基因进行相同的操作
qc.nexprs2 <- isOutlier(df$detected, log=TRUE, type="lower")
isOutlier() 还会在输出向量的属性中返回每个指标的确切过滤阈值。 这些对于检查自动选择的阈值是否合适非常有用。
attr(qc.lib2, "thresholds")
## lower higher
## 434083 Inf
attr(qc.nexprs2, "thresholds")
## lower higher
## 5231 Inf
我们为具有相同功能的基于比例的指标识别异常值。 这些分布经常显示出较重的右尾，但是与前两个指标不同，正是右尾本身包含了假定的低质量细胞。 因此，我们不执行任何变换来缩小尾巴-而是希望将尾巴中的细胞识别为较大的离群值。 （虽然理论上有可能获得高于100％的毫无意义的阈值，但这很少见，因此不会引起实际关注。）
qc.spike2 <- isOutlier(df$altexps_ERCC_percent, type="higher")
attr(qc.spike2, "thresholds")
## lower higher
## -Inf 14.15
qc.mito2 <- isOutlier(df$subsets_Mito_percent, type="higher")
attr(qc.mito2, "thresholds")
## lower higher
## -Inf 11.92
对于这些指标中的任何一个异常的细胞都被认为质量低下并被丢弃。
discard2 <- qc.lib2 | qc.nexprs2 | qc.spike2 | qc.mito2
# Summarize the number of cells removed for each reason.
DataFrame(LibSize=sum(qc.lib2), NExprs=sum(qc.nexprs2),
SpikeProp=sum(qc.spike2), MitoProp=sum(qc.mito2), Total=sum(discard2))
## DataFrame with 1 row and 5 columns
## LibSize NExprs SpikeProp MitoProp Total
## <integer> <integer> <integer> <integer> <integer>
## 1 4 0 1 2 6
或者，可以使用quickPerCellQC（）函数将整个过程一步完成。
reasons <- quickPerCellQC(df,
sub.fields=c("subsets_Mito_percent", "altexps_ERCC_percent"))
colSums(as.matrix(reasons))
## low_lib_size low_n_features high_subsets_Mito_percent
## 4 0 2
## high_altexps_ERCC_percent discard
## 1 6
使用此策略，阈值可适应给定指标的值分布的位置和分布。 这使得QC程序可以适应测序深度，cDNA捕获效率，线粒体含量等方面的变化，而无需任何用户干预或事先经验。 但是，它确实需要一些假设，下面将对其进行详细讨论。
离群值检测的假设
离群检测假定大多数细胞具有可接受的质量。这通常是合理的，并且在某些情况下可以通过肉眼检查细胞是否完整（例如在微孔板上）进行实验支持。如果大多数的细胞质量（不可接受）低，则自适应阈值显然会失败，因为它们无法删除大多数细胞。当然，在旁观者眼中，可接受的与否是随情境改变的——例如，众所周知，神经元很难分解，而我们通常会在某个QC指标下将细胞保留在神经元scRNA-seq数据集的，而在更严格的条件下如胚胎干细胞这个QC指标是不可接受的。
前面提到的另一个假设是，质控指标与每个细胞的生物学状态无关。这在高度异质性细胞群体中最有可能被违反，其中某些细胞类型天然具有更少的总RNA或更多的线粒体。即使没有捕获或测序方面的任何技术问题，也可能将此类细胞视为异常值并删除。通过考虑QC指标中的生物变异性，MAD的使用在某种程度上减轻了该问题。异类种群在高质量细胞之间的指标应具有较高的可变性，从而增加了MAD并减少了错误删除特定细胞类型的机会（以降低去除低质量细胞的能力为代价）。
通常，这些假设是合理的，或者它们的违反对下游结论影响很小。尽管如此，在解释结果时记住它们还是有帮助的。
考虑实验因素
更复杂的研究可能涉及使用不同实验参数（例如测序深度）生成的细胞批次。在这种情况下，自适应策略应分别应用于每个批次。从包含多个批次样品的混合物分布中计算中位数和MAD几乎没有意义。例如，如果一个批次中的测序覆盖率比其他批次低，则它将拖低中位数并使MAD膨胀。这将降低自适应阈值对其他批次的适用性。
如果每个批次都由其自己的SingleCellExperiment表示，则isOutlier（）函数可以直接应用于每个批次，如上所示。但是，如果所有批次中的细胞已合并到单个SingleCellExperiment中，则应该使用batch =参数来确保在每个批次中都识别出异常值。这允许isOutlier（）适应批次之间质量控制指标的系统差异。
我们将再次使用416B数据集进行说明，该数据集包含两个实验因素-原始和癌基因诱导状态。我们将这些因素结合在一起，并通过quickPerCellQC（）在isOutlier（）的batch =参数中使用它。这将导致去除更多的细胞，因为（i）批次之间测序深度的系统差异和（ii）致癌基因诱导表达的基因数量差异不再使MAD膨胀。
batch <- paste0(sce.416b$phenotype, "-", sce.416b$block)
batch.reasons <- quickPerCellQC(df, batch=batch,
sub.fields=c("subsets_Mito_percent", "altexps_ERCC_percent"))
colSums(as.matrix(batch.reasons))
## low_lib_size low_n_features high_subsets_Mito_percent
## 5 4 2
## high_altexps_ERCC_percent discard
## 6 9
也就是说，batch =的使用包含一个更强的假设，即每个批次中的大多数单元都是高质量的。 如果整个批次失败，则异常检测将无法充当该批次的适当QC筛选器。 例如，Grun等人中有两个批次。人类胰腺数据集包含相当一部分推定的受损细胞，其ERCC含量高于其他批次（图1）。 这会使这些批次中的中位数和MAD膨胀，导致无法删除假定的低质量细胞。 在这种情况下，最好从其他批次中计算出一个共享的中位数和MAD，并使用这些估计值来为有问题的批次中的细胞获取适当的过滤阈值，如下所示。
library(scRNAseq)
sce.grun <- GrunPancreasData()
sce.grun <- addPerCellQC(sce.grun)
# First attempt with batch-specific thresholds.
discard.ercc <- isOutlier(sce.grun$altexps_ERCC_percent,
type="higher", batch=sce.grun$donor)
with.blocking <- plotColData(sce.grun, x="donor", y="altexps_ERCC_percent",
colour_by=I(discard.ercc))
# Second attempt, sharing information across batches
# to avoid dramatically different thresholds for unusual batches.
discard.ercc2 <- isOutlier(sce.grun$altexps_ERCC_percent,
type="higher", batch=sce.grun$donor,
subset=sce.grun$donor %in% c("D17", "D2", "D7"))
without.blocking <- plotColData(sce.grun, x="donor", y="altexps_ERCC_percent",
colour_by=I(discard.ercc2))
gridExtra::grid.arrange(with.blocking, without.blocking, ncol=2)
为了确定有问题的批次，一个有用的经验法则是找到质量控制阈值与其他批次的阈值相比异常的批次。 这里的假设是，大多数批次由大多数高质量细胞组成，因此阈值应遵循“典型”批次中的某些单峰分布。 如果我们观察到一个批处理具有极高的阈值，我们可能会怀疑它包含大量会使每个批处理的MAD膨胀的低质量细胞。 我们在下面以Grun等人数据演示此过程。
ercc.thresholds <- attr(discard.ercc, "thresholds")["higher",]
ercc.thresholds
## D10 D17 D2 D3 D7
## 73.611 7.600 6.011 113.106 15.217
names(ercc.thresholds)[isOutlier(ercc.thresholds, type="higher")]
## [1] "D10" "D3"
如果我们不能假设大多数批次都包含大多数高质量细胞，那么所有猜测都将消失； 我们必须恢复选择任意阈值并希望达到最佳效果的方法。
其他方法
另一种策略是基于每个细胞的QC指标来识别高维空间中的异常值。 我们使用来自robustbase的方法基于它们的QC指标来量化每个细胞的“异常”，然后使用isOutlier（）来识别表现出异常高水平异常的低质量细胞。
stats <- cbind(log10(df$sum), log10(df$detected),
df$subsets_Mito_percent, df$altexps_ERCC_percent)
library(robustbase)
outlying <- adjOutlyingness(stats, only.outlyingness = TRUE)
multi.outlier <- isOutlier(outlying, type = "higher")
summary(multi.outlier)
## Mode FALSE TRUE
## logical 179 13
这种方法和相关方法（例如基于PCA的异常检测和支持向量机）可以提供更高的能力，以区分低质量的细胞与高质量的细胞，因为它们可以利用许多QC指标中的模式。 但是，这需要付出一些可解释性的代价，因为删除给定细胞的原因可能并不总是很明显。
为了完整起见，我们注意到也可以从基因表达谱而非质量控制指标中识别异常值。 我们认为这是一种冒险的策略，因为它可以去除稀有细胞群中的高质量细胞。
检查诊断图
检查QC指标的分布是个好习惯（图2），可以发现可能的问题。 在最理想的情况下，我们将看到正态分布可以证明离群值检测中使用的3 MAD阈值是合理的。 其他模型的细胞表明QC指标可能与某种生物学状态相关，有可能导致过滤过程中不同细胞类型的损失。 或细胞亚群的文库制备存在不一致，这在基于plate的实验方案中很常见。 然后，可以快速识别出任何指标的系统差值批次，以进行进一步的故障排除或彻底删除，就像上面的图1一样。
colData(sce.416b) <- cbind(colData(sce.416b), df)
sce.416b$block <- factor(sce.416b$block)
sce.416b$phenotype <- ifelse(grepl("induced", sce.416b$phenotype),
"induced", "wild type")
sce.416b$discard <- reasons$discard
gridExtra::grid.arrange(
plotColData(sce.416b, x="block", y="sum", colour_by="discard",
other_fields="phenotype") + facet_wrap(~phenotype) +
scale_y_log10() + ggtitle("Total count"),
plotColData(sce.416b, x="block", y="detected", colour_by="discard",
other_fields="phenotype") + facet_wrap(~phenotype) +
scale_y_log10() + ggtitle("Detected features"),
plotColData(sce.416b, x="block", y="subsets_Mito_percent",
colour_by="discard", other_fields="phenotype") +
facet_wrap(~phenotype) + ggtitle("Mito percent"),
plotColData(sce.416b, x="block", y="altexps_ERCC_percent",
colour_by="discard", other_fields="phenotype") +
facet_wrap(~phenotype) + ggtitle("ERCC percent"),
ncol=1
)
另一个有用的诊断方法是针对其他一些QC指标绘制线粒体计数比例。 目的是要确认没有总计数和线粒体计数都很高的细胞，以确保我们不会无意中去除代谢活跃的高质量细胞（例如肝细胞）。 我们证明了使用来自涉及老鼠大脑的更大实验的数据； 在这种情况下，我们在图3的右上角没有观察到任何可能与新陈代谢活跃且未受损的细胞相对应的点。
#--- loading ---#
library(scRNAseq)
sce.zeisel <- ZeiselBrainData()
library(scater)
sce.zeisel <- aggregateAcrossFeatures(sce.zeisel,
id=sub("_loc[0-9]+$", "", rownames(sce.zeisel)))
#--- gene-annotation ---#
library(org.Mm.eg.db)
rowData(sce.zeisel)$Ensembl <- mapIds(org.Mm.eg.db,
keys=rownames(sce.zeisel), keytype="SYMBOL", column="ENSEMBL")
sce.zeisel <- addPerCellQC(sce.zeisel,
subsets=list(Mt=rowData(sce.zeisel)$featureType=="mito"))
qc <- quickPerCellQC(colData(sce.zeisel),
sub.fields=c("altexps_ERCC_percent", "subsets_Mt_percent"))
sce.zeisel$discard <- qc$discard
plotColData(sce.zeisel, x="sum", y="subsets_Mt_percent", colour_by="discard")
我们看到所有这些指标相互之间显示出较弱的相关性，大概是细胞损伤的共同潜在作用的体现。 相关性的弱点促使人们使用多种指标来捕获技术质量的不同方面。 当然，不利的一面是，这些指标还可能代表生物学的不同方面，增加了丢弃整个细胞类型的风险。
移除劣质细胞
识别出低质量的细胞后，我们可以选择删除它们或对其进行标记。 删除是最简单的选项，可以通过按列设置SingleCellExperiment来实现。
＃保留我们不想丢弃的列。
filtered <- sce.416b[,!reasons$discard]
质量控制期间最大的实际问题是是否会无意中丢弃整个细胞群。 由于QC指标永远不会完全独立于生物学状态，因此始终存在发生这种情况的风险。 我们可以通过寻找丢弃细胞和保留细胞之间基因表达的系统差异来判断细胞类型是否有丢失。 为了证明这一点，我们计算了416B数据集中废弃池和保留池的平均计数，并计算了池平均值之间的对数倍变化。
# Using the 'discard' vector for demonstration purposes,
# as it has more cells for stable calculation of 'lost'.
lost <- calculateAverage(counts(sce.416b)[,!discard])
kept <- calculateAverage(counts(sce.416b)[,discard])
library(edgeR)
logged <- cpm(cbind(lost, kept), log=TRUE, prior.count=2)
logFC <- logged[,1] - logged[,2]
abundance <- rowMeans(logged)
如果丢弃的池中富集了某种细胞类型，则应观察到相应标记基因的表达增加。 在图4中，丢弃池中没有明显的基因上调系统，这表明QC步骤并未无意间过滤掉416B数据集中的细胞类型。
plot(abundance, logFC, xlab="Average count", ylab="Log-FC (lost/kept)", pch=16)
points(abundance[is.mito], logFC[is.mito], col="dodgerblue", pch=16)
为了进行比较，让我们考虑来自10X Genomics的PBMC数据集的QC步骤。 我们将对库大小应用任意固定的阈值来过滤细胞，而不是使用任何基于异常值的方法。 具体来说，我们会删除所有库大小小于500的库。
#--- loading ---#
library(DropletTestFiles)
raw.path <- getTestFile("tenx-2.1.0-pbmc4k/1.0.0/raw.tar.gz")
out.path <- file.path(tempdir(), "pbmc4k")
untar(raw.path, exdir=out.path)
library(DropletUtils)
fname <- file.path(out.path, "raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
#--- gene-annotation ---#
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)
library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID,
column="SEQNAME", keytype="GENEID")
#--- cell-detection ---#
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]
discard <- colSums(counts(sce.pbmc)) < 500
lost <- calculateAverage(counts(sce.pbmc)[,discard])
kept <- calculateAverage(counts(sce.pbmc)[,!discard])
logged <- edgeR::cpm(cbind(lost, kept), log=TRUE, prior.count=2)
logFC <- logged[,1] - logged[,2]
abundance <- rowMeans(logged)
在图5中，废弃池中存在一个独特的种群，这是一组在丢失中强烈上调的基因。 这包括PF4，PPBP和SDPR，它们（（扰流板警报！）表示存在被alt.discard丢弃的血小板群。
plot(abundance, logFC, xlab="Average count", ylab="Log-FC (lost/kept)", pch=16)
platelet <- c("PF4", "PPBP", "SDPR")
points(abundance[platelet], logFC[platelet], col="orange", pch=16)
如果我们怀疑细胞类型已被我们的质量控制程序错误地丢弃，则最直接的解决方案是放宽质量控制过滤器，以获取与真正的生物学差异相关的指标。例如，可以通过在
isOutlier（）调用中增加nmads =来放宽异常检测。当然，这增加了保留更多低质量细胞的风险。顺便说一句，值得一提的是，细胞的真正技术质量也可能与其类型相关。 （这与细胞类型和QC指标之间的相关性不同，因为后者是我们对质量的不完美代理。）如果某些细胞类型在scRNA-seq协议期间不适合解离或微流体处理，则可能出现这种情况。在这种情况下，如果QC的所有细胞都损坏，则可以“正确”丢弃整个细胞类型。确实，与实验方案中的损失相比，对QC期间细胞类型的计算去除的关注可能很小。
6.6标记劣质细胞
另一种选择是照这样标记低质量的细胞，并将其保留在下游分析中。 此处的目的是允许形成低质量细胞的群，然后在解释结果时识别并忽略这些群。 这种方法避免了丢弃质量控制指标值较差的细胞类型，从而为用户提供了机会来决定此类细胞的群是否代表真正的生物学状态。
marked <- sce.416b
marked$discard <- batch.reasons$discard
缺点是它将QC的负担转移到了细胞群的解释上，这已经是scRNA-seq数据分析的瓶颈。确实，如果我们不相信QC指标，我们将不得不仅基于标记基因来区分真正的细胞类型和低质量的细胞，由于后者倾向于“表达”有趣的基因，因此这并不总是那么容易。保留低质量细胞也会损害方差建模的准确性，例如，需要使用更多的PC来抵消早期PC受低质量细胞和其他细胞之间的差异驱动的事实。
对于常规分析，我们建议默认执行清除操作，以避免低质量细胞引起的并发问题。这样一来，大多数群结构的特征都无需担心或至少不会担心其有效性。完成初始分析后，如果对丢弃的细胞类型有任何疑问，可以在仅标记低质量细胞的情况下进行更彻底的重新分析。这样可以恢复具有低RNA含量，高线粒体比例等的细胞类型，仅在初始分析中“填补空白”的情况下才需要对其进行解释。
