3.1 质控

刘小泽写于19.12.3-4 ，更新于2020-06-26
1 质控必要性
前言
scRNA-seq数据的低质量文库可能来自于：细胞分选环节的破坏、文库制备失误（没有足够的反转录或PCR次数）... 表现在：细胞总表达量低、基本没有表达的基因、高线粒体或spike-in占比。
它带来的下游分析问题有：
聚类问题：这些低质量细胞肯定是最相似的了，“近朱者赤近墨者黑”，它们也会聚集成一群，对结果的解释造成干扰，因为从这群细胞中得不到什么有用的信息，但是它的确也是一群。这种现象产生的原因有可能是：细胞破坏以后，线粒体或核RNAs占比升高。 最差的情况就是：不同类型的低质量细胞，也能聚在一起，因为相比于固有的生物差异，更相似的低质量让它们相依相偎。除此以外，本来非常小的细胞文库也能聚成一群，因为log转换后它们的平均表达量会发生很大的变化(A. Lun 2018).
这篇文章举了个例子： 
​
另外还说：
Artifacts are most pronounced when the counts are low and there is large variation in the size factors across cells
Quality control is typically performed to remove cells with small library sizes. This reduces the variation in the size factors and the potential for artificial differences upon log-transformation
Low-abundance genes can be filtered out, which provides further protection as genes with low counts are most affected by the log-transformation
差异估算或主成分选择：首先在PCA分析时，低质量和高质量之间的差异相比于生物学差异会体现更明显，占据主要的成分，减少降维结果的可靠性。另外，某个基因可能在两个细胞之间没什么表达差异，但是如果所处环境差异很大（一个细胞质量很低，另一个细胞质量正常），那么在差异估算过程中，就会把这个差异也会被当成差异表达基因。例如：一个低质量细胞文库的总表达量非常低（接近0），但恰巧还存在一个基因有表达量，那么这个基因的表达量在后续的文库归一化过程中就会尤为突出
意外发现奇怪的转录本上调：实验难免会混入外源的污染转录本，但这个量很少并且在所有细胞中都应该是差不多水平的。但如果某个细胞质量低，文库小，那么在校正文库差异过程中（就像上面👆说的一样），其中的污染转录本表达量就会“突飞猛进”，看起来是一些明显上调的“奇怪基因”。实际上，这些奇怪的基因依然在其他细胞中存在，并且真实的表达量差不多，并且是不应该占据主体地位的
为了最大程度避免上面一种或多种情况的发生，应该在归一化之前去掉这些低质量的细胞，这个过程就是细胞的质控  （尽管对于droplet-based数据来讲，一个液滴/文库未必是一个细胞，因为存在没细胞dropout和双细胞doublet，甚至是多细胞的情况；但是这里为了方便介绍，不对”cell“和”library“名词进行区分）
使用的测试数据
下面将会使用A. T. L. Lun et al. (2017)的小型scRNA数据（未进行QC）
# BiocManager::install("scRNAseq")
library(scRNAseq)
sce.416b <- LunSpikeInData(which="416b")
需要注意的是，在下载数据的时候，它需要在/home下新建一个目录，问你yes/no 
​
如果你对/home目录没有操作权限，那么就会报错：
当然，如果不能新建也没有关系，我已经把这个数据替大家保存成了RData并放到了网盘，直接下载后加载就好
sce.416b链接：https://share.weiyun.com/5pQf0jg 密码：vab6wu
load('sce.416b.RData')
# 看到有两个96孔板的细胞
> table(sce.416b$block)
20160113 20160325 
      96       96 
