質(zhì)控

scRNA-seq數(shù)據(jù)中的低質(zhì)量文庫可能有多種來源寇壳，例如解離過程中的細(xì)胞損傷或文庫制備失敺戮！（例如無效的逆轉(zhuǎn)錄或PCR擴(kuò)增）。這些通常表現(xiàn)為總計(jì)數(shù)低怎抛，表達(dá)的基因很少，線粒體或spike-in比例高的“細(xì)胞”芽淡。這些低質(zhì)量的庫是有問題的马绝，因?yàn)樗鼈兛赡軐?dǎo)致下游分析中的誤導(dǎo)性結(jié)果：

1.它們形成了自己獨(dú)特的細(xì)胞群，使對(duì)結(jié)果的解釋變得復(fù)雜挣菲。最明顯的原因是細(xì)胞損傷后線粒體比例增加或核RNA富集富稻。在最壞的情況下，由不同細(xì)胞類型產(chǎn)生的低質(zhì)量文庫可以基于損傷誘導(dǎo)表達(dá)譜中的相似性聚集在一起白胀，從而在其他不同亞群之間形成人為的中間狀態(tài)或軌跡椭赋。此外，由于轉(zhuǎn)換后平均值的變化或杠，非常小的庫可以形成自己的集群哪怔。
2.他們?cè)诜讲罟烙?jì)或主成分分析過程中扭曲了集群異質(zhì)性的特征。前幾個(gè)主要成分將捕獲質(zhì)量差異而不是生物學(xué)差異向抢，從而降低降維效果认境。同樣，差異最大的基因?qū)⒂傻唾|(zhì)量細(xì)胞與高質(zhì)量細(xì)胞之間的差異驅(qū)動(dòng)挟鸠。最明顯的例子：計(jì)數(shù)非常低的低質(zhì)量文庫元暴，其中歸一化放大了那些庫中恰好具有非零計(jì)數(shù)的基因的表觀變異。
3.它們包含的基因似乎由于主動(dòng)縮放以針對(duì)小文庫進(jìn)行標(biāo)準(zhǔn)化而被強(qiáng)烈“上調(diào)”兄猩。這對(duì)于污染的以低但恒定水平存在于所有文庫中的轉(zhuǎn)錄本（例如茉盏，來自環(huán)境溶液）最成問題。在低質(zhì)量文庫中增加的縮放比例會(huì)以較大的標(biāo)準(zhǔn)化表達(dá)值將這些轉(zhuǎn)錄物的計(jì)數(shù)轉(zhuǎn)換為少量枢冤，從而導(dǎo)致與其他細(xì)胞相比明顯的上調(diào)鸠姨。這可能會(huì)產(chǎn)生誤導(dǎo)，因?yàn)槭苡绊懙幕蛲ǔ＞哂猩飳W(xué)敏感性淹真，但實(shí)際上在另一個(gè)亞群中表達(dá)讶迁。
為了避免（或至少減輕）這些問題，我們需要在分析開始時(shí)刪除這些細(xì)胞核蘸。此步驟通常稱為細(xì)胞層面的質(zhì)量控制（QC）巍糯。 （這里我們將互換使用“庫”和“細(xì)胞”，盡管在處理基于液滴的數(shù)據(jù)時(shí)他們的區(qū)別將變得很重要客扎。）我們將演示使用A. T. L. Lun等人的小型scRNA-seq數(shù)據(jù)集祟峦。 該版本沒有事先進(jìn)行質(zhì)量控制，因此我們可以應(yīng)用徙鱼。

#--- loading ---#
library(scRNAseq)
sce.416b <- LunSpikeInData(which="416b") 
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):
## altExpNames(2): ERCC SIRV

質(zhì)控指標(biāo)的選擇

我們使用幾種常見的QC指標(biāo)宅楞，根據(jù)表達(dá)譜來鑒定低質(zhì)量的細(xì)胞针姿。這些指標(biāo)將在下面針對(duì)SMART-seq2數(shù)據(jù)的reads進(jìn)行質(zhì)控，但相同的過程適用于由其他技術(shù)（如MARS-seq和基于液滴的方法）生成的UMI數(shù)據(jù)厌衙。

