單細胞交響樂3-scRNA的細胞質(zhì)控

劉小澤寫于19.12.3-4 铃将，更新于2020-06-26
為何取名叫“交響樂”择懂？因為單細胞分析就像一個大樂團运敢，需要各個流程的協(xié)同配合

單細胞交響樂1-理解scRNA常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細胞交響樂2-scRNA從實驗到下游簡介

1 質(zhì)控必要性

前言

scRNA-seq數(shù)據(jù)的低質(zhì)量文庫可能來自于：細胞分選環(huán)節(jié)的破壞、文庫制備失誤（沒有足夠的反轉(zhuǎn)錄或PCR次數(shù)）... 表現(xiàn)在：細胞總表達量低穆端、基本沒有表達的基因袱贮、高線粒體或spike-in占比。

它帶來的下游分析問題有：

聚類問題：這些低質(zhì)量細胞肯定是最相似的了体啰，“近朱者赤近墨者黑”攒巍，它們也會聚集成一群，對結(jié)果的解釋造成干擾荒勇，因為從這群細胞中得不到什么有用的信息柒莉，但是它的確也是一群。這種現(xiàn)象產(chǎn)生的原因有可能是：細胞破壞以后沽翔，線粒體或核RNAs占比升高兢孝。
最差的情況就是：不同類型的低質(zhì)量細胞，也能聚在一起仅偎，因為相比于固有的生物差異跨蟹，更相似的低質(zhì)量讓它們相依相偎。除此以外橘沥，本來非常小的細胞文庫也能聚成一群窗轩，因為log轉(zhuǎn)換后它們的平均表達量會發(fā)生很大的變化(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分析時，低質(zhì)量和高質(zhì)量之間的差異相比于生物學(xué)差異會體現(xiàn)更明顯威恼，占據(jù)主要的成分品姓，減少降維結(jié)果的可靠性寝并。另外箫措，某個基因可能在兩個細胞之間沒什么表達差異，但是如果所處環(huán)境差異很大（一個細胞質(zhì)量很低衬潦，另一個細胞質(zhì)量正常）斤蔓，那么在差異估算過程中，就會把這個差異也會被當成差異表達基因镀岛。例如：一個低質(zhì)量細胞文庫的總表達量非常低（接近0）弦牡，但恰巧還存在一個基因有表達量友驮，那么這個基因的表達量在后續(xù)的文庫歸一化過程中就會尤為突出
意外發(fā)現(xiàn)奇怪的轉(zhuǎn)錄本上調(diào)：實驗難免會混入外源的污染轉(zhuǎn)錄本，但這個量很少并且在所有細胞中都應(yīng)該是差不多水平的驾锰。但如果某個細胞質(zhì)量低卸留，文庫小，那么在校正文庫差異過程中（就像上面??說的一樣）椭豫，其中的污染轉(zhuǎn)錄本表達量就會“突飛猛進”耻瑟，看起來是一些明顯上調(diào)的“奇怪基因”。實際上赏酥，這些奇怪的基因依然在其他細胞中存在喳整，并且真實的表達量差不多，并且是不應(yīng)該占據(jù)主體地位的

為了最大程度避免上面一種或多種情況的發(fā)生裸扶，應(yīng)該在歸一化之前去掉這些低質(zhì)量的細胞框都，這個過程就是**細胞的質(zhì)控 ** （盡管對于droplet-based數(shù)據(jù)來講，一個液滴/文庫未必是一個細胞呵晨，因為存在沒細胞dropout和雙細胞doublet魏保，甚至是多細胞的情況；但是這里為了方便介紹何荚，不對”cell“和”library“名詞進行區(qū)分）

使用的測試數(shù)據(jù)

下面將會使用A. T. L. Lun et al. (2017)的小型scRNA數(shù)據(jù)（未進行QC）

# BiocManager::install("scRNAseq")
library(scRNAseq)
sce.416b <- LunSpikeInData(which="416b")

需要注意的是囱淋，在下載數(shù)據(jù)的時候，它需要在/home下新建一個目錄餐塘，問你yes/no

如果你對/home目錄沒有操作權(quán)限妥衣，那么就會報錯：

當然，如果不能新建也沒有關(guān)系戒傻，我已經(jīng)把這個數(shù)據(jù)替大家保存成了RData并放到了網(wǎng)盤税手，直接下載后加載就好

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 質(zhì)控的指標

前言

鑒定細胞是否是低質(zhì)量的，需要用到幾個指標需纳。雖然下面這些指標是使用SMART-seq2數(shù)據(jù)進行展示的芦倒，但依然適用于UMI數(shù)據(jù)（比如MARS-seq、droplet-based技術(shù)）

文庫大胁霍妗（library size） ：指的是每個細胞中所有基因的表達量之和兵扬。細胞的文庫如果很小，說明在文庫制備過程中存在RNA的損失口蝠，要么是由于細胞裂解器钟，要么是cDNA捕獲不足導(dǎo)致后續(xù)的擴增量少
每個細胞中表達基因的數(shù)目（number of expressed features in each cell）：指的是細胞中非0表達量的基因數(shù)量。如果細胞中基本沒有基因表達妙蔗，很可能是轉(zhuǎn)錄本捕獲失敗
比對到spike-in轉(zhuǎn)錄本的reads比例（proportion of reads mapped to spike-in transcripts）：計算方法是：spike-in counts / all features (including spike-ins) for each cell傲霸。每個細胞都應(yīng)該加入等量的外源轉(zhuǎn)錄本（spike-in），如果哪個細胞的spike-in比例提高了，說明它的內(nèi)源RNA含量減少（比如在細胞分選階段出現(xiàn)的細胞裂解或者RNA降解）
比對到線粒體基因組的reads比例（proportion of reads mapped to genes in the mitochondrial genome） ：如果沒有spike-in昙啄，那么使用線粒體指標也是能說明問題的(Islam et al. 2014; Ilicic et al. 2016)穆役。比對到線粒體的reads增多，說明細胞質(zhì)中的RNA減少梳凛，可能存在細胞穿孔的情況耿币，而這個孔的大小，可能只是將細胞質(zhì)中存在的mRNA流出去韧拒，但線粒體的體積超過了孔的大小掰读，所以還留在了細胞中，造成一定程度的富集叭莫，導(dǎo)致線粒體RNA占比升高蹈集。
另外對于單核轉(zhuǎn)錄組測序（snRNA, single-nuclei RNA-seq），高線粒體占比也能反映細胞質(zhì)是否被成功剝離雇初。
下面關(guān)于snRNA的介紹來自：http://www.ebiotrade.com/newsf/2019-12/2019129172619642.htm
- 2018年拢肆，僅僅單核RNA-seq（snRNA-seq）就有160多篇論文發(fā)表。這是一種相對較新的方法靖诗，分析的是細胞核郭怪，而不是完整的細胞
- 在有些情況下，細胞是很難完整分離的刊橘，比如長期保存的組織鄙才，或者大腦細胞和脂肪細胞。在分離細胞時所用的酶和破壞力往往會影響其他細胞區(qū)室的內(nèi)容促绵。此時攒庵，單核RNA測序就顯得特別重要。
- snRNA-seq采用微流體技術(shù)败晴，將帶有條形碼的微珠和單個細胞核一起裝入微小的液滴內(nèi)浓冒。這些液滴作為納升級的反應(yīng)管，實現(xiàn)了成千上萬個細胞核分離和建庫尖坤，從而大幅降低了建庫的成本
- 單核RNA測序通過揭示復(fù)雜腫瘤的細胞類型稳懒、狀態(tài)、遺傳多樣性和相互作用慢味，拓寬了單細胞研究的能力场梆，特別是長期保存或冷凍保存的樣品。它克服了解離偏倚的問題纯路，也適用于結(jié)構(gòu)復(fù)雜的組織

對于每個細胞或油，可以用scater包的perCellQCMetrics()函數(shù)進行計算，其中sum這一列表示每個細胞的文庫大懈兄纭装哆；detected這一列表示檢測到的基因數(shù)量罐脊；subsets_Mito_percent這一列表示比對到線粒體基因組的reads占比定嗓；altexps_ERCC_percent表示比對到ERCC spike-in的reads占比

上手操作之統(tǒng)計線粒體信息

通過上面??查看sce.416b發(fā)現(xiàn)蜕琴，其中有spike-in的信息，但卻沒有線粒體相關(guān)的信息宵溅。這個沒有關(guān)系凌简，其實可以自己統(tǒng)計

第一種方法：

官網(wǎng)教程提供的是AnnotationHub()，這個方法需要下載80多M的數(shù)據(jù)庫文件恃逻，對網(wǎng)速要求很高，動不動就失敗。而且即使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個線粒體基因

第三種方法：

如果可以調(diào)用rowRanges來獲取基因組坐標數(shù)據(jù)的話，就可以結(jié)合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中，方便后面調(diào)取

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統(tǒng)計的前提假設(shè)是：每個細胞的qc指標都是獨立于生物學(xué)狀態(tài)的歌憨。也就是說，qc指標不會那么智能地識別細胞類型然后進行判斷墩衙。它會認為（如文庫太小务嫡、高線粒體占比）都是由技術(shù)誤差引起的，而非生物因素漆改。

那么問題來了植袍，如果某些細胞類型本身的RNA含量就很低，或者它們本來就含有很多的線粒體轉(zhuǎn)錄本呢籽懦？再根據(jù)這個指標進行過濾于个，會不會損失一些細胞類型呢？
【這個問題會在下面??第4部分和第6部分進行說明】

3 鑒定低質(zhì)量細胞

3.1 使用固定的閾值--簡單粗暴

這是最簡單粗暴的方法暮顺，例如設(shè)定文庫低于10萬reads的細胞是低質(zhì)量的厅篓，或者表達基因數(shù)量少于5000個，spike-in或線粒體占比高于10%...

# 基于之前perCellQCMetrics計算的QC結(jié)果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
# 都設(shè)定好以后傳給discard
discard <- qc.lib | qc.nexprs | qc.spike | qc.mito

還可以統(tǒng)計一下每種過濾掉多少細胞

> 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%捶码，是不能對細胞進行過濾的羽氮，因為我們根據(jù)TxDb.Mmusculus.UCSC.mm10.ensGene一共才找到了13個線粒體基因；而作者的教程中可能找到的線粒體基因比較多惫恼，因此他最終過濾掉了14個細胞

# 另外看到Total這里并不是簡單地將四項指標相加档押，這是因為可能一個細胞同時符合兩項或多項過濾條件，最終只計算一次
# 比如下面就是查看：既符合qc.lib過濾又符合qc.spike過濾條件的細胞=》第191個細胞
> intersect(which(qc.lib == 1),which(qc.spike == 1))
[1] 191

雖然看起來簡單，但使用這種方法需要豐富的經(jīng)驗令宿，了解實驗設(shè)計和細胞狀態(tài)叼耙；另外即使使用同一種文庫制備方法，但由于細胞捕獲效率和測序深度的不同粒没，這個閾值依然需要適時調(diào)整

3.2 使用”相對“的閾值

3.2.1 鑒定離群點

先假設(shè)大部分細胞都是高質(zhì)量的筛婉，然后去找離群點作為低質(zhì)量。

那么按照什么方法找離群點呢癞松？常用的一個函數(shù)isOutlier使用的是MAD指標（絕對中位差來估計方差,先計算出數(shù)據(jù)與它們的中位數(shù)之間的偏差爽撒，然后這些偏差的絕對值的中位數(shù)就是MAD，median absolute deviation）

函數(shù)認為超過中位數(shù)3倍MAD的值就是離群值

使用isOutlier時响蓉，如果要相減（例如：df$sum - 3* MAD）硕勿，就用type="lower" ，此時一般還要做個log轉(zhuǎn)換log=TRUE 枫甲，保證得到的結(jié)果不是負數(shù)

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")

看下得到的結(jié)果

> 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

注意到：最后的線粒體部分首尼，官方教程計算的結(jié)果是：
##  lower higher 
##   -Inf  11.92

用相對閾值過濾的細胞數(shù)量統(tǒng)計：

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

最后過濾掉的細胞數(shù)量與官方教程一致，根據(jù)線粒體過濾的細胞都是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

小結(jié)：使用”相對“的閾值一個好處就是可以根據(jù)測序深度查排、cDNA捕獲效率、線粒體含量等等進行閾值的調(diào)整抄沮，在經(jīng)驗不是足夠豐富的時候跋核，可以輔助判斷。但仍需要注意的是叛买，使用相對閾值是有前提的砂代，那就是：認為大部分細胞都是高質(zhì)量的

3.2.2 離群點檢測的假設(shè)

第一點：上面提到，離群點的檢測的假設(shè)是大部分細胞的質(zhì)量都不錯率挣。這一點假設(shè)可以通過實驗去驗證（比如肉眼檢查細胞完整性）刻伊。但當大部分細胞質(zhì)量都很低，使用相對閾值結(jié)果就對大量的低質(zhì)量細胞無計可施椒功。因為它使用了MAD值捶箱，和中位數(shù)有關(guān)系，那么可以試想：如果一堆數(shù)都不合格动漾，那么它們的中位數(shù)也不合格丁屎，因此原來正確的值，其實在這群不合格的數(shù)值中旱眯，就是“離群”的晨川，因此它只能從總體中剔除個別证九，而不能為了個別犧牲總體【簡單理解：“誰人多勢眾相對閾值就聽誰”】
另一個假設(shè)就是：每個細胞的qc指標都是獨立于生物學(xué)狀態(tài)的。也就是說共虑，qc指標不會那么智能地識別細胞類型然后進行判斷愧怜。在異質(zhì)性很高的組織中，就是存在內(nèi)源RNA含量低看蚜，而線粒體基因占比高的細胞。如果使用”一刀切“的固定閾值赔桌，它們就很可能會被當成離群點被過濾供炎。而是用MAD計算方法檢測的結(jié)果可能就是：雖然一堆細胞的某個qc指標差異很大，但中位數(shù)也在變疾党，隨之變化的還有MAD值音诫，這樣最后的結(jié)果不至于和真實生物學(xué)情況差太多

因此，絕對閾值和相對閾值各有千秋雪位，如果沒有對”一刀切“方法有強烈的需求（比如當大部分細胞質(zhì)量都很低時竭钝，就可以考慮一刀切），一般情況還是使用相對閾值為好

另外雹洗，這些假設(shè)知道即可香罐，不會影響后續(xù)的分析，也不用太多在意时肿，先完成后完美

3.2.3 檢測離群點還需要考慮實驗設(shè)計（批次信息）

很多大型的實驗都包含多個細胞系庇茫，而且可能采用的實驗方法不同（比如選用不同的測序深度），這就產(chǎn)生了實驗的不同批次螃成。這種情況下旦签，應(yīng)該對每個批次分別進行檢測。而不要混在一起再去檢測MAD寸宏，那樣勢必有一批會被”矯枉過正“宁炫，一批會”放任自流“。

如果每個批次都是一個SingleCellExperiment對象氮凝，那么isOutlier()函數(shù)可以按上面的方法直接使用羔巢；但是如果不同批次的細胞已經(jīng)混合成一整個SingleCellExperiment對象，那么batch=這個參數(shù)就派上用場了罩阵。

還是以這個416B數(shù)據(jù)集為例朵纷，我們之前看到包含了兩組實驗信息

# 一個是實驗批次，兩個細胞板
> table(sce.416b$block)
20160113 20160325 
      96       96 
# 一個是細胞表型永脓，兩種生物狀態(tài)
> table(sce.416b$phenotype)
induced CBFB-MYH11 oncogene expression  96            
wild type phenotype  96

而它們現(xiàn)在是混合在一起的袍辞，于是就可以設(shè)定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參數(shù)不是萬能的常摧，之前也說過搅吁，這種相對閾值需要一個假設(shè)：每個批次的大部分細胞都是高質(zhì)量的威创。當某個批次的細胞整體質(zhì)量偏低，離群點檢測很有可能失敗谎懦。舉個例子肚豺，Grun et al. (2016)做了一個人類胰腺癌的數(shù)據(jù)集，里面總共有5個donor的數(shù)據(jù)界拦，但事實上有兩個donor的實驗是失敗的吸申。它們的ERCC含量相對其他批次高，增加了中位數(shù)和MAD值享甸，導(dǎo)致過濾低質(zhì)量細胞失敗截碴。因此這種情況下，可以先算其他幾個批次的中位數(shù)和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)