​
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):
spikeNames(0):
altExpNames(2): ERCC SIRV
2 质控的指标
前言
鉴定细胞是否是低质量的，需要用到几个指标。虽然下面这些指标是使用SMART-seq2数据进行展示的，但依然适用于UMI数据（比如MARS-seq、droplet-based技术）
文库大小（library size） ：指的是每个细胞中所有基因的表达量之和。细胞的文库如果很小，说明在文库制备过程中存在RNA的损失，要么是由于细胞裂解，要么是cDNA捕获不足导致后续的扩增量少
每个细胞中表达基因的数目（number of expressed features in each cell）：指的是细胞中非0表达量的基因数量。如果细胞中基本没有基因表达，很可能是转录本捕获失败
比对到spike-in转录本的reads比例（proportion of reads mapped to spike-in transcripts）：计算方法是：spike-in counts / all features (including spike-ins) for each cell。每个细胞都应该加入等量的外源转录本（spike-in），如果哪个细胞的spike-in比例提高了，说明它的内源RNA含量减少（比如在细胞分选阶段出现的细胞裂解或者RNA降解）
比对到线粒体基因组的reads比例（proportion of reads mapped to genes in the mitochondrial genome） ：如果没有spike-in，那么使用线粒体指标也是能说明问题的(Islam et al. 2014; Ilicic et al. 2016)。比对到线粒体的reads增多，说明细胞质中的RNA减少，可能存在细胞穿孔的情况，而这个孔的大小，可能只是将细胞质中存在的mRNA流出去，但线粒体的体积超过了孔的大小，所以还留在了细胞中，造成一定程度的富集，导致线粒体RNA占比升高。
另外对于单核转录组测序（snRNA, single-nuclei RNA-seq），高线粒体占比也能反映细胞质是否被成功剥离。 下面是关于snRNA的介绍：
2018年，仅仅单核RNA-seq（snRNA-seq）就有160多篇论文发表。这是一种相对较新的方法，分析的是细胞核，而不是完整的细胞
在有些情况下，细胞是很难完整分离的，比如长期保存的组织，或者大脑细胞和脂肪细胞。在分离细胞时所用的酶和破坏力往往会影响其他细胞区室的内容。此时，单核RNA测序就显得特别重要。
snRNA-seq采用微流体技术，将带有条形码的微珠和单个细胞核一起装入微小的液滴内。这些液滴作为纳升级的反应管，实现了成千上万个细胞核分离和建库，从而大幅降低了建库的成本
单核RNA测序通过揭示复杂肿瘤的细胞类型、状态、遗传多样性和相互作用，拓宽了单细胞研究的能力，特别是长期保存或冷冻保存的样品。它克服了解离偏倚的问题，也适用于结构复杂的组织
对于每个细胞，可以用scater包的perCellQCMetrics()函数进行计算，其中sum这一列表示每个细胞的文库大小；detected这一列表示检测到的基因数量；subsets_Mito_percent这一列表示比对到线粒体基因组的reads占比；altexps_ERCC_percent表示比对到ERCC spike-in的reads占比
上手操作之统计线粒体信息
通过上面👆查看sce.416b发现，其中有spike-in的信息，但却没有线粒体相关的信息。这个没有关系，其实可以自己统计
第一种方法：
官网教程提供的是AnnotationHub()，这个方法需要下载80多M的数据库文件，对网速要求很高，动不动就失败。而且即使AnnotationHub()成功，后面的提取ah[["AH73905"]] 也很悬。不管怎样，把这个方法先放在这里，以供参考即可。
library(AnnotationHub)
ah <- AnnotationHub()
ens.mm.v97 <- ah[["AH73905"]]
location <- mapIds(ens.mm.v97, keys=rownames(sce.416b),
                   keytype="GENEID", column="SEQNAME")
is.mito <- which(location=="MT")
第二种方法：
library(TxDb.Mmusculus.UCSC.mm10.ensGene)
location <- mapIds(TxDb.Mmusculus.UCSC.mm10.ensGene, 
                   keys=rownames(sce.416b), 
                   keytype="GENEID",column="CDSCHROM")
is.mito <- which(location=="MT")
summary(location=="chrM")
> summary(location=="chrM")
   Mode   FALSE    TRUE    NA's 