1.庫大小定義為每個(gè)細(xì)胞的所有相關(guān)feature的總計(jì)數(shù)之和距淫。在這里，我們將相關(guān)feature視為內(nèi)源基因婶希。具有小文庫的細(xì)胞質(zhì)量低下榕暇，因?yàn)樵谖膸熘苽溥^程中某個(gè)時(shí)刻由于細(xì)胞裂解或無效的cDNA捕獲和擴(kuò)增導(dǎo)致RNA丟失。
2.每個(gè)細(xì)胞中表達(dá)feature的數(shù)量定義為該細(xì)胞計(jì)數(shù)為非零的內(nèi)源基因的數(shù)量喻杈。由于尚未成功捕獲多樣化的轉(zhuǎn)錄物種群彤枢，任何具有很少表達(dá)基因的細(xì)胞都可能質(zhì)量較差。
3.計(jì)算“映射到spike-in轉(zhuǎn)錄本的reads數(shù)量”相對(duì)于“每個(gè)細(xì)胞所有feature（包括spike-in）reads的總數(shù)”的比例奕塑。由于應(yīng)該向每個(gè)細(xì)胞中添加相同量的spike-in RNA堂污，因此spike-in計(jì)數(shù)的任何富集都是內(nèi)源RNA丟失的征兆家肯。因此龄砰，高比例表明劣質(zhì)細(xì)胞，可能由于部分細(xì)胞裂解或解離過程中的RNA降解而丟失了內(nèi)源RNA讨衣。

4. 在沒有spike-in轉(zhuǎn)錄物的情況下换棚，可以使用定位到線粒體基因組中基因的reads的比例。高比例表明細(xì)胞質(zhì)量較差反镇，可能是由于穿孔細(xì)胞的細(xì)胞質(zhì)RNA丟失固蚤。理由是，在存在適度損害的情況下歹茶，細(xì)胞膜上的孔允許單個(gè)轉(zhuǎn)錄物分子外排夕玩，但過小而無法使線粒體逸出，從而導(dǎo)致線粒體轉(zhuǎn)錄物相對(duì)富集惊豺。對(duì)于單核RNA-seq實(shí)驗(yàn)燎孟，高比例也很有用，因?yàn)樗鼈兛梢詷?biāo)記尚未成功剝離細(xì)胞質(zhì)的細(xì)胞尸昧。
   對(duì)于每個(gè)細(xì)胞揩页，我們使用來自scater包的perCellQCMetrics（）函數(shù)來計(jì)算這些QC指標(biāo)。 sum列包含每個(gè)細(xì)胞的總數(shù)烹俗，detected列包含檢測(cè)到的基因數(shù)爆侣。 subsets_Mito_percent列包含映射到線粒體轉(zhuǎn)錄本的reads的百分比。 （出于演示目的幢妄，我們展示了兩種確定每個(gè)轉(zhuǎn)錄本的基因組位置的方法兔仰。）最后，altexps_ERCC_percent列包含映射到ERCC轉(zhuǎn)錄本的reads的百分比蕉鸳。

# Retrieving the mitochondrial transcripts using genomic locations included in
# the row-level annotation for the SingleCellExperiment.
location <- rowRanges(sce.416b)
is.mito <- any(seqnames(location)=="MT")

# 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)
df <- perCellQCMetrics(sce.416b, subsets=list(Mito=is.mito))
df

## DataFrame with 192 rows and 12 columns
##                                             sum  detected subsets_Mito_sum
##                                       <numeric> <numeric>        <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1     865936      7618            78790
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1    1076277      7521            98613
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1    1180138      8306           100341
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1    1342593      8143           104882
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1    1668311      7154           129559
## ...                                         ...       ...              ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1    776622      8174            48126
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1   1299950      8956           112225
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1   1800696      9530           135693
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1     46731      6649             3505
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1   1866692     10964           150375
##                                       subsets_Mito_detected
##                                                   <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1                     20
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1                     20
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1                     19
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1                     20
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1                     22
## ...                                                     ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1                    20
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1                    25
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1                    23
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1                    16
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1                    29
##                                       subsets_Mito_percent altexps_ERCC_sum
##                                                  <numeric>        <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1               9.09882            65278
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1               9.16242            74748
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1               8.50248            60878
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1               7.81190            60073
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1               7.76588           136810
## ...                                                    ...              ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1              6.19684            61575
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1              8.63302            94982
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1              7.53559           113707
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1              7.50037             7580
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1              8.05569            48664
##                                       altexps_ERCC_detected
##                                                   <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1                     39
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1                     40
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1                     42
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1                     42
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1                     44
## ...                                                     ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1                    39
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1                    41
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1                    40
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1                    44
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1                    39
##                                       altexps_ERCC_percent altexps_SIRV_sum
##                                                  <numeric>        <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1               6.80658            27828
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1               6.28030            39173
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1               4.78949            30058
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1               4.18567            32542
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1               7.28887            71850
## ...                                                    ...              ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1              7.17620            19848
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1              6.65764            31729
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1              5.81467            41116
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1             13.48898             1883
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1              2.51930            16289
##                                       altexps_SIRV_detected
##                                                   <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1                      7
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1                      7
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1                      7
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1                      7
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1                      7
## ...                                                     ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1                     7
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1                     7
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1                     7
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1                     7
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1                     7
##                                       altexps_SIRV_percent     total
##                                                  <numeric> <numeric>
## SLX-9555.N701_S502.C89V9ANXX.s_1.r_1               2.90165    959042
## SLX-9555.N701_S503.C89V9ANXX.s_1.r_1               3.29130   1190198
## SLX-9555.N701_S504.C89V9ANXX.s_1.r_1               2.36477   1271074
## SLX-9555.N701_S505.C89V9ANXX.s_1.r_1               2.26741   1435208
## SLX-9555.N701_S506.C89V9ANXX.s_1.r_1               3.82798   1876971
## ...                                                    ...       ...
## SLX-11312.N712_S505.H5H5YBBXX.s_8.r_1             2.313165    858045
## SLX-11312.N712_S506.H5H5YBBXX.s_8.r_1             2.224004   1426661
## SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1             2.102562   1955519
## SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1             3.350892     56194
## SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1             0.843271   1931645

