【原】OSCA单细胞数据分析笔记-5 Quality control

对应原版教程第6章  http:///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格式（小鼠），所以需要注释下基因所处的染色体序列信息。

# 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

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数据集为例

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等单独指标

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的点图(剔除细胞是否聚集在左上角)

plotColData(sce.416b, x="sum", y="subsets_Mito_percent", colour_by="discard")


如下图因为细胞数少，所以聚集不太明显

（3）线粒体基因比例和spike-in 比例的点图(保留细胞是否聚集在左下角)

plotColData(sce.416b, x="altexps_ERCC_percent", y="subsets_Mito_percent", colour_by="discard")


4.2 核实剔除细胞是否包含一种特定的细胞类型

绘制剔除细胞与保留细胞的average count--logFoldchange的基因点图。

标注出我们感兴趣的细胞类型的marker gene是否显著的分布在上方(高表达)，如果有那么就有可能不小心误删某种细胞；如果没有那么就符合我们的预期。

# 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 剔除低质量细胞

在做完上述的步骤之后，就可以真正剔除细胞了。这一步还是非常简单的。一步到位。

# Keeping the columns we DON'T want to discard.
filtered <- sce.416b[,!reasons$discard]


在质控完成后，就是对原始count数据的标准化，会在下一节进行学习。如有问题，欢迎留言指出~~