logical   22428      13   24163 
# 看到只有13个线粒体基因
第三种方法：
如果可以调用rowRanges来获取基因组坐标数据的话，就可以结合seqnames找到MT
location <- rowRanges(sce.416b)
is.mito <- any(seqnames(location)=="MT")
上手操作之计算各种qc指标
library(scater)
df <- perCellQCMetrics(sce.416b, subsets=list(Mito=is.mito))
df
另外，还可以使用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"                 "percent_top_50"          
## [13] "percent_top_100"          "percent_top_200"         
## [15] "percent_top_500"          "subsets_Mito_sum"        
## [17] "subsets_Mito_detected"    "subsets_Mito_percent"    
## [19] "altexps_ERCC_sum"         "altexps_ERCC_detected"   
## [21] "altexps_ERCC_percent"     "altexps_SIRV_sum"        
## [23] "altexps_SIRV_detected"    "altexps_SIRV_percent"    
## [25] "total"
当然，这里做QC统计的前提假设是：每个细胞的qc指标都是独立于生物学状态的。也就是说，qc指标不会那么智能地识别细胞类型然后进行判断。它会认为（如文库太小、高线粒体占比）都是由技术误差引起的，而非生物因素。
那么问题来了，如果某些细胞类型本身的RNA含量就很低，或者它们本来就含有很多的线粒体转录本呢？再根据这个指标进行过滤，会不会损失一些细胞类型呢？ 【这个问题会在下面👇第4部分和第6部分进行说明】
3 鉴定低质量细胞
3.1 使用固定的阈值--简单粗暴
这是最简单粗暴的方法，例如设定文库低于10万reads的细胞是低质量的，或者表达基因数量少于5000个，spike-in或线粒体占比高于10%...
# 基于之前perCellQCMetrics计算的QC结果df
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
discard <- qc.lib | qc.nexprs | qc.spike | qc.mito
还可以统计一下每种过滤掉多少细胞
> 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         0        21
# 这里看到按照线粒体部分按照阈值为10%，是不能对细胞进行过滤的，因为我们根据TxDb.Mmusculus.UCSC.mm10.ensGene一共才找到了13个线粒体基因；而作者的教程中可能找到的线粒体基因比较多，因此他最终过滤掉了14个细胞
​
# 另外看到Total这里并不是简单地将四项指标相加，这是因为可能一个细胞同时符合两项或多项过滤条件，最终只计算一次
# 比如下面就是查看：既符合qc.lib过滤又符合qc.spike过滤条件的细胞=》第191个细胞
> intersect(which(qc.lib == 1),which(qc.spike == 1))
[1] 191
虽然看起来简单，但使用这种方法需要丰富的经验，了解实验设计和细胞状态；另外即使使用同一种文库制备方法，但由于细胞捕获效率和测序深度的不同，这个阈值依然需要适时调整
3.2 使用”相对“的阈值
3.2.1 鉴定离群点
先假设大部分细胞都是高质量的，然后去找离群点作为低质量。
那么按照什么方法找离群点呢？常用的一个函数isOutlier使用的是MAD指标（绝对中位差来估计方差,先计算出数据与它们的中位数之间的偏差，然后这些偏差的绝对值的中位数就是MAD，median absolute deviation）
函数认为超过中位数3倍MAD的值就是离群值
使用isOutlier时，如果要相减（例如：df$sum - 3* MAD），就用type="lower" ，此时一般还要做个log转换log=TRUE ，保证得到的结果不是负数
qc.lib2 <- isOutlier(df$sum, log=TRUE, type="lower")
qc.nexprs2 <- isOutlier(df$detected, log=TRUE, type="lower")
qc.spike2 <- isOutlier(df$altexps_ERCC_percent, type="higher")
qc.mito2 <- isOutlier(df$subsets_Mito_percent, type="higher")
看下得到的结果
> attr(qc.lib2, "thresholds")
   lower   higher 
434082.9      Inf 
> attr(qc.nexprs2, "thresholds")
   lower   higher 