可以看到，左圖是按照每個批次分別鑒定的離群點栅盲；右圖是用質(zhì)量好的批次計算的閾值湘今，然后運用到差的批次上，群策群力得到的離群點

如何鑒別有問題的批次呢剪菱？

可以先將每個批次分別計算結(jié)果摩瞎，然后比較它們的閾值，如果比同類批次超出太多孝常，就要警覺

ercc.thresholds <- attr(discard.ercc, "thresholds")["higher",]
ercc.thresholds
##        D10        D17         D2         D3         D7 
##  73.610696   7.599947   6.010975 113.105828  15.216956

從數(shù)值就能看到D10旗们、D3的閾值就超過其他批次太多

names(ercc.thresholds)[isOutlier(ercc.thresholds, type="higher")]
## [1] "D10" "D3"

但還是那句話，這么做的前提都是：我們認為批次中的細胞質(zhì)量整體還不錯构灸。

但如果我們保證不了細胞質(zhì)量上渴，那么這種計算相對閾值的方法就不成立，還是要使用絕對閾值喜颁，手動去除

3.3 另辟蹊徑的檢測方法

利用robustbase 包稠氮，將細胞放在高維空間，根據(jù)他們的QC指標（文庫大小半开、表達基因數(shù)隔披、spike-in比例、線粒體比例）