或者斋陪，有的人可能更喜歡使用addPerCellQC（）函數(shù)。 這將計(jì)算每個(gè)細(xì)胞的QC統(tǒng)計(jì)信息并將其附加到SingleCellExperiment對(duì)象的colData中，從而使我們能夠?qū)⑺邢嚓P(guān)信息保留在單個(gè)對(duì)象中无虚，以供以后處理缔赠。

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"                 "subsets_Mito_sum"        
## [13] "subsets_Mito_detected"    "subsets_Mito_percent"    
## [15] "altexps_ERCC_sum"         "altexps_ERCC_detected"   
## [17] "altexps_ERCC_percent"     "altexps_SIRV_sum"        
## [19] "altexps_SIRV_detected"    "altexps_SIRV_percent"    
## [21] "total"

此處的關(guān)鍵假設(shè)是QC指標(biāo)與每個(gè)細(xì)胞的生物學(xué)狀態(tài)無關(guān)。 差異（例如友题，低文庫大小嗤堰，高線粒體比例）被認(rèn)為是由技術(shù)因素而不是生物學(xué)過程驅(qū)動(dòng)的，這意味著隨后去除的細(xì)胞不會(huì)在下游分析中歪曲生物學(xué)過程度宦。 嚴(yán)重違反此假設(shè)可能會(huì)導(dǎo)致細(xì)胞類型丟失踢匣，例如該實(shí)驗(yàn)體系本身具有較低的RNA含量或大量的線粒體。 我們可以使用其他診斷工具來檢查此類現(xiàn)象（后續(xù)高級(jí)分析中講解）戈抄。

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

識(shí)別低質(zhì)量單元的最簡(jiǎn)單方法是在QC指標(biāo)上應(yīng)用閾值离唬。 例如，如果庫大小小于100,000reads划鸽，我們可能會(huì)認(rèn)為這些單元的質(zhì)量較低输莺； 表達(dá)少于5,000個(gè)基因； spike-in率超過10％; 或線粒體比例超過10％裸诽。

c.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

雖然簡(jiǎn)單嫂用，但此策略需要大量經(jīng)驗(yàn)才能確定每種實(shí)驗(yàn)方案和生物系統(tǒng)的合適閾值。 基于reads計(jì)數(shù)的數(shù)據(jù)的閾值根本不適用于基于UMI的數(shù)據(jù)丈冬，反之亦然嘱函。 線粒體活性或總RNA含量的差異要求分別針對(duì)不同的生物系統(tǒng)不斷調(diào)整線粒體閾值和spike-in閾值。 確實(shí)埂蕊，即使使用相同的方法和系統(tǒng)往弓，由于cDNA捕獲效率和每細(xì)胞測(cè)序深度的差異，適當(dāng)?shù)拈撝狄部赡芤蜻\(yùn)行而異蓄氧。

自適應(yīng)閾值

識(shí)別異常值

為了獲得合適的閾值函似，我們假設(shè)大多數(shù)數(shù)據(jù)集由高質(zhì)量細(xì)胞組成。然后匀们，我們根據(jù)所有細(xì)胞中每個(gè)指標(biāo)的中位數(shù)的中位數(shù)絕對(duì)偏差（MAD）缴淋，確定對(duì)于各種QC指標(biāo)而言異常的細(xì)胞。具體來說泄朴，如果某個(gè)值在“問題”方向上距離中位數(shù)多于3個(gè)MAD重抖，則認(rèn)為該值是異常值。這種過濾器將保留遵循正態(tài)分布的99％的非異常值祖灰。

