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

对应原版教程第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格式(小鼠),所以需要注释下基因所处的染色体序列信息。

# 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换为detectedsubsets_Mito_percentaltexps_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数据的标准化,会在下一节进行学习。如有问题,欢迎留言指出~~

(0)

相关推荐