它的用法是：For an n * p data matrix (or data frame) x, compute the “outlyingness” of all n observations

先做一個行是細胞寂拆，列是QC指標的數(shù)據(jù)框

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 其他

除此以外，有時還可以根據(jù)基因表達量進行過濾鬓长，不過在bulk轉(zhuǎn)錄組中用的比較多谒拴，但是在scRNA中這樣操作可能會損失一些本身數(shù)量就比較少的高質(zhì)量細胞群體（比如一些罕見細胞，本身基因表達量就不是很高）

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")
# 增加之前過濾的結(jié)果
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
)

第二張圖線粒體占比與文庫大小

為了確保不存在這樣的細胞：雖然細胞文庫大英上，但它的線粒體占比也高。另外也是為了避免意外過濾掉本來是高質(zhì)量但同時具有高代謝活性（即高線粒體占比）的細胞（如肝臟細胞）

主要關(guān)注右上角（代表大文庫啤覆，高線粒體含量）苍日，如果真的存在，就要看與這些細胞的生物類型是不是相符

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和線粒體比例

有時會存在這樣的情況：細胞質(zhì)量低城侧，文庫大小低易遣，它的線粒體占比也低彼妻，但是spike-in占比很高

它可能是因為細胞破壞非常嚴重嫌佑，導(dǎo)致細胞質(zhì)組分基本丟失，最后建庫只捕獲了外源的RNA