對(duì)于416B數(shù)據(jù)钟沛，我們確定對(duì)數(shù)轉(zhuǎn)換的庫大小比中位數(shù)低3個(gè)MAD的細(xì)胞。對(duì)數(shù)轉(zhuǎn)換用type =“ lower”時(shí)以較小的值提高分辨率局扶。具體而言恨统，它保證閾值不是負(fù)值叁扫，這對(duì)于非負(fù)矩陣而言將毫無意義。此外畜埋，文庫大小的分布表現(xiàn)出較重的右尾并不罕見莫绣。對(duì)數(shù)轉(zhuǎn)換可避免MAD膨脹，從而可能損害左尾的異常檢測(cè)悠鞍。 （更一般而言对室，它證明上述99％的理由是合理的。）

qc.lib2 <- isOutlier(df$sum, log=TRUE, type="lower")

我們對(duì)經(jīng)對(duì)數(shù)轉(zhuǎn)化的表達(dá)基因進(jìn)行相同的操作

qc.nexprs2 <- isOutlier(df$detected, log=TRUE, type="lower")

isOutlier() 還會(huì)在輸出向量的屬性中返回每個(gè)指標(biāo)的確切過濾閾值咖祭。 這些對(duì)于檢查自動(dòng)選擇的閾值是否合適非常有用掩宜。

attr(qc.lib2, "thresholds")

##  lower higher 
## 434083    Inf

attr(qc.nexprs2, "thresholds")

##  lower higher 
##   5231    Inf

我們?yōu)榫哂邢嗤δ艿幕诒壤闹笜?biāo)識(shí)別異常值。 這些分布經(jīng)常顯示出較重的右尾么翰，但是與前兩個(gè)指標(biāo)不同牺汤，正是右尾本身包含了假定的低質(zhì)量細(xì)胞。 因此浩嫌，我們不執(zhí)行任何變換來縮小尾巴-而是希望將尾巴中的細(xì)胞識(shí)別為較大的離群值檐迟。 （雖然理論上有可能獲得高于100％的毫無意義的閾值，但這很少見固该，因此不會(huì)引起實(shí)際關(guān)注锅减。）

qc.spike2 <- isOutlier(df$altexps_ERCC_percent, type="higher")
attr(qc.spike2, "thresholds")

##  lower higher 
##   -Inf  14.15

qc.mito2 <- isOutlier(df$subsets_Mito_percent, type="higher")
attr(qc.mito2, "thresholds")

##  lower higher 
##   -Inf  11.92

對(duì)于這些指標(biāo)中的任何一個(gè)異常的細(xì)胞都被認(rèn)為質(zhì)量低下并被丟棄糖儡。

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

或者伐坏，可以使用quickPerCellQC（）函數(shù)將整個(gè)過程一步完成。

reasons <- quickPerCellQC(df, 
    sub.fields=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

使用此策略握联，閾值可適應(yīng)給定指標(biāo)的值分布的位置和分布桦沉。 這使得QC程序可以適應(yīng)測(cè)序深度，cDNA捕獲效率金闽，線粒體含量等方面的變化纯露，而無需任何用戶干預(yù)或事先經(jīng)驗(yàn)。 但是代芜，它確實(shí)需要一些假設(shè)埠褪，下面將對(duì)其進(jìn)行詳細(xì)討論。

離群值檢測(cè)的假設(shè)

離群檢測(cè)假定大多數(shù)細(xì)胞具有可接受的質(zhì)量挤庇。這通常是合理的钞速，并且在某些情況下可以通過肉眼檢查細(xì)胞是否完整（例如在微孔板上）進(jìn)行實(shí)驗(yàn)支持。如果大多數(shù)的細(xì)胞質(zhì)量（不可接受）低嫡秕，則自適應(yīng)閾值顯然會(huì)失敗渴语，因?yàn)樗鼈儫o法刪除大多數(shù)細(xì)胞。當(dāng)然昆咽，在旁觀者眼中驾凶，可接受的與否是隨情境改變的——例如牙甫，眾所周知，神經(jīng)元很難分解调违，而我們通常會(huì)在某個(gè)QC指標(biāo)下將細(xì)胞保留在神經(jīng)元scRNA-seq數(shù)據(jù)集的窟哺，而在更嚴(yán)格的條件下如胚胎干細(xì)胞這個(gè)QC指標(biāo)是不可接受的。

