对应原版教程第6章 http://bioconductor.org/books/release/OSCA/overview.html 在单细胞数据分析中的第一步质控往往是剔除不合格的细胞。本小节主要介绍了如何鉴定出低质量细胞的相关知识点。
笔记要点
- 1、为什么要质控
- 1.1 低质量细胞文库的特点
- 1.2 低质量细胞文库的影响
- 2、质控指标的选择与计算
- 2.1 质控指标选择的假设
- 2.2 常见的指标
- 2.3 计算方法
- 3 确定低质量细胞的指标阈值
- 3.1 固定阈值
- 3.2 outliner确定异常细胞
- 3.3 遇到batch批次效应的解决方法
- 4 根据阈值筛选细胞
- 4.1 可视化阈值指标
- 4.2 核实剔除细胞是否包含一种特定的细胞类型
- 4.3 剔除低质量细胞
1、为什么要质控
1.1 低质量细胞文库的特点
- scRNA-seq的结果是col为细胞标签，row为基因名的表达矩阵。一列即代表一个细胞的文库测序情况。
- 低质量的细胞文库一般有3个特征： （1）总counts数很少，即列总和较低； （2）非零基因很少，即相当部分基因的count为0； （3）spike-in或者线粒体基因表达相对较高。
- 低质量细胞的测序结果主要是实验过程中细胞发生裂解，低效的逆转录、PCR扩增等因素造成的。
1.2 低质量细胞文库的影响
（1）无意义的cluster
这些低质量细胞因为“低质量的相似性”在聚类时聚成一个单独的cluster，但这个cluster本身是没有任何生物意义的；反而会干扰后续的差异分析、细胞注释、拟时序分析等步骤
（2）异常的异质性
在后续的挑选高变基因、PCA主成分分析等。这些细胞的基因表达会低于其它细胞的平均水平，“贡献”出非生物学意义的异质性(高质量细胞与低质量细胞间的异质性)，从而挑出来的高变基因、计算的Top前几的主成分捕捉的生物学水平的异质性占比下降。
（3）冒牌的高表达基因
细胞或者说一群细胞的高表达基因是鉴定细胞类型的marker gene。但如果存在一些基因自身表达在所有细胞中均维持在较低水平，但受文库大小的标准化处理后，处于低质量细胞的这些基因表达相对其它正常细胞就会被无端放大。
2、质控指标的选择与计算
2.1 质控指标选择的假设
- 对于所选择用于质控的若干指标，我们其实就假设了细胞的这些指标的波动variance主要由于实验过程造成的，而非生物学水平的因素造成的。
- 基于这个假设，我们可以认为根据这些指标删除的是低质量细胞。如果违反了这一假设(例如某一类细胞的线粒体表达就是相对较高)，就可能误删了部分合格的细胞。
2.2 常见的指标
即根据上面1.1低质量细胞文库的三个特点对应的3/4个指标
- （1）library size 细胞文库大小，即对于表达矩阵的每一列的求和。该值越小，越有可能是低质量细胞。原因如1.1点的第一小点
- （2）the number of expressed gene 细胞表达基因数目，即对于表达矩阵的每一列的非0值的基因数目。该值越小，越有可能是低质量细胞。因为未能捕获到细胞的转录水平的多样性。
- （3）spike-in/线粒体表达相对比例
- 比例越高，越有可能是低质量细胞。
- 对于spike-in外参转录本，由于初始每个细胞加的量都是一样的。如果某些细胞的spike-in表达相比其它细胞异常增大，可能就意味着细胞的内源基因表达相对偏低。
- 对于线粒体基因，在穿孔裂解的细胞里，细胞质RNA较容易游离到胞外，而由于线粒体体积较大，不易“逃出”胞体，从而线粒体基因含量相对稳定。
2.3 计算方法
使用scater包的addPerCellQC函数，可自动计算上述所有指标
library(scRNAseq)
sce.416b <- LunSpikeInData(which="416b")
sce.416b$block <- factor(sce.416b$block)
sce.416b
由于需要计算线粒体基因表达比例，并且基因ID为ensembl格式（小鼠），所以需要注释下基因所处的染色体序列信息。
代码语言：javascript复制# 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)
sce.416b <- addPerCellQC(sce.416b, subsets=list(Mito=is.mito))
colData(sce.416b)
colnames(colData(sce.416b))
# [1] "Source Name" "cell line" "cell type"
# [4] "single cell well quality" "genotype" "phenotype"
# [7] "strain" "spike-in addition" "block"
# [10] "sum" "detected" "percent_top_50"
# [13] "percent_top_100" "percent_top_200" "percent_top_500"
# [16] "subsets_Mito_sum" "subsets_Mito_detected" "subsets_Mito_percent"
# [19] "altexps_ERCC_sum" "altexps_ERCC_detected" "altexps_ERCC_percent"
# [22] "altexps_SIRV_sum" "altexps_SIRV_detected" "altexps_SIRV_percent"
# [25] "total"
如上，最重要的4列分别是sum列即为library size；detected即为表达基因数；subsets_Mito_percent即为线粒体基因表达比例；altexps_ERCC_percent即为spike-in转录本表达比例。
percent_top_x的若干列表示按该细胞的基因表达从高到低，前top x的基因表达总counts占细胞文库的比例大小。
3 确定低质量细胞的指标阈值
3.1 固定阈值
手动设定低质量细胞的指标阈值虽然是最简单、最直接的方法，但需要对所做的细胞类型、状态、实验参数等都需要有准确的把控。
例如针对上述数据，把文库大小阈值设为100000(这是smart-seq，与10X不可比)；表达基因数设为5000,；spike-in以及线粒体基因比例均设为10%（仅作为举例，不可作为标准参考~）。如下结果，会剔除33个cell
代码语言：javascript复制qc.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
3.2 outliner确定异常细胞
(1) outliner的假设
- 首先还是这些使用的QC指标仅表示这些细胞的实验测序水平的误差。
- 其次采用outliner方法，默认大部分细胞是合格的，只有少部分是低质量细胞，即寻找这些少数的离群点。
(2) 计算方法
- 单独计算
qc.lib2 <- isOutlier(df$sum, log=TRUE, type="lower")
attr(qc.lib2, "thresholds")
qc.nexprs2 <- isOutlier(df$detected, log=TRUE, type="lower")
attr(qc.nexprs2, "thresholds")
qc.spike2 <- isOutlier(df$altexps_ERCC_percent, type="higher")
attr(qc.spike2, "thresholds")
qc.mito2 <- isOutlier(df$subsets_Mito_percent, type="higher")
attr(qc.mito2, "thresholds")
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
如上代码分别计算了每个指标里的离群点(低质量细胞)。其中的type参数表示离群点方向，对于前2者，取低；对于后两者取高。至于前两个指标为何取log，如下图所示数据点log前后的变化，在log变换后可发现左边出现较明显的尾巴，从而符合之前大部分细胞正常（同一分布趋势），小部分细胞异常的假设。对于后两者，因为作图发现右边存在尾巴，符合筛选预期，故不做变动
- scater包的quickPerCellQC函数可自动完成上述的计算
reasons <- quickPerCellQC(colData(sce.416b),
percent_subsets=c("subsets_Mito_percent", "altexps_ERCC_percent"))
colSums(as.matrix(reasons))
3.3 遇到batch批次效应的解决方法
- 如果每个batch分别在是一个sce对象，那就分别计算离群点即可；
- 如果一个sce对象里包含多个批次，那么可以拆成多个同批次的sce，也可以使用
scater包的quickPerCellQC函数时设置batch参数 - 例如sce.416b的数据集，有两个批次的细胞，可如下处理
table(sce.416b$block)
batch <- paste0(sce.416b$phenotype, "-", sce.416b$block)
batch.reasons <- quickPerCellQC(df, batch=batch,
percent_subsets=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
如上方法与结果有两个注意点
- （1）首先从筛选结果来看，共剔除9个细胞，高于上面不设batch的6个细胞的结果。原因是考虑了batch因素的影响，降低了相对于原来MAD的各个batch的MAD值，从而在每个batch中会出现更多偏离3倍MAD的细胞指标，即剔除更多细胞。
- （2）使用batch参数的前提假设是每个batch的大部分细胞都是高质量细胞。如果某一个batch出现问题，按照此假设也只能检测出少量的outliner。如何确定存在有问题的批次，作者建议，在存在多个批次的情况下，并且假设大部分batch都是合格的，那么如果少数几个batch的阈值明显偏离其它大部分batch的阈值，那么很有可能就是bad batch。对于存在问题的batch的解决方法，可以通过参考其它合格batch的阈值进行鉴别outliner。相关代码如下
library(scRNAseq)
sce.grun <- GrunPancreasData()
#这个sce里有5个batch，其中有两个是有问题的(ERCC占比过高)
sce.grun <- addPerCellQC(sce.grun)
discard.ercc <- isOutlier(sce.grun$altexps_ERCC_percent,
type="higher", batch=sce.grun$donor)
ercc.thresholds <- attr(discard.ercc, "thresholds")["higher",]
ercc.thresholds #如下结果 D10、D3的阈值存在明显异常
# D10 D17 D2 D3 D7
# 73.610696 7.599947 6.010975 113.105828 15.216956
# 解决办法：综合参考其它正常batch的阈值
discard.ercc2 <- isOutlier(sce.grun$altexps_ERCC_percent,
type="higher", batch=sce.grun$donor,
subset=sce.grun$donor %in% c("D17", "D2", "D7"))
attr(discard.ercc2, "thresholds")["higher",]
# D10 D17 D2 D3 D7
# 9.475052 7.599947 6.010975 9.475052 15.216956
4 根据阈值筛选细胞
4.1 可视化阈值指标
主要用来进一步核实确定我们的指标是否合理，然后就可以剔除细胞了。
还是以上面的sce.416b数据集为例
代码语言：javascript复制sce.416b
sce.416b$block <- factor(sce.416b$block)
#注意下这个数据集只有block(20160113 20160325)里注释的两个批次，如下是实验分组信息。
sce.416b$phenotype <- ifelse(grepl("induced", sce.416b$phenotype),
"induced", "wild type")
sce.416b$discard <- reasons$discard
（1）library size等单独指标
代码语言：javascript复制plotColData(sce.416b, x="block", y="sum", colour_by="discard",
other_fields="phenotype") facet_wrap(~phenotype)
scale_y_log10() ggtitle("Total count")
如下图，橙色的点即为鉴定为不合格的细胞。把上述绘图代码里的sum换为detected、subsets_Mito_percent、altexps_ERCC_percent即可绘制对应的图。
（2）线粒体基因比例和library size的点图(剔除细胞是否聚集在左上角)
代码语言：javascript复制plotColData(sce.416b, x="sum", y="subsets_Mito_percent", colour_by="discard")
如下图因为细胞数少，所以聚集不太明显
（3）线粒体基因比例和spike-in 比例的点图(保留细胞是否聚集在左下角)
代码语言：javascript复制plotColData(sce.416b, x="altexps_ERCC_percent", y="subsets_Mito_percent", colour_by="discard")
4.2 核实剔除细胞是否包含一种特定的细胞类型
绘制剔除细胞与保留细胞的average count--logFoldchange的基因点图。
标注出我们感兴趣的细胞类型的marker gene是否显著的分布在上方(高表达)，如果有那么就有可能不小心误删某种细胞；如果没有那么就符合我们的预期。
代码语言：javascript复制# 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)
plot(abundance, logFC, xlab="Average count", ylab="Log-FC (lost/kept)", pch=16)
cell_type_marker_gene <- c("gene1", "gene2", "gene3",.....)
points(abundance[cell_type_marker_gene], logFC[cell_type_marker_gene], col="orange", pch=16)
举一个其它数据集的例子。该数据集，确定阈值指标后，发现血小板的3个marker gene在剔除细胞中高表达，表明可能剔除了血小板细胞。
4.3 剔除低质量细胞
在做完上述的步骤之后，就可以真正剔除细胞了。这一步还是非常简单的。一步到位。
代码语言：javascript复制# Keeping the columns we DON'T want to discard.
filtered <- sce.416b[,!reasons$discard]
在质控完成后，就是对原始count数据的标准化，会在下一节进行学习。如有问题，欢迎留言指出~~
往期回顾
单细胞分析十八般武艺4：velocyto
clustree—聚类可视化利器
肺的正常上皮细胞可以分成这5群
明码标价之10X转录组原始测序数据的cellranger流程
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
- 数据挖掘（GEO,TCGA,单细胞）2021第2期
- 生信爆款入门-2021第2期
- 96核心384G内存的超级服务器(共享)使用权一年
