質(zhì)控
scRNA-seq數(shù)據(jù)中的低質(zhì)量文庫(kù)可能有多種來(lái)源，例如解離過(guò)程中的細(xì)胞損傷或文庫(kù)制備失?。ɡ鐭o(wú)效的逆轉(zhuǎn)錄或PCR擴(kuò)增）。這些通常表現(xiàn)為總計(jì)數(shù)低，表達(dá)的基因很少，線(xiàn)粒體或spike-in比例高的“細(xì)胞”。這些低質(zhì)量的庫(kù)是有問(wèn)題的，因?yàn)樗鼈兛赡軐?dǎo)致下游分析中的誤導(dǎo)性結(jié)果：
1.它們形成了自己獨(dú)特的細(xì)胞群，使對(duì)結(jié)果的解釋變得復(fù)雜。最明顯的原因是細(xì)胞損傷后線(xiàn)粒體比例增加或核RNA富集。在最壞的情況下，由不同細(xì)胞類(lèi)型產(chǎn)生的低質(zhì)量文庫(kù)可以基于損傷誘導(dǎo)表達(dá)譜中的相似性聚集在一起，從而在其他不同亞群之間形成人為的中間狀態(tài)或軌跡。此外，由于轉(zhuǎn)換后平均值的變化，非常小的庫(kù)可以形成自己的集群。
2.他們?cè)诜讲罟烙?jì)或主成分分析過(guò)程中扭曲了集群異質(zhì)性的特征。前幾個(gè)主要成分將捕獲質(zhì)量差異而不是生物學(xué)差異，從而降低降維效果。同樣，差異最大的基因?qū)⒂傻唾|(zhì)量細(xì)胞與高質(zhì)量細(xì)胞之間的差異驅(qū)動(dòng)。最明顯的例子：計(jì)數(shù)非常低的低質(zhì)量文庫(kù)，其中歸一化放大了那些庫(kù)中恰好具有非零計(jì)數(shù)的基因的表觀(guān)變異。
3.它們包含的基因似乎由于主動(dòng)縮放以針對(duì)小文庫(kù)進(jìn)行標(biāo)準(zhǔn)化而被強(qiáng)烈“上調(diào)”。這對(duì)于污染的以低但恒定水平存在于所有文庫(kù)中的轉(zhuǎn)錄本（例如，來(lái)自環(huán)境溶液）最成問(wèn)題。在低質(zhì)量文庫(kù)中增加的縮放比例會(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á)。
為了避免（或至少減輕）這些問(wèn)題，我們需要在分析開(kāi)始時(shí)刪除這些細(xì)胞。此步驟通常稱(chēng)為細(xì)胞層面的質(zhì)量控制（QC）。 （這里我們將互換使用“庫(kù)”和“細(xì)胞”，盡管在處理基于液滴的數(shù)據(jù)時(shí)他們的區(qū)別將變得很重要。）我們將演示使用A. T. L. Lun等人的小型scRNA-seq數(shù)據(jù)集。 該版本沒(méi)有事先進(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)的選擇
我們使用幾種常見(jiàn)的QC指標(biāo)，根據(jù)表達(dá)譜來(lái)鑒定低質(zhì)量的細(xì)胞。這些指標(biāo)將在下面針對(duì)SMART-seq2數(shù)據(jù)的reads進(jìn)行質(zhì)控，但相同的過(guò)程適用于由其他技術(shù)（如MARS-seq和基于液滴的方法）生成的UMI數(shù)據(jù)。
1.庫(kù)大小定義為每個(gè)細(xì)胞的所有相關(guān)feature的總計(jì)數(shù)之和。在這里，我們將相關(guān)feature視為內(nèi)源基因。具有小文庫(kù)的細(xì)胞質(zhì)量低下，因?yàn)樵谖膸?kù)制備過(guò)程中某個(gè)時(shí)刻由于細(xì)胞裂解或無(wú)效的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ì)胞裂解或解離過(guò)程中的RNA降解而丟失了內(nèi)源RNA。
- 在沒(méi)有spike-in轉(zhuǎn)錄物的情況下，可以使用定位到線(xiàn)粒體基因組中基因的reads的比例。高比例表明細(xì)胞質(zhì)量較差，可能是由于穿孔細(xì)胞的細(xì)胞質(zhì)RNA丟失。理由是，在存在適度損害的情況下，細(xì)胞膜上的孔允許單個(gè)轉(zhuǎn)錄物分子外排，但過(guò)小而無(wú)法使線(xiàn)粒體逸出，從而導(dǎo)致線(xiàn)粒體轉(zhuǎn)錄物相對(duì)富集。對(duì)于單核RNA-seq實(shí)驗(yàn)，高比例也很有用，因?yàn)樗鼈兛梢詷?biāo)記尚未成功剝離細(xì)胞質(zhì)的細(xì)胞。
對(duì)于每個(gè)細(xì)胞，我們使用來(lái)自scater包的perCellQCMetrics（）函數(shù)來(lái)計(jì)算這些QC指標(biāo)。sum列包含每個(gè)細(xì)胞的總數(shù)，detected列包含檢測(cè)到的基因數(shù)。subsets_Mito_percent列包含映射到線(xiàn)粒體轉(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)無(wú)關(guān)。 差異（例如，低文庫(kù)大小，高線(xiàn)粒體比例）被認(rèn)為是由技術(shù)因素而不是生物學(xué)過(guò)程驅(qū)動(dòng)的，這意味著隨后去除的細(xì)胞不會(huì)在下游分析中歪曲生物學(xué)過(guò)程。 嚴(yán)重違反此假設(shè)可能會(huì)導(dǎo)致細(xì)胞類(lèi)型丟失，例如該實(shí)驗(yàn)體系本身具有較低的RNA含量或大量的線(xiàn)粒體。 我們可以使用其他診斷工具來(lái)檢查此類(lèi)現(xiàn)象（后續(xù)高級(jí)分析中講解）。
鑒定低質(zhì)量細(xì)胞
識(shí)別低質(zhì)量單元的最簡(jiǎn)單方法是在QC指標(biāo)上應(yīng)用閾值。 例如，如果庫(kù)大小小于100,000reads，我們可能會(huì)認(rèn)為這些單元的質(zhì)量較低； 表達(dá)少于5,000個(gè)基因； spike-in率超過(guò)10％; 或線(xiàn)粒體比例超過(guò)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ù)，反之亦然。 線(xiàn)粒體活性或總RNA含量的差異要求分別針對(duì)不同的生物系統(tǒng)不斷調(diào)整線(xiàn)粒體閾值和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ì)胞。具體來(lái)說(shuō)，如果某個(gè)值在“問(wèn)題”方向上距離中位數(shù)多于3個(gè)MAD，則認(rèn)為該值是異常值。這種過(guò)濾器將保留遵循正態(tài)分布的99％的非異常值。
對(duì)于416B數(shù)據(jù)，我們確定對(duì)數(shù)轉(zhuǎn)換的庫(kù)大小比中位數(shù)低3個(gè)MAD的細(xì)胞。對(duì)數(shù)轉(zhuǎn)換用type =“ lower”時(shí)以較小的值提高分辨率。具體而言，它保證閾值不是負(fù)值，這對(duì)于非負(fù)矩陣而言將毫無(wú)意義。此外，文庫(kù)大小的分布表現(xiàn)出較重的右尾并不罕見(jià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)的確切過(guò)濾閾值。 這些對(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í)行任何變換來(lái)縮小尾巴-而是希望將尾巴中的細(xì)胞識(shí)別為較大的離群值。 （雖然理論上有可能獲得高于100％的毫無(wú)意義的閾值，但這很少見(jiàn)，因此不會(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è)過(guò)程一步完成。
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捕獲效率，線(xiàn)粒體含量等方面的變化，而無(wú)需任何用戶(hù)干預(yù)或事先經(jīng)驗(yàn)。 但是，它確實(shí)需要一些假設(shè)，下面將對(duì)其進(jìn)行詳細(xì)討論。
離群值檢測(cè)的假設(shè)
離群檢測(cè)假定大多數(shù)細(xì)胞具有可接受的質(zhì)量。這通常是合理的，并且在某些情況下可以通過(guò)肉眼檢查細(xì)胞是否完整（例如在微孔板上）進(jìn)行實(shí)驗(yàn)支持。如果大多數(shù)的細(xì)胞質(zhì)量（不可接受）低，則自適應(yīng)閾值顯然會(huì)失敗，因?yàn)樗鼈儫o(wú)法刪除大多數(shù)細(xì)胞。當(dāng)然，在旁觀(guān)者眼中，可接受的與否是隨情境改變的——例如，眾所周知，神經(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)無(wú)關(guān)。這在高度異質(zhì)性細(xì)胞群體中最有可能被違反，其中某些細(xì)胞類(lèi)型天然具有更少的總RNA或更多的線(xiàn)粒體。即使沒(méi)有捕獲或測(cè)序方面的任何技術(shù)問(wèn)題，也可能將此類(lèi)細(xì)胞視為異常值并刪除。通過(guò)考慮QC指標(biāo)中的生物變異性，MAD的使用在某種程度上減輕了該問(wèn)題。異類(lèi)種群在高質(zhì)量細(xì)胞之間的指標(biāo)應(yīng)具有較高的可變性，從而增加了MAD并減少了錯(cuò)誤刪除特定細(xì)胞類(lèi)型的機(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幾乎沒(méi)有意義。例如，如果一個(gè)批次中的測(cè)序覆蓋率比其他批次低，則它將拖低中位數(shù)并使MAD膨脹。這將降低自適應(yīng)閾值對(duì)其他批次的適用性。
如果每個(gè)批次都由其自己的SingleCellExperiment表示，則isOutlier（）函數(shù)可以直接應(yīng)用于每個(gè)批次，如上所示。但是，如果所有批次中的細(xì)胞已合并到單個(gè)SingleCellExperiment中，則應(yīng)該使用batch =參數(shù)來(lái)確保在每個(gè)批次中都識(shí)別出異常值。這允許isOutlier（）適應(yīng)批次之間質(zhì)量控制指標(biāo)的系統(tǒng)差異。
我們將再次使用416B數(shù)據(jù)集進(jìn)行說(shuō)明，該數(shù)據(jù)集包含兩個(gè)實(shí)驗(yàn)因素-原始和癌基因誘導(dǎo)狀態(tài)。我們將這些因素結(jié)合在一起，并通過(guò)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
也就是說(shuō)，batch =的使用包含一個(gè)更強(qiáng)的假設(shè)，即每個(gè)批次中的大多數(shù)單元都是高質(zhì)量的。 如果整個(gè)批次失敗，則異常檢測(cè)將無(wú)法充當(dāng)該批次的適當(dāng)QC篩選器。 例如，Grun等人中有兩個(gè)批次。人類(lèi)胰腺數(shù)據(jù)集包含相當(dāng)一部分推定的受損細(xì)胞，其ERCC含量高于其他批次（圖1）。 這會(huì)使這些批次中的中位數(shù)和MAD膨脹，導(dǎo)致無(wú)法刪除假定的低質(zhì)量細(xì)胞。 在這種情況下，最好從其他批次中計(jì)算出一個(gè)共享的中位數(shù)和MAD，并使用這些估計(jì)值來(lái)為有問(wèn)題的批次中的細(xì)胞獲取適當(dāng)?shù)倪^(guò)濾閾值，如下所示。
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)
為了確定有問(wèn)題的批次，一個(gè)有用的經(jīng)驗(yàn)法則是找到質(zhì)量控制閾值與其他批次的閾值相比異常的批次。 這里的假設(shè)是，大多數(shù)批次由大多數(shù)高質(zhì)量細(xì)胞組成，因此閾值應(yīng)遵循“典型”批次中的某些單峰分布。 如果我們觀(guān)察到一個(gè)批處理具有極高的閾值，我們可能會(huì)懷疑它包含大量會(huì)使每個(gè)批處理的MAD膨脹的低質(zhì)量細(xì)胞。 我們?cè)谙旅嬉訥run等人數(shù)據(jù)演示此過(guò)程。
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)來(lái)識(shí)別高維空間中的異常值。 我們使用來(lái)自robustbase的方法基于它們的QC指標(biāo)來(lái)量化每個(gè)細(xì)胞的“異?！?，然后使用isOutlier（）來(lái)識(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ì)胞的原因可能并不總是很明顯。
為了完整起見(jiàn)，我們注意到也可以從基因表達(dá)譜而非質(zhì)量控制指標(biāo)中識(shí)別異常值。 我們認(rèn)為這是一種冒險(xiǎn)的策略，因?yàn)樗梢匀コ∮屑?xì)胞群中的高質(zhì)量細(xì)胞。
檢查診斷圖
檢查QC指標(biāo)的分布是個(gè)好習(xí)慣（圖2），可以發(fā)現(xiàn)可能的問(wèn)題。 在最理想的情況下，我們將看到正態(tài)分布可以證明離群值檢測(cè)中使用的3 MAD閾值是合理的。 其他模型的細(xì)胞表明QC指標(biāo)可能與某種生物學(xué)狀態(tài)相關(guān)，有可能導(dǎo)致過(guò)濾過(guò)程中不同細(xì)胞類(lèi)型的損失。 或細(xì)胞亞群的文庫(kù)制備存在不一致，這在基于plate的實(shí)驗(yàn)方案中很常見(jià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
)
另一個(gè)有用的診斷方法是針對(duì)其他一些QC指標(biāo)繪制線(xiàn)粒體計(jì)數(shù)比例。 目的是要確認(rèn)沒(méi)有總計(jì)數(shù)和線(xiàn)粒體計(jì)數(shù)都很高的細(xì)胞，以確保我們不會(huì)無(wú)意中去除代謝活躍的高質(zhì)量細(xì)胞（例如肝細(xì)胞）。 我們證明了使用來(lái)自涉及老鼠大腦的更大實(shí)驗(yàn)的數(shù)據(jù)； 在這種情況下，我們?cè)趫D3的右上角沒(méi)有觀(guān)察到任何可能與新陳代謝活躍且未受損的細(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")
我們看到所有這些指標(biāo)相互之間顯示出較弱的相關(guān)性，大概是細(xì)胞損傷的共同潛在作用的體現(xiàn)。 相關(guān)性的弱點(diǎn)促使人們使用多種指標(biāo)來(lái)捕獲技術(shù)質(zhì)量的不同方面。 當(dāng)然，不利的一面是，這些指標(biāo)還可能代表生物學(xué)的不同方面，增加了丟棄整個(gè)細(xì)胞類(lèi)型的風(fēng)險(xiǎn)。
移除劣質(zhì)細(xì)胞
識(shí)別出低質(zhì)量的細(xì)胞后，我們可以選擇刪除它們或?qū)ζ溥M(jìn)行標(biāo)記。 刪除是最簡(jiǎn)單的選項(xiàng)，可以通過(guò)按列設(shè)置SingleCellExperiment來(lái)實(shí)現(xiàn)。
＃保留我們不想丟棄的列。
filtered <- sce.416b[,!reasons$discard]
質(zhì)量控制期間最大的實(shí)際問(wèn)題是是否會(huì)無(wú)意中丟棄整個(gè)細(xì)胞群。 由于QC指標(biāo)永遠(yuǎn)不會(huì)完全獨(dú)立于生物學(xué)狀態(tài)，因此始終存在發(fā)生這種情況的風(fēng)險(xiǎn)。 我們可以通過(guò)尋找丟棄細(xì)胞和保留細(xì)胞之間基因表達(dá)的系統(tǒng)差異來(lái)判斷細(xì)胞類(lèi)型是否有丟失。 為了證明這一點(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ì)胞類(lèi)型，則應(yīng)觀(guān)察到相應(yīng)標(biāo)記基因的表達(dá)增加。 在圖4中，丟棄池中沒(méi)有明顯的基因上調(diào)系統(tǒng)，這表明QC步驟并未無(wú)意間過(guò)濾掉416B數(shù)據(jù)集中的細(xì)胞類(lèi)型。
plot(abundance, logFC, xlab="Average count", ylab="Log-FC (lost/kept)", pch=16)
points(abundance[is.mito], logFC[is.mito], col="dodgerblue", pch=16)
為了進(jìn)行比較，讓我們考慮來(lái)自10X Genomics的PBMC數(shù)據(jù)集的QC步驟。 我們將對(duì)庫(kù)大小應(yīng)用任意固定的閾值來(lái)過(guò)濾細(xì)胞，而不是使用任何基于異常值的方法。 具體來(lái)說(shuō)，我們會(huì)刪除所有庫(kù)大小小于500的庫(kù)。
#--- 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)?。┍硎敬嬖诒籥lt.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)
如果我們懷疑細(xì)胞類(lèi)型已被我們的質(zhì)量控制程序錯(cuò)誤地丟棄，則最直接的解決方案是放寬質(zhì)量控制過(guò)濾器，以獲取與真正的生物學(xué)差異相關(guān)的指標(biāo)。例如，可以通過(guò)在
isOutlier（）調(diào)用中增加nmads =來(lái)放寬異常檢測(cè)。當(dāng)然，這增加了保留更多低質(zhì)量細(xì)胞的風(fēng)險(xiǎn)。順便說(shuō)一句，值得一提的是，細(xì)胞的真正技術(shù)質(zhì)量也可能與其類(lèi)型相關(guān)。 （這與細(xì)胞類(lèi)型和QC指標(biāo)之間的相關(guān)性不同，因?yàn)楹笳呤俏覀儗?duì)質(zhì)量的不完美代理。）如果某些細(xì)胞類(lèi)型在scRNA-seq協(xié)議期間不適合解離或微流體處理，則可能出現(xiàn)這種情況。在這種情況下，如果QC的所有細(xì)胞都損壞，則可以“正確”丟棄整個(gè)細(xì)胞類(lèi)型。確實(shí)，與實(shí)驗(yàn)方案中的損失相比，對(duì)QC期間細(xì)胞類(lèi)型的計(jì)算去除的關(guān)注可能很小。
6.6標(biāo)記劣質(zhì)細(xì)胞
另一種選擇是照這樣標(biāo)記低質(zhì)量的細(xì)胞，并將其保留在下游分析中。 此處的目的是允許形成低質(zhì)量細(xì)胞的群，然后在解釋結(jié)果時(shí)識(shí)別并忽略這些群。 這種方法避免了丟棄質(zhì)量控制指標(biāo)值較差的細(xì)胞類(lèi)型，從而為用戶(hù)提供了機(jī)會(huì)來(lái)決定此類(lèi)細(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)記基因來(lái)區(qū)分真正的細(xì)胞類(lèi)型和低質(zhì)量的細(xì)胞，由于后者傾向于“表達(dá)”有趣的基因，因此這并不總是那么容易。保留低質(zhì)量細(xì)胞也會(huì)損害方差建模的準(zhǔn)確性，例如，需要使用更多的PC來(lái)抵消早期PC受低質(zhì)量細(xì)胞和其他細(xì)胞之間的差異驅(qū)動(dòng)的事實(shí)。
對(duì)于常規(guī)分析，我們建議默認(rèn)執(zhí)行清除操作，以避免低質(zhì)量細(xì)胞引起的并發(fā)問(wèn)題。這樣一來(lái)，大多數(shù)群結(jié)構(gòu)的特征都無(wú)需擔(dān)心或至少不會(huì)擔(dān)心其有效性。完成初始分析后，如果對(duì)丟棄的細(xì)胞類(lèi)型有任何疑問(wèn)，可以在僅標(biāo)記低質(zhì)量細(xì)胞的情況下進(jìn)行更徹底的重新分析。這樣可以恢復(fù)具有低RNA含量，高線(xiàn)粒體比例等的細(xì)胞類(lèi)型，僅在初始分析中“填補(bǔ)空白”的情況下才需要對(duì)其進(jìn)行解釋。