前面提到的另一個(gè)假設(shè)是技肩，質(zhì)控指標(biāo)與每個(gè)細(xì)胞的生物學(xué)狀態(tài)無關(guān)脏答。這在高度異質(zhì)性細(xì)胞群體中最有可能被違反，其中某些細(xì)胞類型天然具有更少的總RNA或更多的線粒體亩鬼。即使沒有捕獲或測(cè)序方面的任何技術(shù)問題殖告，也可能將此類細(xì)胞視為異常值并刪除。通過考慮QC指標(biāo)中的生物變異性雳锋，MAD的使用在某種程度上減輕了該問題黄绩。異類種群在高質(zhì)量細(xì)胞之間的指標(biāo)應(yīng)具有較高的可變性，從而增加了MAD并減少了錯(cuò)誤刪除特定細(xì)胞類型的機(jī)會(huì)（以降低去除低質(zhì)量細(xì)胞的能力為代價(jià)）玷过。

通常爽丹，這些假設(shè)是合理的，或者它們的違反對(duì)下游結(jié)論影響很小辛蚊。盡管如此粤蝎，在解釋結(jié)果時(shí)記住它們還是有幫助的。

考慮實(shí)驗(yàn)因素

更復(fù)雜的研究可能涉及使用不同實(shí)驗(yàn)參數(shù)（例如測(cè)序深度）生成的細(xì)胞批次袋马。在這種情況下初澎，自適應(yīng)策略應(yīng)分別應(yīng)用于每個(gè)批次。從包含多個(gè)批次樣品的混合物分布中計(jì)算中位數(shù)和MAD幾乎沒有意義虑凛。例如碑宴，如果一個(gè)批次中的測(cè)序覆蓋率比其他批次低，則它將拖低中位數(shù)并使MAD膨脹桑谍。這將降低自適應(yīng)閾值對(duì)其他批次的適用性延柠。

如果每個(gè)批次都由其自己的SingleCellExperiment表示，則isOutlier（）函數(shù)可以直接應(yīng)用于每個(gè)批次锣披，如上所示贞间。但是，如果所有批次中的細(xì)胞已合并到單個(gè)SingleCellExperiment中雹仿，則應(yīng)該使用batch =參數(shù)來確保在每個(gè)批次中都識(shí)別出異常值增热。這允許isOutlier（）適應(yīng)批次之間質(zhì)量控制指標(biāo)的系統(tǒng)差異。

我們將再次使用416B數(shù)據(jù)集進(jìn)行說明盅粪，該數(shù)據(jù)集包含兩個(gè)實(shí)驗(yàn)因素-原始和癌基因誘導(dǎo)狀態(tài)钓葫。我們將這些因素結(jié)合在一起，并通過quickPerCellQC（）在isOutlier（）的batch =參數(shù)中使用它票顾。這將導(dǎo)致去除更多的細(xì)胞础浮，因?yàn)椋╥）批次之間測(cè)序深度的系統(tǒng)差異和（ii）致癌基因誘導(dǎo)表達(dá)的基因數(shù)量差異不再使MAD膨脹帆调。

batch <- paste0(sce.416b$phenotype, "-", sce.416b$block)
batch.reasons <- quickPerCellQC(df, batch=batch,
    sub.fields=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

也就是說，batch =的使用包含一個(gè)更強(qiáng)的假設(shè)豆同，即每個(gè)批次中的大多數(shù)單元都是高質(zhì)量的番刊。 如果整個(gè)批次失敗，則異常檢測(cè)將無法充當(dāng)該批次的適當(dāng)QC篩選器影锈。 例如舞竿，Grun等人中有兩個(gè)批次嗜价。人類胰腺數(shù)據(jù)集包含相當(dāng)一部分推定的受損細(xì)胞叫倍，其ERCC含量高于其他批次（圖1）布疙。 這會(huì)使這些批次中的中位數(shù)和MAD膨脹，導(dǎo)致無法刪除假定的低質(zhì)量細(xì)胞辆床。 在這種情況下佳晶，最好從其他批次中計(jì)算出一個(gè)共享的中位數(shù)和MAD，并使用這些估計(jì)值來為有問題的批次中的細(xì)胞獲取適當(dāng)?shù)倪^濾閾值讼载，如下所示轿秧。

library(scRNAseq)
sce.grun <- GrunPancreasData()
sce.grun <- addPerCellQC(sce.grun)

# First attempt with batch-specific thresholds.
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))