5231.468      Inf 
> attr(qc.spike2, "thresholds")
   lower   higher 
    -Inf 14.15371 
> attr(qc.mito2, "thresholds")
   lower   higher 
    -Inf 6.472705
注意到：最后的线粒体部分，官方教程计算的结果是：
##  lower higher 
##   -Inf  11.92
用相对阈值过滤的细胞数量统计：
discard2 <- qc.lib2 | qc.nexprs2 | qc.spike2 | qc.mito2
​
> 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
最后过滤掉的细胞数量与官方教程一致，根据线粒体过滤的细胞都是2个。
除此以外，还有一种更快的计算方法，一步整合了上面的操作：
# 看帮助文档得知，需要提供额外的percent信息：
# quickPerCellQC(df, lib_size = "sum", n_features = "detected",
#   percent_subsets = NULL, ...)
reasons <- quickPerCellQC(df, percent_subsets=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
小结：使用”相对“的阈值一个好处就是可以根据测序深度、cDNA捕获效率、线粒体含量等等进行阈值的调整，在经验不是足够丰富的时候，可以辅助判断。但仍需要注意的是，使用相对阈值是有前提的，那就是：认为大部分细胞都是高质量的
3.2.2 离群点检测的假设
第一点：上面提到，离群点的检测的假设是大部分细胞的质量都不错。这一点假设可以通过实验去验证（比如肉眼检查细胞完整性）。但当大部分细胞质量都很低，使用相对阈值结果就对大量的低质量细胞无计可施。因为它使用了MAD值，和中位数有关系，那么可以试想：如果一堆数都不合格，那么它们的中位数也不合格，因此原来正确的值，其实在这群不合格的数值中，就是“离群”的，因此它只能从总体中剔除个别，而不能为了个别牺牲总体【简单理解：“谁人多势众相对阈值就听谁”】
另一个假设就是：每个细胞的qc指标都是独立于生物学状态的。也就是说，qc指标不会那么智能地识别细胞类型然后进行判断。在异质性很高的组织中， 就是存在内源RNA含量低，而线粒体基因占比高的细胞。如果使用”一刀切“的固定阈值，它们就很可能会被当成离群点被过滤。而是用MAD计算方法检测的结果可能就是：虽然一堆细胞的某个qc指标差异很大，但中位数也在变，随之变化的还有MAD值，这样最后的结果不至于和真实生物学情况差太多
因此，绝对阈值和相对阈值各有千秋，如果没有对”一刀切“方法有强烈的需求（比如当大部分细胞质量都很低时，就可以考虑一刀切），一般情况还是使用相对阈值为好
另外，这些假设知道即可， 不会影响后续的分析，也不用太多在意，先完成后完美
3.2.3 检测离群点还需要考虑实验设计（批次信息）
很多大型的实验都包含多个细胞系，而且可能采用的实验方法不同（比如选用不同的测序深度），这就产生了实验的不同批次。这种情况下， 应该对每个批次分别进行检测。而不要混在一起再去检测MAD，那样势必有一批会被”矫枉过正“，一批会”放任自流“。
如果每个批次都是一个SingleCellExperiment对象，那么isOutlier()函数可以按上面的方法直接使用；但是如果不同批次的细胞已经混合成一整个SingleCellExperiment对象，那么batch=这个参数就派上用场了。
还是以这个416B数据集为例，我们之前看到包含了两组实验信息
# 一个是实验批次，两个细胞板
> table(sce.416b$block)
20160113 20160325 
      96       96 
# 一个是细胞表型，两种生物状态
> table(sce.416b$phenotype)
induced CBFB-MYH11 oncogene expression  96            
wild type phenotype  96
而它们现在是混合在一起的，于是就可以设定batch信息
batch <- paste0(sce.416b$phenotype, "-", sce.416b$block)
# 得到2x2=4个批次信息
table(batch)
# induced CBFB-MYH11 oncogene expression-20160113 
# 48 
# induced CBFB-MYH11 oncogene expression-20160325 
# 48 
# wild type phenotype-20160113 
# 48 
# wild type phenotype-20160325 
# 48 
​
batch.reasons <- quickPerCellQC(df, percent_subsets=c("subsets_Mito_percent",
                                                      "altexps_ERCC_percent"), batch=batch)
​
> colSums(as.matrix(batch.reasons))
##              low_lib_size            low_n_features high_subsets_Mito_percent 
##                         4                         2                         2 
## high_altexps_ERCC_percent                   discard 
##                         4                         7
但是，batch参数不是万能的，之前也说过，这种相对阈值需要一个假设：每个批次的大部分细胞都是高质量的。当某个批次的细胞整体质量偏低，离群点检测很有可能失败。举个例子，Grun et al. (2016)做了一个人类胰腺癌的数据集，里面总共有5个donor的数据，但事实上有两个donor的实验是失败的。它们的ERCC含量相对其他批次高，增加了中位数和MAD值，导致过滤低质量细胞失败。因此这种情况下，可以先算其他几个批次的中位数和MAD值，然后参考这些值去对有问题的批次进行过滤
library(scRNAseq)
sce.grun <- GrunPancreasData()
# 批次信息
> table(sce.grun$donor)
D10 D17  D2  D3  D7 
288 480  96 480 384
​
sce.grun <- addPerCellQC(sce.grun)
# 分批次进行检测，看看效果
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))
with.blocking
图中每个点表示一个细胞。可以看到，D10和D3的检测是失败的，因为它们的大部分细胞ERCC含量本身就很高，所以最后并没有检测到太多的离群点（黄色点）
这时，就可以先算其他几个批次的离群点：
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)
可以看到，左图是按照每个批次分别鉴定的离群点；右图是用质量好的批次计算的阈值，然后运用到差的批次上，群策群力得到的离群点
如何鉴别有问题的批次呢？
可以先将每个批次分别计算结果，然后比较它们的阈值，如果比同类批次超出太多，就要警觉
ercc.thresholds <- attr(discard.ercc, "thresholds")["higher",]
ercc.thresholds
##        D10        D17         D2         D3         D7 
##  73.610696   7.599947   6.010975 113.105828  15.216956
从数值就能看到D10、D3的阈值就超过其他批次太多
names(ercc.thresholds)[isOutlier(ercc.thresholds, type="higher")]
## [1] "D10" "D3"
但还是那句话，这么做的前提都是：我们认为批次中的细胞质量整体还不错。
但如果我们保证不了细胞质量，那么这种计算相对阈值的方法就不成立，还是要使用绝对阈值，手动去除
3.3 另辟蹊径的检测方法
利用robustbase 包，将细胞放在高维空间，根据他们的QC指标（文库大小、表达基因数、spike-in比例、线粒体比例）
它的用法是：For an n * p data matrix (or data frame) x, compute the “outlyingness” of all n observations
先做一个行是细胞，列是QC指标的数据框
stats <- cbind(log10(df$sum), log10(df$detected),
               df$subsets_Mito_percent, df$altexps_ERCC_percent)