而另一種相反的情況：高線粒體占比侨歉，低spike-in占比屋摇，則有可能是說明細胞沒有損傷，只是它的代謝活性很強幽邓。

這種分析思路同樣適用于單核RNA-Seq（snRNA-Seq）炮温，但關(guān)注點從細胞文庫質(zhì)量轉(zhuǎn)移到了細胞核文庫質(zhì)量

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"))

如果在右下角存在許多點時就要注意了，可能細胞受到損傷牵舵；而左上角的點還需要進一步結(jié)合細胞生物學(xué)類型進行判斷

5 針對Droplet數(shù)據(jù)的細胞判斷（例如10X）

5.1 前言

由于這種建庫方法的獨特性柒啤，我們沒辦法事先知道某一個細胞文庫（比如一個cell barcode）是真正包含一個細胞還是只是一個空的液滴（droplet）。因此畸颅，第一步是需要根據(jù)得到的表達量信息担巩，來推斷液滴是不是空的。但僅僅根據(jù)表達量判斷還是不太靠譜（是不是很復(fù)雜的問題没炒？）涛癌，因為還有可能在空的液滴中依然包含外源RNA，最后的細胞文庫依舊不為0送火，但確實不包含細胞拳话。

這里為了說明這個問題，使用了一個未過濾的10X PBMC數(shù)據(jù)