# Second attempt, sharing information across batches
# to avoid dramatically different thresholds for unusual batches.
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)

圖1：在Grun胰腺數(shù)據(jù)集的每個(gè)供體中ERCC轉(zhuǎn)錄物比例的分布。 每個(gè)點(diǎn)代表一個(gè)細(xì)胞咨堤，并根據(jù)是在每個(gè)批次中將其識(shí)別為異常值（左）還是使用公共閾值（右）來著色

為了確定有問題的批次菇篡，一個(gè)有用的經(jīng)驗(yàn)法則是找到質(zhì)量控制閾值與其他批次的閾值相比異常的批次。 這里的假設(shè)是一喘，大多數(shù)批次由大多數(shù)高質(zhì)量細(xì)胞組成驱还，因此閾值應(yīng)遵循“典型”批次中的某些單峰分布。 如果我們觀察到一個(gè)批處理具有極高的閾值津滞，我們可能會(huì)懷疑它包含大量會(huì)使每個(gè)批處理的MAD膨脹的低質(zhì)量細(xì)胞铝侵。 我們?cè)谙旅嬉訥run等人數(shù)據(jù)演示此過程灼伤。

ercc.thresholds <- attr(discard.ercc, "thresholds")["higher",]
ercc.thresholds

##     D10     D17      D2      D3      D7 
##  73.611   7.600   6.011 113.106  15.217

names(ercc.thresholds)[isOutlier(ercc.thresholds, type="higher")]

## [1] "D10" "D3"

如果我們不能假設(shè)大多數(shù)批次都包含大多數(shù)高質(zhì)量細(xì)胞触徐，那么所有猜測(cè)都將消失； 我們必須恢復(fù)選擇任意閾值并希望達(dá)到最佳效果的方法狐赡。

其他方法

另一種策略是基于每個(gè)細(xì)胞的QC指標(biāo)來識(shí)別高維空間中的異常值撞鹉。 我們使用來自robustbase的方法基于它們的QC指標(biāo)來量化每個(gè)細(xì)胞的“異常”颖侄，然后使用isOutlier（）來識(shí)別表現(xiàn)出異常高水平異常的低質(zhì)量細(xì)胞鸟雏。

stats <- cbind(log10(df$sum), log10(df$detected),
    df$subsets_Mito_percent, df$altexps_ERCC_percent)

library(robustbase)
outlying <- adjOutlyingness(stats, only.outlyingness = TRUE)
multi.outlier <- isOutlier(outlying, type = "higher")
summary(multi.outlier)

##    Mode   FALSE    TRUE 
## logical     179      13

這種方法和相關(guān)方法（例如基于PCA的異常檢測(cè)和支持向量機(jī)）可以提供更高的能力，以區(qū)分低質(zhì)量的細(xì)胞與高質(zhì)量的細(xì)胞览祖，因?yàn)樗鼈兛梢岳迷S多QC指標(biāo)中的模式孝鹊。 但是，這需要付出一些可解釋性的代價(jià)展蒂，因?yàn)閯h除給定細(xì)胞的原因可能并不總是很明顯又活。

為了完整起見苔咪，我們注意到也可以從基因表達(dá)譜而非質(zhì)量控制指標(biāo)中識(shí)別異常值。 我們認(rèn)為這是一種冒險(xiǎn)的策略柳骄，因?yàn)樗梢匀コ∮屑?xì)胞群中的高質(zhì)量細(xì)胞团赏。

檢查診斷圖

檢查QC指標(biāo)的分布是個(gè)好習(xí)慣（圖2），可以發(fā)現(xiàn)可能的問題耐薯。 在最理想的情況下舔清，我們將看到正態(tài)分布可以證明離群值檢測(cè)中使用的3 MAD閾值是合理的。 其他模型的細(xì)胞表明QC指標(biāo)可能與某種生物學(xué)狀態(tài)相關(guān)曲初，有可能導(dǎo)致過濾過程中不同細(xì)胞類型的損失体谒。 或細(xì)胞亞群的文庫制備存在不一致，這在基于plate的實(shí)驗(yàn)方案中很常見臼婆。 然后营密，可以快速識(shí)別出任何指標(biāo)的系統(tǒng)差值批次，以進(jìn)行進(jìn)一步的故障排除或徹底刪除目锭，就像上面的圖1一樣评汰。

colData(sce.416b) <- cbind(colData(sce.416b), df)
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=1
)

圖2：416B數(shù)據(jù)集中每個(gè)批次和表型的質(zhì)量控制指標(biāo)分布。 每個(gè)點(diǎn)代表一個(gè)細(xì)胞痢虹，并根據(jù)是否被丟棄分別著色被去。