> dim(stats)
[1] 192   4
然后j就像做PCA分析一样，剔除共性，找出最有差异的部分
library(robustbase)
outlying <- adjOutlyingness(stats, only.outlyingness = TRUE)
multi.outlier <- isOutlier(outlying, type = "higher")
​
> attr(multi.outlier, "thresholds")
   lower   higher 
    -Inf 2.146625 
​
> summary(multi.outlier)
   Mode   FALSE    TRUE 
logical     182      10
3.4 其他
除此以外，有时还可以根据基因表达量进行过滤，不过在bulk转录组中用的比较多，但是在scRNA中这样操作可能会损失一些本身数量就比较少的高质量细胞群体（比如一些罕见细胞，本身基因表达量就不是很高）
4 画图检查
就是对QC指标的分布进行可视化，来检查是否存在问题
第一张图 不同批次的QC指标
首先 把QC指标和原来的sce.416b细胞信息整合起来
> dim(colData(sce.416b))
[1] 192   9
> dim(df)
[1] 192  16
​
colData(sce.416b) <- cbind(colData(sce.416b), df)
> names(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"                 "percent_top_50"          
[13] "percent_top_100"          "percent_top_200"         
[15] "percent_top_500"          "subsets_Mito_sum"        
[17] "subsets_Mito_detected"    "subsets_Mito_percent"    
[19] "altexps_ERCC_sum"         "altexps_ERCC_detected"   
[21] "altexps_ERCC_percent"     "altexps_SIRV_sum"        
[23] "altexps_SIRV_detected"    "altexps_SIRV_percent"    
[25] "total"
修改一下整合后的信息
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=2
)
第二张图 线粒体占比与文库大小
为了确保不存在这样的细胞：虽然细胞文库大，但它的线粒体占比也高。另外也是为了避免意外过滤掉本来是高质量但同时具有高代谢活性（即高线粒体占比）的细胞（如肝脏细胞）
主要关注右上角（代表大文库，高线粒体含量），如果真的存在，就要看与这些细胞的生物类型是不是相符
plotColData(sce.416b, x="sum", y="subsets_Mito_percent", 
    colour_by="discard", other_fields=c("block", "phenotype")) +
    facet_grid(block~phenotype) +
    theme(panel.border = element_rect(color = "grey"))