# 學(xué)習一下數(shù)據(jù)下載的方式（新建一個目錄种吸，然后把鏈接放進函數(shù)進行下載）
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"))
# 把數(shù)據(jù)解壓到當前目錄
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)

當然弃衍，如果你想省去下載這一步，還可以用我保存的數(shù)據(jù)【sce.pbmc的網(wǎng)盤鏈接：鏈接：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只展示一個點（去重復(fù)）
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的文庫大小有高有低笨鸡，并且相差較大姜钳，因此可能對應(yīng)著真實存在的細胞和空液滴。當然最簡單的辦法還是”一刀切“形耗，留下那些文庫較大的細胞哥桥，不過還是有損失真實細胞類型的風險。

5.2 【作了解】檢測空的液滴

主要使用emptyDrops()函數(shù)激涤，檢查每個barcode的表達量是不是和外源RNA的表達量有顯著差異(Lun et al. 2019)拟糕。如果存在顯著差異，就說明barcode中更有可能是一個細胞倦踢。這種方法可以幫助區(qū)分：測序質(zhì)量好的空液滴和包含細胞但RNA含量較少的液滴送滞。盡管它們總體的表達量可能很相似，但本質(zhì)不同辱挥，還是要區(qū)分的犁嗅。

這個函數(shù)假設(shè)總體UMI含量低就是空的液滴，使用Monte Carlo統(tǒng)計模擬計算p值晤碘，如果要重復(fù)結(jié)果褂微，需要設(shè)置隨機種子。另外設(shè)置 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