另一個(gè)有用的診斷方法是針對(duì)其他一些QC指標(biāo)繪制線粒體計(jì)數(shù)比例。 目的是要確認(rèn)沒有總計(jì)數(shù)和線粒體計(jì)數(shù)都很高的細(xì)胞奖唯，以確保我們不會(huì)無意中去除代謝活躍的高質(zhì)量細(xì)胞（例如肝細(xì)胞）惨缆。 我們證明了使用來自涉及老鼠大腦的更大實(shí)驗(yàn)的數(shù)據(jù)； 在這種情況下丰捷，我們?cè)趫D3的右上角沒有觀察到任何可能與新陳代謝活躍且未受損的細(xì)胞相對(duì)應(yīng)的點(diǎn)坯墨。

#--- loading ---#
library(scRNAseq)
sce.zeisel <- ZeiselBrainData()

library(scater)
sce.zeisel <- aggregateAcrossFeatures(sce.zeisel, 
    id=sub("_loc[0-9]+$", "", rownames(sce.zeisel)))

#--- gene-annotation ---#
library(org.Mm.eg.db)
rowData(sce.zeisel)$Ensembl <- mapIds(org.Mm.eg.db, 
    keys=rownames(sce.zeisel), keytype="SYMBOL", column="ENSEMBL")

sce.zeisel <- addPerCellQC(sce.zeisel, 
    subsets=list(Mt=rowData(sce.zeisel)$featureType=="mito"))

qc <- quickPerCellQC(colData(sce.zeisel), 
    sub.fields=c("altexps_ERCC_percent", "subsets_Mt_percent"))
sce.zeisel$discard <- qc$discard

plotColData(sce.zeisel, x="sum", y="subsets_Mt_percent", colour_by="discard")

圖3：在Zeisel腦數(shù)據(jù)集中分配給線粒體轉(zhuǎn)錄本的UMI的百分比，相對(duì)于UMI的總數(shù)進(jìn)行繪制（頂部）病往。 每個(gè)點(diǎn)都代表一個(gè)細(xì)胞捣染，并根據(jù)是否被視為劣質(zhì)并被丟棄對(duì)其進(jìn)行著色。

我們看到所有這些指標(biāo)相互之間顯示出較弱的相關(guān)性停巷，大概是細(xì)胞損傷的共同潛在作用的體現(xiàn)耍攘。 相關(guān)性的弱點(diǎn)促使人們使用多種指標(biāo)來捕獲技術(shù)質(zhì)量的不同方面。 當(dāng)然畔勤，不利的一面是蕾各，這些指標(biāo)還可能代表生物學(xué)的不同方面，增加了丟棄整個(gè)細(xì)胞類型的風(fēng)險(xiǎn)庆揪。

移除劣質(zhì)細(xì)胞

識(shí)別出低質(zhì)量的細(xì)胞后式曲，我們可以選擇刪除它們或?qū)ζ溥M(jìn)行標(biāo)記。 刪除是最簡(jiǎn)單的選項(xiàng)缸榛，可以通過按列設(shè)置SingleCellExperiment來實(shí)現(xiàn)吝羞。

＃保留我們不想丟棄的列始鱼。
filtered <- sce.416b[,!reasons$discard]

質(zhì)量控制期間最大的實(shí)際問題是是否會(huì)無意中丟棄整個(gè)細(xì)胞群。 由于QC指標(biāo)永遠(yuǎn)不會(huì)完全獨(dú)立于生物學(xué)狀態(tài)脆贵，因此始終存在發(fā)生這種情況的風(fēng)險(xiǎn)医清。 我們可以通過尋找丟棄細(xì)胞和保留細(xì)胞之間基因表達(dá)的系統(tǒng)差異來判斷細(xì)胞類型是否有丟失。 為了證明這一點(diǎn)卖氨，我們計(jì)算了416B數(shù)據(jù)集中廢棄池和保留池的平均計(jì)數(shù)会烙，并計(jì)算了池平均值之間的對(duì)數(shù)倍變化。

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

如果丟棄的池中富集了某種細(xì)胞類型筒捺，則應(yīng)觀察到相應(yīng)標(biāo)記基因的表達(dá)增加柏腻。 在圖4中，丟棄池中沒有明顯的基因上調(diào)系統(tǒng)系吭，這表明QC步驟并未無意間過濾掉416B數(shù)據(jù)集中的細(xì)胞類型五嫂。