第三张图 比较ERCC和线粒体比例
有时会存在这样的情况：细胞质量低，文库大小低，它的线粒体占比也低，但是spike-in占比很高
它可能是因为细胞破坏非常严重，导致细胞质组分基本丢失，最后建库只捕获了外源的RNA
而另一种相反的情况：高线粒体占比，低spike-in占比，则有可能是说明细胞没有损伤，只是它的代谢活性很强。
这种分析思路同样适用于单核RNA-Seq（snRNA-Seq），但关注点从细胞文库质量转移到了细胞核文库质量
plotColData(sce.416b, x="altexps_ERCC_percent", y="subsets_Mito_percent",
    colour_by="discard", other_fields=c("block", "phenotype")) +
    facet_grid(block~phenotype) + 
    theme(panel.border = element_rect(color = "grey"))
如果在右下角存在许多点时就要注意了，可能细胞受到损伤；而左上角的点还需要进一步结合细胞生物学类型进行判断
5 针对Droplet数据的细胞判断（例如10X）
5.1 前言
由于这种建库方法的独特性，我们没办法事先知道某一个细胞文库（比如一个cell barcode）是真正包含一个细胞还是只是一个空的液滴（droplet）。因此，第一步是需要根据得到的表达量信息，来推断液滴是不是空的。但仅仅根据表达量判断还是不太靠谱（是不是很复杂的问题？），因为还有可能在空的液滴中依然包含外源RNA，最后的细胞文库依旧不为0，但确实不包含细胞。
这里为了说明这个问题，使用了一个未过滤的10X PBMC数据
# 学习一下数据下载的方式（新建一个目录，然后把链接放进函数进行下载）
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
raw.path <- bfcrpath(bfc, file.path("http://cf.10xgenomics.com/samples",
                                    "cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"))