將最后的結(jié)果保存：

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處理之后的）扰楼，而不會是幾十萬的液滴（包含空液滴）呀癣。當然，如果用已過濾的數(shù)據(jù)再去檢測空細胞弦赖，結(jié)果一樣是不準確的项栏。

5.3 與其他質(zhì)控指標的結(jié)合

雖然上一步emptyDrops檢測了空的液滴，但是依然不能知道有細胞的液滴中細胞質(zhì)量如何腾节。因為即使包含細胞忘嫉，還是有可能是死細胞或損傷的細胞。

假設(shè)我們知道這些細胞中沒有代謝活性很高的細胞（也就是說可以根據(jù)線粒體占比去除細胞）：

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")

最后的過濾依然是需要自己對數(shù)據(jù)的把握：

一方面考慮這些篩選出來的”不太合格“的細胞對下游分析的影響到底有多大案腺；
二是考慮如果去掉這些細胞庆冕，會不會損失一些細胞類型？

6 去掉低質(zhì)量細胞

一旦細胞質(zhì)控完成劈榨，可以選擇是直接去掉這些細胞访递，還是先標記出來給個”緩刑“的機會

直接去除很簡單，直接對sce對象的列取子集即可同辣，例如：

filtered <- sce.416b[,!reasons$discard]

不過拷姿，質(zhì)控過程的最大難題是不小心丟掉了某些細胞類型惭载，畢竟質(zhì)控指標總是和生物狀態(tài)有著千絲萬縷的關(guān)系。

好响巢，下面先給那些不符合要求的細胞一個”緩刑“的機會：

計算了舍棄組和保留組細胞的平均表達量以及舍棄組與保留組之間的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)

如果說目前被丟棄的細胞真的是某個特別的細胞類型描滔，那么這個細胞類型應(yīng)該有屬于自己的marker基因，那么這些基因應(yīng)該在被丟棄的這群中顯著高表達踪古『ぃ看到圖中縱軸0以上應(yīng)該是舍棄組比保留組中高表達的基因，但是呢伏穆，這部分基因在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)

發(fā)現(xiàn)這時就會有很多基因在丟棄組中上調(diào)表達，且表達量還不低

根據(jù)先驗知識枕扫，推測有一些血小板相關(guān)基因在丟棄組中高表達（例如：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)

看圖片，果然血小板相關(guān)的基因在丟棄組中高表達

如果感覺自己的過濾太過于嚴格烟瞧，而擔心損失一些細胞诗鸭，那么可以嘗試調(diào)高isOutlier()的nmads=參數(shù)

7 【看看就好】如果不去除，只是標記呢燕刻？

低質(zhì)量細胞還可以只是標記出來只泼，然后繼續(xù)進行下游分析剖笙，這樣為了防止自己隨便過濾而造成細胞類型丟失卵洗。

這樣其實可以在后面的聚類分析中，能看到標記出的這部分細胞可以單獨聚在一起（而我們心里也會清楚弥咪，它們?yōu)楹螘墼谝黄穑?/p>

如何進行標記过蹂？

marked <- sce.416b
marked$discard <- batch.reasons$discard

雖然這樣做很保險，但會增加質(zhì)控結(jié)果的解釋難度聚至。當QC解讀起來太難酷勺，可能又會到后期根據(jù)marker基因去分辨低質(zhì)量和正常的細胞。另外會增加后續(xù)的計算量扳躬。

總結(jié)

作者依然是建議提前去除的脆诉，長痛不如短痛，不要把低質(zhì)量的信息加到后面更加復(fù)雜的計算中贷币。

如果對去除的標準沒有信心击胜，可以在走完一遍常規(guī)流程后，從頭回來再去進行標記役纹，走一下另一條路偶摔，看看和常規(guī)的結(jié)果有什么不同，是不是真的多了一些”容易忽略“的細胞類型促脉？

歡迎關(guān)注我們的公眾號～_～　　
我們是兩個農(nóng)轉(zhuǎn)生信的小碩辰斋，打造生信星球策州，想讓它成為一個不拽術(shù)語、通俗易懂的生信知識平臺宫仗。需要幫助或提出意見請后臺留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!