plot(abundance, logFC, xlab="Average count", ylab="Log-FC (lost/kept)", pch=16)
points(abundance[is.mito], logFC[is.mito], col="dodgerblue", pch=16)

圖4：與416B數(shù)據(jù)集中的保留細(xì)胞相比，丟棄的細(xì)胞中表達(dá)的對(duì)數(shù)倍變化肯尺。 每個(gè)點(diǎn)代表一個(gè)帶有藍(lán)色線粒體轉(zhuǎn)錄本的基因沃缘。

為了進(jìn)行比較，讓我們考慮來自10X Genomics的PBMC數(shù)據(jù)集的QC步驟则吟。 我們將對(duì)庫大小應(yīng)用任意固定的閾值來過濾細(xì)胞槐臀，而不是使用任何基于異常值的方法。 具體來說氓仲，我們會(huì)刪除所有庫大小小于500的庫水慨。

#--- loading ---#
library(DropletTestFiles)
raw.path <- getTestFile("tenx-2.1.0-pbmc4k/1.0.0/raw.tar.gz")
out.path <- file.path(tempdir(), "pbmc4k")
untar(raw.path, exdir=out.path)

library(DropletUtils)
fname <- file.path(out.path, "raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)

#--- gene-annotation ---#
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
    rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)

library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID, 
    column="SEQNAME", keytype="GENEID")

#--- cell-detection ---#
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]

discard <- colSums(counts(sce.pbmc)) < 500
lost <- calculateAverage(counts(sce.pbmc)[,discard])
kept <- calculateAverage(counts(sce.pbmc)[,!discard])

logged <- edgeR::cpm(cbind(lost, kept), log=TRUE, prior.count=2)
logFC <- logged[,1] - logged[,2]
abundance <- rowMeans(logged)

在圖5中，廢棄池中存在一個(gè)獨(dú)特的種群敬扛，這是一組在丟失中強(qiáng)烈上調(diào)的基因晰洒。 這包括PF4，PPBP和SDPR啥箭，它們（（擾流板警報(bào)５骸）表示存在被alt.discard丟棄的血小板群。

plot(abundance, logFC, xlab="Average count", ylab="Log-FC (lost/kept)", pch=16)
platelet <- c("PF4", "PPBP", "SDPR")
points(abundance[platelet], logFC[platelet], col="orange", pch=16)

圖5

6.6標(biāo)記劣質(zhì)細(xì)胞

另一種選擇是照這樣標(biāo)記低質(zhì)量的細(xì)胞，并將其保留在下游分析中粪牲。 此處的目的是允許形成低質(zhì)量細(xì)胞的群古瓤，然后在解釋結(jié)果時(shí)識(shí)別并忽略這些群。 這種方法避免了丟棄質(zhì)量控制指標(biāo)值較差的細(xì)胞類型腺阳，從而為用戶提供了機(jī)會(huì)來決定此類細(xì)胞的群是否代表真正的生物學(xué)狀態(tài)落君。

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

缺點(diǎn)是它將QC的負(fù)擔(dān)轉(zhuǎn)移到了細(xì)胞群的解釋上，這已經(jīng)是scRNA-seq數(shù)據(jù)分析的瓶頸亭引。確實(shí)叽奥，如果我們不相信QC指標(biāo)，我們將不得不僅基于標(biāo)記基因來區(qū)分真正的細(xì)胞類型和低質(zhì)量的細(xì)胞痛侍，由于后者傾向于“表達(dá)”有趣的基因朝氓，因此這并不總是那么容易。保留低質(zhì)量細(xì)胞也會(huì)損害方差建模的準(zhǔn)確性主届，例如赵哲，需要使用更多的PC來抵消早期PC受低質(zhì)量細(xì)胞和其他細(xì)胞之間的差異驅(qū)動(dòng)的事實(shí)。

對(duì)于常規(guī)分析君丁，我們建議默認(rèn)執(zhí)行清除操作枫夺，以避免低質(zhì)量細(xì)胞引起的并發(fā)問題。這樣一來绘闷，大多數(shù)群結(jié)構(gòu)的特征都無需擔(dān)心或至少不會(huì)擔(dān)心其有效性橡庞。完成初始分析后，如果對(duì)丟棄的細(xì)胞類型有任何疑問印蔗，可以在僅標(biāo)記低質(zhì)量細(xì)胞的情況下進(jìn)行更徹底的重新分析扒最。這樣可以恢復(fù)具有低RNA含量，高線粒體比例等的細(xì)胞類型华嘹，僅在初始分析中“填補(bǔ)空白”的情況下才需要對(duì)其進(jìn)行解釋吧趣。