# 把数据解压到当前目录
untar(raw.path, exdir=file.path(gwtwd(), "pbmc4k"))
​
library(DropletUtils)
library(Matrix)
fname <- file.path(gwtwd(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
当然，如果你想省去下载这一步，还可以用我保存的数据【sce.pbmc的网盘链接：链接：https://share.weiyun.com/5LRF8Jn 密码：sz7wmv】
load("sce.pbmc.Rdata")
# 298M的对象
> sce.pbmc
class: SingleCellExperiment 
dim: 33694 737280 
metadata(1): Samples
assays(1): counts
rownames(33694): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
  ENSG00000268674
rowData names(2): ID Symbol
colnames(737280): AAACCTGAGAAACCAT-1 AAACCTGAGAAACCGC-1 ... TTTGTCATCTTTAGTC-1
  TTTGTCATCTTTCCTC-1
colData names(2): Sample Barcode
reducedDimNames(0):
spikeNames(0):
altExpNames(0):
看看不同的barcodes（不一定都是真实的细胞）的文库大小分布：
bcrank <- barcodeRanks(counts(sce.pbmc))
# 为了作图速度，每个rank只展示一个点（去重复）
uniq <- !duplicated(bcrank$rank)
plot(bcrank$rank[uniq], bcrank$total[uniq], log="xy",
     xlab="Rank", ylab="Total UMI count", cex.lab=1.2)
​
abline(h=metadata(bcrank)$inflection, col="darkgreen", lty=2)
abline(h=metadata(bcrank)$knee, col="dodgerblue", lty=2)
​
legend("bottomleft", legend=c("Inflection", "Knee"), 
       col=c("darkgreen", "dodgerblue"), lty=2, cex=1.2)
看到各个barcodes的文库大小有高有低，并且相差较大，因此可能对应着真实存在的细胞和空液滴。当然最简单的办法还是”一刀切“，留下那些文库较大的细胞，不过还是有损失真实细胞类型的风险。
5.2 【作了解】检测空的液滴
主要使用emptyDrops()函数，检查每个barcode的表达量是不是和外源RNA的表达量有显著差异(Lun et al. 2019)。如果存在显著差异，就说明barcode中更有可能是一个细胞。这种方法可以帮助区分：测序质量好的空液滴 和 包含细胞但RNA含量较少的液滴。尽管它们总体的表达量可能很相似，但本质不同，还是要区分的。
这个函数假设总体UMI含量低就是空的液滴，使用Monte Carlo统计模拟计算p值，如果要重复结果，需要设置随机种子。另外设置 false discovery rate (FDR)为0.1%，意味着不超过0.1%的barcodes是空的
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
> summary(e.out$FDR <= 0.001)
   Mode   FALSE    TRUE    NA's 
logical    1056    4233  731991
将最后的结果保存：
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]
​
> ncol(counts(sce.pbmc))
[1] 737280
> length(which(e.out$FDR <= 0.001))
[1] 4233
可以看到从73万个液滴中，最后只得到了4000多个带细胞的液滴
这个过程，自己很有可能是用不到的，CellRanger V3提供了一套检测细胞的算法，和emptyDrops的功能相似。因此我们一般读进来直接是几千个细胞（CellRanger处理之后的），而不会是几十万的液滴（包含空液滴）。当然，如果用已过滤的数据再去检测空细胞，结果一样是不准确的。
5.3 与其他质控指标的结合
虽然上一步emptyDrops检测了空的液滴，但是依然不能知道有细胞的液滴中细胞质量如何。因为即使包含细胞，还是有可能是死细胞或损伤的细胞。
假设我们知道这些细胞中没有代谢活性很高的细胞（也就是说可以根据线粒体占比去除细胞）：
is.mito <- grep("^MT-", rowData(sce.pbmc)$Symbol)
pbmc.qc <- perCellQCMetrics(sce.pbmc, subsets=list(MT=is.mito))
discard.mito <- isOutlier(pbmc.qc$subsets_MT_percent, type="higher")
> summary(discard.mito)
   Mode   FALSE    TRUE 
logical    3922     311
作图检查：
plot(pbmc.qc$sum, pbmc.qc$subsets_MT_percent, log="x",
    xlab="Total count", ylab='Mitochondrial %')
abline(h=attr(discard.mito, "thresholds")["higher"], col="red")
最后的过滤依然是需要自己对数据的把握：
一方面考虑这些筛选出来的”不太合格“的细胞对下游分析的影响到底有多大；
二是考虑如果去掉这些细胞，会不会损失一些细胞类型？
6 去掉低质量细胞
一旦细胞质控完成，可以选择是直接去掉这些细胞，还是先标记出来给个”缓刑“的机会
直接去除很简单，直接对sce对象的列取子集即可，例如：
filtered <- sce.416b[,!reasons$discard]
不过，质控过程的最大难题是不小心丢掉了某些细胞类型，毕竟质控指标总是和生物状态有着千丝万缕的关系。
好，下面先给那些不符合要求的细胞一个”缓刑“的机会：
计算了舍弃组和保留组细胞的平均表达量以及舍弃组与保留组之间的logFC
# 首先计算每个feature的平均表达量
lost <- calculateAverage(counts(sce.416b)[,!discard])
kept <- calculateAverage(counts(sce.416b)[,discard])
​
library(edgeR)
# 对lost和kept分别计算log(cpm + 2)
logged <- cpm(cbind(lost, kept), log=TRUE, prior.count=2)
​
> head(logged)
                       lost     kept
ENSMUSG00000102693 1.013482 1.013482
ENSMUSG00000064842 1.013482 1.013482
ENSMUSG00000051951 1.013482 1.013482
ENSMUSG00000102851 1.013482 1.013482
ENSMUSG00000103377 1.019356 1.013482
ENSMUSG00000104017 1.013482 1.013482
​
abundance <- rowMeans(logged)
logFC <- logged[,1] - logged[,2]
做个图：
plot(abundance, logFC, xlab="Average count", ylab="Log-FC (lost/kept)", pch=16)
points(abundance[is.mito], logFC[is.mito], col="dodgerblue", pch=16)
如果说目前被丢弃的细胞真的是某个特别的细胞类型，那么这个细胞类型应该有属于自己的marker基因，那么这些基因应该在被丢弃的这群中显著高表达。看到图中纵轴0以上应该是舍弃组比保留组中高表达的基因，但是呢，这部分基因在X轴上展示的平均表达量又没有特别高（主要在2-10之间）
只看一个例子还不够，再看一个反例：
假如这里用固定阈值随意去除细胞（那么肯定有真实的细胞类型被去除）
# 就是这么随意去除
alt.discard <- colSums(counts(sce.pbmc)) < 500
# 下面再次用👆的计算方法
lost <- calculateAverage(counts(sce.pbmc)[,alt.discard])
kept <- calculateAverage(counts(sce.pbmc)[,!alt.discard])
​
logged <- edgeR::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)
发现这时就会有很多基因在丢弃组中上调表达，且表达量还不低
根据先验知识，推测有一些血小板相关基因在丢弃组中高表达（例如：PF4, PPBP and CAVIN2），进行可视化
platelet <- c("PF4", "PPBP", "CAVIN2")
library(org.Hs.eg.db)
ids <- mapIds(org.Hs.eg.db, keys=platelet, column="ENSEMBL", keytype="SYMBOL")
points(abundance[ids], logFC[ids], col="orange", pch=16)
看图片，果然血小板相关的基因在丢弃组中高表达
如果感觉自己的过滤太过于严格，而担心损失一些细胞，那么可以尝试调高isOutlier()的nmads=参数
7 【看看就好】如果不去除，只是标记呢？
低质量细胞还可以只是标记出来，然后继续进行下游分析，这样为了防止自己随便过滤而造成细胞类型丢失。
这样其实可以在后面的聚类分析中，能看到标记出的这部分细胞可以单独聚在一起（而我们心里也会清楚，它们为何会聚在一起）
如何进行标记？
marked <- sce.416b
marked$discard <- batch.reasons$discard
虽然这样做很保险，但会增加质控结果的解释难度。当QC解读起来太难，可能又会到后期根据marker基因去分辨低质量和正常的细胞。另外会增加后续的计算量。
总结
作者依然是建议提前去除的，长痛不如短痛，不要把低质量的信息加到后面更加复杂的计算中。
如果对去除的标准没有信心，可以在走完一遍常规流程后，从头回来再去进行标记，走一下另一条路，看看和常规的结果有什么不同，是不是真的多了一些”容易忽略“的细胞类型？