胞轉(zhuǎn)錄組質(zhì)控全解析：用Seurat精準(zhǔn)過濾低質(zhì)量細(xì)胞)
1. 我為什么要把質(zhì)控這件事單獨寫一篇如果你跟我一樣接手過幾套單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)那你大概率體會過這種感覺數(shù)據(jù)下載完暴力跑完Seurat全流程滿心歡喜打開UMAP圖發(fā)現(xiàn)細(xì)胞團全部糊成一團marker基因的熱圖也畫不出層次感。這個時候你第一反應(yīng)是這個樣本生物學(xué)就是這么復(fù)雜其實大概率是剛開始的質(zhì)控步驟沒做對。三年前我第一次處理10X平臺的單細(xì)胞數(shù)據(jù)時就是這副慘樣。當(dāng)時項目催得急我拿著默認(rèn)參數(shù)就直接開跑過濾后的對象里有12000多個細(xì)胞看起來數(shù)量挺漂亮但整個數(shù)據(jù)怎么調(diào)都出不了清晰的分群。后來靜下心來重新看數(shù)據(jù)發(fā)現(xiàn)光質(zhì)控這一關(guān)就濾掉了接近三分之一——細(xì)胞數(shù)從12000掉到8000出頭但重新跑出來的UMAP分群干凈利落marker基因的表達模式終于跟文獻對上了。那次返工讓我徹底明白了一個道理單細(xì)胞數(shù)據(jù)分析的成敗一半決定在質(zhì)控。不是說后面的聚類、差異分析不重要而是如果進來的細(xì)胞本身就有大量的破損細(xì)胞、雙重體、空滴后面所有結(jié)果都會受到污染。這篇文章我不打算講整套單細(xì)胞分析流程而是把用Seurat過濾低質(zhì)量細(xì)胞這一件事拆開揉碎從指標(biāo)含義、閾值設(shè)定到完整代碼、常見翻車場景一次性講透。這篇文章適合剛開始接觸單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析的研究生、想優(yōu)化自己質(zhì)控流程的生信從業(yè)者也適合那些跑完流程但結(jié)果總是不對勁、懷疑自己哪里做錯了的同學(xué)。2. 數(shù)據(jù)到手先別急先理解Seurat對象里的三個質(zhì)控指標(biāo)很多新手拿到10X的filtered_feature_bc_matrix文件夾第一件事就是趕緊讀進Seurat開始跑這其實是在跳過最關(guān)鍵的一步你得先知道你手里的數(shù)據(jù)是什么樣子的、哪些細(xì)胞其實根本不能用。在Seurat里做質(zhì)控翻來覆去就三個核心指標(biāo)nFeature_RNA、nCount_RNA、percent.mt。下面一個一個說清楚。2.1 nFeature_RNA每個細(xì)胞檢測到的基因數(shù)nFeature_RNA直譯過來就是檢測到的基因數(shù)量代表這個細(xì)胞里到底有多少種基因被捕捉到了。它是判斷細(xì)胞狀態(tài)最敏感的指標(biāo)之一。一個健康完整的細(xì)胞在10X平臺上通常能檢測到幾百到幾千個基因。如果檢測到的基因數(shù)特別低比如只有二三十個那這個所謂的細(xì)胞大概率是一個被抽破了、mRNA幾乎漏光了的空泡或者是一個即將凋亡、轉(zhuǎn)錄本已經(jīng)降解得差不多的垂死細(xì)胞。反過來如果nFeature_RNA高得離譜比如超過一萬那它很可能是兩個甚至多個細(xì)胞被同一個GEM微滴包裹了也就是所謂雙重體doublet或多細(xì)胞聚合體。這里有一個很容易被忽略的細(xì)節(jié)nFeature_RNA的判斷標(biāo)準(zhǔn)跟細(xì)胞類型是強相關(guān)的。中性粒細(xì)胞、某些靜止的T細(xì)胞亞群它們本身轉(zhuǎn)錄本就很低nFeature只有幾百是正常的而腫瘤細(xì)胞或高度活化的B細(xì)胞nFeature高也說得通。所以你不能一概而論地認(rèn)為nFeature小于500就是垃圾細(xì)胞也不能認(rèn)為nFeature大于6000就一定是雙重體。質(zhì)控的精髓恰恰在于結(jié)合你的樣本類型和數(shù)據(jù)分布去判斷。2.2 nCount_RNA每個細(xì)胞的UMI總數(shù)nCount_RNA代表這個細(xì)胞里所有基因的UMI計數(shù)總和可以粗略理解成這個細(xì)胞被捕獲到的轉(zhuǎn)錄本總量。它和nFeature_RNA是高度相關(guān)的——轉(zhuǎn)錄本總量越高自然檢出的基因種類也會更多。你在做質(zhì)控的時候通常會看到nCount_RNA低得離譜的那一批往往跟nFeature低的細(xì)胞高度重合也就是同一條低質(zhì)量線。但兩邊的信息不完全冗余有些細(xì)胞盡管nFeature看著還行但nCount極低說明這個細(xì)胞整體的轉(zhuǎn)錄活躍度已經(jīng)很低了很可能是處于應(yīng)激或者瀕死狀態(tài)。Seurat官方的質(zhì)控教程里也特別強調(diào)nFeature和nCount要組合起來看這就是為什么要畫FeatureScatter的原因——兩個指標(biāo)做散點圖能很直觀地看出哪些點偏離了正常的相關(guān)性區(qū)間。2.3 線粒體基因比例判斷細(xì)胞狀態(tài)的重要線索percent.mt指的是線粒體基因的UMI計數(shù)占總UMI的比例。這個指標(biāo)算是質(zhì)控里的黃金指標(biāo)因為它在生物學(xué)機制上特別說得通一個健康細(xì)胞擁有完整的細(xì)胞膜和細(xì)胞質(zhì)細(xì)胞質(zhì)里的mRNA整體是保留完好的而一個細(xì)胞膜破損、胞質(zhì)丟失的垂死細(xì)胞細(xì)胞質(zhì)里的mRNA大量流失但線粒體因為存在于相對穩(wěn)定的細(xì)胞器結(jié)構(gòu)中往往留在細(xì)胞碎片里于是線粒體轉(zhuǎn)錄本占總轉(zhuǎn)錄本的比例就會異常升高。用人話翻譯一下線粒體比例高說明這個細(xì)胞的漿漏了大概率處于凋亡或機械破損狀態(tài)。這個指標(biāo)在過濾時極其有效因為它是獨立于nFeature和nCount的另一個維度能把前面兩個指標(biāo)看不出來的漏漿細(xì)胞揪出來。這也是為什么Seurat官方教程里會算出一個percent.mt來幾乎沒有一個單細(xì)胞項目會跳過它。Seurat里有一行固定操作pbmc[[percent.mt]] - PercentageFeatureSet(pbmc, pattern ^MT-)。注意這行的人源數(shù)據(jù)是^MT-小鼠則是^mt-大小寫不要搞錯否則算出來的線粒體比例永遠(yuǎn)是0。這個低級錯誤我在多個同學(xué)的腳本里見過血的教訓(xùn)。3. 完整代碼從Cell Ranger輸出到一個干凈的單細(xì)胞對象理解了指標(biāo)含義之后就該上手操作了。下面這套代碼是我在實際項目中打磨過的版本可以直接復(fù)制修改使用但我不建議無腦跑因為每個步驟背后的參數(shù)含義你需要看懂。3.1 準(zhǔn)備環(huán)境與加載數(shù)據(jù)質(zhì)控這一步對R和Seurat的版本要求并不苛刻R 4.2以上、Seurat 4.x或5.x都行。我目前主要用Seurat 5但下面這套代碼在4.x上同樣能跑。需要的基礎(chǔ)依賴包括dplyr和Matrix通常裝Seurat的時候會自動帶上來。library(Seurat) library(dplyr) library(Matrix) # 讀取Cell Ranger的filtered_feature_bc_matrix目錄 # 目錄下應(yīng)該有barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz三個文件 counts - Read10X(data.dir filtered_feature_bc_matrix/) # 如果拿到的是h5文件這樣讀 # counts - Read10X_h5(filtered_feature_bc_matrix.h5)這里有個小知識點Read10X()讀進來的是一個稀疏矩陣行是基因列是細(xì)胞barcode。Cell Ranger在得到這個矩陣之前其實做了一步空滴過濾EmptyDrops把很大一部分沒有細(xì)胞但背景RNA污染產(chǎn)生的barcode去掉了。但即便經(jīng)過了這一步矩陣?yán)锸Ｏ碌腷arcode也仍然混雜著大量低質(zhì)量細(xì)胞——Cell Ranger的空滴過濾只能篩掉絕對空篩不掉破了但還有殘留轉(zhuǎn)錄本的細(xì)胞。這就是為什么我們還需要Seurat來再做一輪質(zhì)控。3.2 創(chuàng)建Seurat對象并計算質(zhì)控指標(biāo)# 創(chuàng)建Seurat對象 sce - CreateSeuratObject(counts counts, project my_project, min.cells 3, min.features 200)min.cells 3的意思是一個基因至少在3個細(xì)胞里檢測到才被保留。這是為了過濾掉那些極低表達、基本屬于背景噪聲的基因一般設(shè)在3~10之間太小吃不干凈噪聲太大會誤傷一些真實但低表達的基因。min.features 200的意思是一個細(xì)胞至少檢測到200個基因才被保留這一步做的是一次非常粗的預(yù)過濾先把那些一看就是空泡的垃圾barcode剔除掉。注意這個200只是一個保底的初始值真正的精細(xì)過濾在下一步。然后是計算線粒體比例# 計算線粒體基因比例人源用 ^MT-小鼠用 ^mt- sce[[percent.mt]] - PercentageFeatureSet(sce, pattern ^MT-)PercentageFeatureSet()的本質(zhì)就是把匹配了^MT-這些基因的總UMI計數(shù)求和再除以這個細(xì)胞所有基因的總UMI最后乘以100。所以這個指標(biāo)是一個百分比數(shù)值范圍從0到100。3.3 質(zhì)控指標(biāo)可視化拿到三個指標(biāo)之后先別急著設(shè)閾值一定要可視化。不看圖就設(shè)閾值約等于閉著眼開車。下面這幾行代碼就能幫你快速了解數(shù)據(jù)長什么樣# 3個質(zhì)控指標(biāo)的分布小提琴圖 VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3, pt.size 0.01) # nCount與percent.mt的散點圖 plot1 - FeatureScatter(sce, feature1 nCount_RNA, feature2 percent.mt) # nCount與nFeature的散點圖 plot2 - FeatureScatter(sce, feature1 nCount_RNA, feature2 nFeature_RNA) plot1 plot2小提琴圖能讓你一眼看到三個指標(biāo)的分布形態(tài)和異常尾巴。散點圖則能幫你檢查這些指標(biāo)的聯(lián)動關(guān)系是否合理。正常情況下nCount和nFeature應(yīng)該呈明顯的正相關(guān)如果你看到一大團細(xì)胞nCount很低但nFeature卻很高或者反過來那都是信號異常需要格外注意。3.4 核心過濾操作看完了分布圖才輪到核心的過濾代碼。以下是我在多數(shù)外周血和腫瘤組織樣本上驗證過的一套初始參數(shù)# 過濾 # 閾值需要結(jié)合你的樣本類型和數(shù)據(jù)分布靈活調(diào)整下面只是常見起始值 sce_filtered - subset(sce, subset nFeature_RNA 200 nFeature_RNA 5000 percent.mt 15)注意subset()函數(shù)會把不符合條件的細(xì)胞直接從對象中移除。過濾完之后建議再看一次小提琴圖和散點圖確認(rèn)過濾后的數(shù)據(jù)分布是否合理細(xì)胞數(shù)變化多少。這一步是所有質(zhì)控中最核心的一環(huán)也是最容易被忽視的一環(huán)——很多人過濾完就直接進入NormalizeData了根本不做驗證。3.5 過濾后的輸出與常用檢查# 查看過濾前后細(xì)胞數(shù)變化 ncol(sce) ncol(sce_filtered) # 查看過濾后對象的詳細(xì)信息 sce_filtered # 查看每個細(xì)胞的質(zhì)控指標(biāo)統(tǒng)計 summary(sce_filtered$nFeature_RNA) summary(sce_filtered$nCount_RNA) summary(sce_filtered$percent.mt)如果過濾比例適中一般濾掉10%~30%說明數(shù)據(jù)質(zhì)量還可以。如果一次性濾掉了一半以上就要回頭想想是組織解離過程出了問題還是閾值設(shè)得太狠還是這個樣本本身就比較臟這些都需要結(jié)合實際情況判斷而不是簡單調(diào)低閾值把細(xì)胞硬留下來。4. 閾值不是拍腦袋基于分布特征確定質(zhì)控參數(shù)上面的代碼里給的nFeature_RNA 200 nFeature_RNA 5000 percent.mt 15是一個起步參考值離精準(zhǔn)還有距離。為什么這么說因為這個參數(shù)組合是我在特定樣本類型上調(diào)出來的換一個樣本類型、換一個測序平臺它就可能完全不適用。Seurat官方教程里的經(jīng)典參數(shù)是線粒體比例小于5%那是基于PBMC數(shù)據(jù)集給出的示例很多初學(xué)者直接把它套到自己的腫瘤或組織樣本上結(jié)果濾掉了大量真實細(xì)胞。4.1 默認(rèn)參數(shù)為什么經(jīng)常翻車單細(xì)胞轉(zhuǎn)錄組的本質(zhì)是抽樣測序抽樣的效率跟建庫批次、10X試劑盒版本v2、v3、Chromium X、測序深度都有關(guān)系。同樣的細(xì)胞類型在不同批次里測出來的nFeature中位數(shù)可能相差上千。如果你拿一個中位nFeature只有800的冷凍腫瘤樣本去套用PBMC的200~2500區(qū)間那你的前半段正常細(xì)胞就會被無情誤殺。更讓人頭疼的是細(xì)胞類型差異。我處理過一個中性粒細(xì)胞比例很高的樣本中性粒細(xì)胞nFeature中位數(shù)不到500按nFeature必須大于600的常規(guī)思路過濾整個中性粒細(xì)胞群幾乎被團滅?？墒侨思抑行粤＜?xì)胞本來就是短命、轉(zhuǎn)錄本少的細(xì)胞你把它們?nèi)珓h了后續(xù)分析還怎么做這就是精準(zhǔn)兩個字的分量過濾不是要套公式而是要理解你的數(shù)據(jù)里有哪些細(xì)胞、它們的指標(biāo)形態(tài)是怎樣的。4.2 用MAD法則和分布形態(tài)來定閾值那怎么定閾值才算靠譜我自己的習(xí)慣是兩步走先看分布再用分位數(shù)輔助判斷。第一步用核密度圖看nFeature和nCount的整體形態(tài)特別注意有沒有一個矮子峰# 畫核密度圖觀察nFeature_RNA的分布 plot(density(sce$nFeature_RNA), main nFeature density, col steelblue) # 觀察percent.mt的分布 plot(density(sce$percent.mt), main percent.mt density, col firebrick)如果nFeature_RNA的密度圖在左邊低基因數(shù)區(qū)域出現(xiàn)一個獨立的小峰或者峰型很扁那說明數(shù)據(jù)里確實混了一批低質(zhì)量的細(xì)胞。而percent.mt的密度圖如果出現(xiàn)一個拖得很長的尾巴或者主峰明顯右偏那線粒體高比例的細(xì)胞一定不在少數(shù)。第二步就是更定量的做法。很多人會采用中位數(shù)±3倍絕對中位差MAD來幫助確定閾值。它的邏輯跟平均值±3倍標(biāo)準(zhǔn)差類似但MAD對異常值更不敏感更適合單細(xì)胞這種厚尾分布的數(shù)據(jù)# 用中位數(shù)±3倍MAD估計合理范圍僅供初步參考 lower_bound - median(sce$nFeature_RNA) - 3 * mad(sce$nFeature_RNA) upper_bound - median(sce$nFeature_RNA) 3 * mad(sce$nFeature_RNA) cat(nFeature_RNA 參考范圍:, lower_bound, -, upper_bound, \n)這套方法的好處是閾值完全由你自己的數(shù)據(jù)驅(qū)動而不是生搬硬套別人的經(jīng)驗。但它同樣有局限性——如果你的數(shù)據(jù)里低質(zhì)量細(xì)胞占比太高MAD本身就會被它們拉著走導(dǎo)致算出來的范圍過寬。所以我的建議是MAD算出來的范圍用來做參考最終閾值還是回到密度圖上人工確認(rèn)一下看看切的位置是否正好落在主峰和低質(zhì)量小峰之間的洼地上。4.3 不同樣本類型的初始閾值參考為了讓你有個大概的手感我把自己處理過的幾類樣本的初始閾值整理成了表格。注意這只是第一輪粗過濾的參考值不代表所有項目通用樣本類型nFeature_RNA 參考下限nFeature_RNA 參考上限percent.mt 參考上限PBMC 外周血200~4004000~500010~15腫瘤組織新鮮解離300~5006000~800015~20冷凍組織200~3004000~600020~25中性粒細(xì)胞富集樣本100~2003000~400010~15snRNA-seq核測序500~8008000~100001~5表格里的數(shù)值跨度比較大因為每批數(shù)據(jù)的質(zhì)量波動實在太大。我一直強調(diào)一個觀點表格里的數(shù)字只是用來入門的不是你交論文時的最終答案。最終答案永遠(yuǎn)要靠你對自己樣本的理解和數(shù)據(jù)分布的實際觀察來回答。這也是精準(zhǔn)質(zhì)控和應(yīng)付差事式質(zhì)控的核心差別。5. 過濾完別急著開心怎么驗證質(zhì)控效果過濾代碼跑通了細(xì)胞數(shù)也減少了并不代表質(zhì)控做完了。很多人在這時候馬上進入NormalizeData、PCA、UMAP的下一個階段結(jié)果最后發(fā)現(xiàn)分群有問題才想起回頭找質(zhì)控的問題。我更建議你在過濾之后用幾分鐘做幾項驗證把異常在早期就暴露出來。5.1 過濾前后指標(biāo)的對比最直接的驗證就是對比過濾前后的指標(biāo)分布# 過濾前后指標(biāo)分布對比 VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3) VlnPlot(sce_filtered, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3)過濾后你應(yīng)該看到nFeature_RNA和nCount_RNA的低值尾被切掉了整體分布更集中percent.mt的高值尾巴被削平大部分細(xì)胞集中在一個相對正常的區(qū)間小提琴圖整體變得干凈不再有明顯的離群點。如果過濾后你的nFeature主峰反而被砍掉一大半或者percent.mt的分布幾乎沒有變化那說明閾值設(shè)得有問題要么誤傷了大群正常細(xì)胞要么壓根沒把低質(zhì)量細(xì)胞濾除干凈。5.2 檢查高表達基因列表是否合理這是一個非常容易被忽略但很聰明的檢查方法。正常情況下一個健康單細(xì)胞數(shù)據(jù)集中高表達基因應(yīng)當(dāng)是各類細(xì)胞類型的功能基因比如PBMC里應(yīng)該是MALAT1、B2M、各類核糖體蛋白基因等。但如果你的數(shù)據(jù)里混入了大量垂死細(xì)胞或破裂細(xì)胞高表達基因排行榜會被線粒體基因霸榜——因為線粒體基因在破損細(xì)胞里占比太高了。我有個習(xí)慣性操作質(zhì)控后跑一下每個細(xì)胞的top基因看看排名靠前的是不是有大量MT-基因。具體可以這樣# 檢查每個細(xì)胞中最顯著高表達的基因 # 取前500個高表達基因 top_genes - rownames(sce_filtered)[ order(Matrix::rowSums(sce_filteredassays$RNA$counts), decreasing TRUE) ] head(top_genes, 30)如果你看到head出來的30個基因里有大量MT-開頭的基因那就要警惕線粒體基因過濾可能還不夠干凈尤其在有核紅細(xì)胞、心肌組織等線粒體天然豐富的樣本里。至于要不要進一步往上提高percent.mt閾值同樣還是回到分布圖和生物學(xué)實際情況來決定。5.3 雙重體到底要不要處理過濾驗證階段最容易被忽略的其實是雙重體doublets。Seurat的質(zhì)控指標(biāo)可以篩掉低質(zhì)量細(xì)胞但對雙重體的效果非常有限——因為一個雙重體的nFeature通常并不是特別離譜它可能是兩個中等基因數(shù)的細(xì)胞粘在一起生成的nFeature值恰好落在正常范圍的高值區(qū)蒙混過關(guān)。要嚴(yán)格去雙就得用專門的工具。目前最常用的是DoubletFinder# 這個步驟建議在質(zhì)控完成、數(shù)據(jù)標(biāo)準(zhǔn)化之后運行 # 簡單用法示意 # library(DoubletFinder) # sce_filtered - NormalizeData(sce_filtered) # sce_filtered - FindVariableFeatures(sce_filtered) # sce_filtered - ScaleData(sce_filtered) # sce_filtered - RunPCA(sce_filtered) # sce_filtered - FindNeighbors(sce_filtered, dims 1:20) # sce_filtered - FindClusters(sce_filtered, resolution 0.5) # sce_filtered - RunUMAP(sce_filtered, dims 1:20) # # 然后基于聚類結(jié)果預(yù)測doublet # sweep.res - paramSweep(sce_filtered, PCs 1:20, sct FALSE) # sweep.stats - summarizeSweep(sweep.res, GT FALSE) # bcmvn - find.pK(sweep.stats) # ...DoubletFinder的原理是通過人工模擬doublet看哪些細(xì)胞在特征空間上與模擬的doublet最相似從而打分預(yù)測。這個工具能篩掉相當(dāng)一部分漏網(wǎng)的doublet但它也有誤殺風(fēng)險而且是基于聚類結(jié)果預(yù)測的如果前面聚類本身就不好預(yù)測的可靠性也會受影響。所以我通常建議在分群結(jié)果出來后、確定marker基因之前結(jié)合doublet注釋結(jié)果做一次排查盡量不要在只有幾十個細(xì)胞的支持下強行去雙。5.4 與后續(xù)分析流程的正確銜接質(zhì)控是起點不是終點。過濾干凈之后后續(xù)步驟必須銜接好# 完整的下游分析銜接 sce_filtered - NormalizeData(sce_filtered, normalization.method LogNormalize, scale.factor 10000) sce_filtered - FindVariableFeatures(sce_filtered, selection.method vst, nfeatures 2000) sce_filtered - ScaleData(sce_filtered) sce_filtered - RunPCA(sce_filtered, npcs 30) sce_filtered - RunUMAP(sce_filtered, dims 1:20) sce_filtered - FindNeighbors(sce_filtered, dims 1:20) sce_filtered - FindClusters(sce_filtered, resolution 0.5) # 后續(xù)差異分析 markers - FindAllMarkers(sce_filtered, only.pos TRUE, min.pct 0.25, logfc.threshold 0.25)這里我不展開講每個步驟的細(xì)節(jié)但要注意的是質(zhì)控越靠前后面全部步驟都能受益如果你在整合多個樣本最好先分別對每個樣本做各自的質(zhì)控再進入整合流程。把不同質(zhì)量的樣本混在一起再統(tǒng)一過濾往往會出現(xiàn)好的樣本被拉到差的樣本水平這種尷尬處境。6. 實戰(zhàn)中的特殊樣本那些容易翻車的質(zhì)控場景寫代碼部分的質(zhì)控大家都能跑通真正拉開差距的是處理那些非典型樣本。我在多個項目里踩過坑下面幾個場景幾乎每個做單細(xì)胞的人都可能遇到專門拿出來說。6.1 免疫細(xì)胞樣本的線粒體比例天然偏高如果你是做免疫細(xì)胞、尤其是組織浸潤免疫細(xì)胞的一定不要迷信Seurat官方教程的percent.mt 5。5%這個閾值來自PBMC外周血的樣本背景干凈、細(xì)胞狀態(tài)比較均一。但組織里的T細(xì)胞、巨噬細(xì)胞在解離過程中經(jīng)歷了機械剪切、酶消化應(yīng)激狀態(tài)下線粒體活性會升高percent.mt就很容易沖到10%~20%甚至更高。我處理過一批腫瘤浸潤T細(xì)胞樣本percent.mt的中位數(shù)就有14%如果按5%切直接砍掉一大半T細(xì)胞后面根本沒法分析。這種情況下我的建議是把閾值放到20%甚至25%同時結(jié)合nFeature和nCount來判斷——一個T細(xì)胞如果percent.mt是20%但nFeature有1200它大概率還是一條好漢但如果percent.mt 20%同時nFeature只有150那它大概率確實快死了。6.2 腫瘤樣本低質(zhì)量與惡性信號混雜的難題腫瘤組織是單細(xì)胞項目里最復(fù)雜的一類樣本。腫瘤細(xì)胞本身存在拷貝數(shù)變異、倍性異常、增殖活躍它們的nFeature和nCount往往很高高到如果用nFeature_RNA 5000一刀切會把真正的腫瘤細(xì)胞當(dāng)成doublet去掉。另一方面腫瘤微環(huán)境里又有大量浸潤的免疫細(xì)胞、基質(zhì)細(xì)胞各類細(xì)胞之間的nFeature跨度極大。在這種樣本里做質(zhì)控我的策略是在解離后先做一個偏保守的粗過濾確保把明顯的垃圾先清掉然后在聚類之后用inferCNV等工具確認(rèn)腫瘤細(xì)胞群單獨看它們的質(zhì)控指標(biāo)分布再做一輪針對性的二次檢視。這一步很關(guān)鍵——如果你在前期就下手太狠可能把一個重要的腫瘤亞群給做沒了后面再怎么注釋都找不回來。6.3 snRNA-seq核測序千萬不要當(dāng)成普通單細(xì)胞處理單細(xì)胞核測序snRNA-seq這幾年特別火尤其在冷凍組織中優(yōu)勢明顯因為核不容易受到解離損傷的影響。但它的質(zhì)控指標(biāo)跟普通單細(xì)胞完全不一樣核里面不含完整的細(xì)胞質(zhì)線粒體基因本來就極少percent.mt通常只有1%以下同時核轉(zhuǎn)錄本的捕獲效率更高nFeature中位數(shù)往往比完整細(xì)胞更高。如果有人拿著普通單細(xì)胞的標(biāo)準(zhǔn)去過濾snRNA-seq數(shù)據(jù)大概率會把數(shù)據(jù)搞壞。snRNA-seq的質(zhì)控重點是排除帶大量細(xì)胞質(zhì)殘留的完整細(xì)胞碎片和空的細(xì)胞核而不是糾結(jié)線粒體比例。另外snRNA-seq中低nFeature并不一定代表低質(zhì)量——有些靜止細(xì)胞的核轉(zhuǎn)錄本本來就少。如果后面還要結(jié)合ATAC數(shù)據(jù)做多組學(xué)分析質(zhì)控更要謹(jǐn)慎。6.4 污染與背景RNA的影響做過組織樣本的朋友應(yīng)該都有體會組織解離之后總有一些細(xì)胞碎片和游離RNA它們會被隨機包裹進GEM微滴。這種背景RNA污染會造成一個現(xiàn)象所有細(xì)胞的nFeature和nCount都被稀釋了——具體來說垃圾細(xì)胞帶著高背景RNA會讓整體數(shù)據(jù)出現(xiàn)一個低質(zhì)量團塊。Seurat的普通質(zhì)控對這種污染有一定的清理效果但未必能完全解決。如果背景污染比較嚴(yán)重我更推薦用CellBender或者SoupX這類專門的去污染工具在質(zhì)控之前先把背景RNA信號扣除掉。這個操作能讓后續(xù)閾值判斷更精準(zhǔn)但要注意這類工具對計算資源和參數(shù)設(shè)置有一定要求不是所有項目都必需。一般建議是如果聚類后所有細(xì)胞類型都帶有同樣一組高表達基因、UMAP圖分群模糊并且高表達基因列表里有明顯的環(huán)境RNA特征比如各種血紅蛋白基因才需要考慮去污染。7. 寫在最后的實操心得文章寫到這兒代碼和理論都講得差不多了。最后分享幾點只有自己上手跑過幾批數(shù)據(jù)才會總結(jié)出來的體會。第一質(zhì)控閾值不要追求一次到位。我自己現(xiàn)在跑一個新數(shù)據(jù)集很少會第一次就定出完美的閾值。通常是先按照分布圖設(shè)一個大概范圍過濾后看看細(xì)胞數(shù)和UMAP效果再回來微調(diào)。這個跑一遍—看結(jié)果—回來調(diào)的循環(huán)在單細(xì)胞項目里極其正常不要覺得麻煩。第二質(zhì)控日志一定要記。具體包括過濾前多少細(xì)胞、過濾后多少細(xì)胞、三個指標(biāo)的閾值各是多少、當(dāng)時看了哪些圖、為什么這么定。不要高估自己的記憶力。一篇論文從投稿到返修可能隔了大半年如果到時候?qū)徃迦藛柲隳愕倪^濾標(biāo)準(zhǔn)是怎么確定的而你連跑時的參數(shù)都忘了那場面會很尷尬。我現(xiàn)在每個項目都會留一個quality_control_notes.md文件隨手記錄當(dāng)時的判斷邏輯。第三不同樣本之間不要套用同一個閾值。哪怕兩個樣本來自同一個病人的同一類組織只要批次不同、建庫時間不同就應(yīng)該各看各的分布再分別定閾值。把兩個樣本合并成一個對象后統(tǒng)一過濾是新手最容易犯的錯誤之一它會直接導(dǎo)致高質(zhì)量的樣本被你誤傷。第四質(zhì)控的上限不只是不濾錯還要考慮下一步需要什么。比如你后面想做擬時序分析那G0/G1期的細(xì)胞就要盡量保留如果后面想做轉(zhuǎn)錄因子活性分析轉(zhuǎn)錄本種類不夠豐富的細(xì)胞可能支撐不起結(jié)論。每一種下游分析對數(shù)據(jù)質(zhì)量的要求都不同。質(zhì)控參數(shù)不是越嚴(yán)格越好而是越契合你的分析目標(biāo)越好。單細(xì)胞轉(zhuǎn)錄組質(zhì)控這條路說簡單也簡單就是看三個指標(biāo)、設(shè)三個閾值說復(fù)雜也復(fù)雜因為每個數(shù)據(jù)集都有它自己的脾氣。如果你把這篇內(nèi)容看懂了再回到自己的數(shù)據(jù)上去實踐幾輪應(yīng)該能明顯感受到精準(zhǔn)質(zhì)控和套參數(shù)跑流程之間的差距。希望這些經(jīng)驗?zāi)軒湍闵僮咭稽c彎路。 單細(xì)胞轉(zhuǎn)錄組質(zhì)控全解析：用Seurat精準(zhǔn)過濾低質(zhì)量細(xì)胞
1. 我為什么要把質(zhì)控這件事單獨寫一篇如果你跟我一樣接手過幾套單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)那你大概率體會過這種感覺數(shù)據(jù)下載完暴力跑完Seurat全流程滿心歡喜打開UMAP圖發(fā)現(xiàn)細(xì)胞團全部糊成一團marker基因的熱圖也畫不出層次感。這個時候你第一反應(yīng)是這個樣本生物學(xué)就是這么復(fù)雜其實大概率是剛開始的質(zhì)控步驟沒做對。三年前我第一次處理10X平臺的單細(xì)胞數(shù)據(jù)時就是這副慘樣。當(dāng)時項目催得急我拿著默認(rèn)參數(shù)就直接開跑過濾后的對象里有12000多個細(xì)胞看起來數(shù)量挺漂亮但整個數(shù)據(jù)怎么調(diào)都出不了清晰的分群。后來靜下心來重新看數(shù)據(jù)發(fā)現(xiàn)光質(zhì)控這一關(guān)就濾掉了接近三分之一——細(xì)胞數(shù)從12000掉到8000出頭但重新跑出來的UMAP分群干凈利落marker基因的表達模式終于跟文獻對上了。那次返工讓我徹底明白了一個道理單細(xì)胞數(shù)據(jù)分析的成敗一半決定在質(zhì)控。不是說后面的聚類、差異分析不重要而是如果進來的細(xì)胞本身就有大量的破損細(xì)胞、雙重體、空滴后面所有結(jié)果都會受到污染。這篇文章我不打算講整套單細(xì)胞分析流程而是把用Seurat過濾低質(zhì)量細(xì)胞這一件事拆開揉碎從指標(biāo)含義、閾值設(shè)定到完整代碼、常見翻車場景一次性講透。這篇文章適合剛開始接觸單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析的研究生、想優(yōu)化自己質(zhì)控流程的生信從業(yè)者也適合那些跑完流程但結(jié)果總是不對勁、懷疑自己哪里做錯了的同學(xué)。2. 數(shù)據(jù)到手先別急先理解Seurat對象里的三個質(zhì)控指標(biāo)很多新手拿到10X的filtered_feature_bc_matrix文件夾第一件事就是趕緊讀進Seurat開始跑這其實是在跳過最關(guān)鍵的一步你得先知道你手里的數(shù)據(jù)是什么樣子的、哪些細(xì)胞其實根本不能用。在Seurat里做質(zhì)控翻來覆去就三個核心指標(biāo)nFeature_RNA、nCount_RNA、percent.mt。下面一個一個說清楚。2.1 nFeature_RNA每個細(xì)胞檢測到的基因數(shù)nFeature_RNA直譯過來就是檢測到的基因數(shù)量代表這個細(xì)胞里到底有多少種基因被捕捉到了。它是判斷細(xì)胞狀態(tài)最敏感的指標(biāo)之一。一個健康完整的細(xì)胞在10X平臺上通常能檢測到幾百到幾千個基因。如果檢測到的基因數(shù)特別低比如只有二三十個那這個所謂的細(xì)胞大概率是一個被抽破了、mRNA幾乎漏光了的空泡或者是一個即將凋亡、轉(zhuǎn)錄本已經(jīng)降解得差不多的垂死細(xì)胞。反過來如果nFeature_RNA高得離譜比如超過一萬那它很可能是兩個甚至多個細(xì)胞被同一個GEM微滴包裹了也就是所謂雙重體doublet或多細(xì)胞聚合體。這里有一個很容易被忽略的細(xì)節(jié)nFeature_RNA的判斷標(biāo)準(zhǔn)跟細(xì)胞類型是強相關(guān)的。中性粒細(xì)胞、某些靜止的T細(xì)胞亞群它們本身轉(zhuǎn)錄本就很低nFeature只有幾百是正常的而腫瘤細(xì)胞或高度活化的B細(xì)胞nFeature高也說得通。所以你不能一概而論地認(rèn)為nFeature小于500就是垃圾細(xì)胞也不能認(rèn)為nFeature大于6000就一定是雙重體。質(zhì)控的精髓恰恰在于結(jié)合你的樣本類型和數(shù)據(jù)分布去判斷。2.2 nCount_RNA每個細(xì)胞的UMI總數(shù)nCount_RNA代表這個細(xì)胞里所有基因的UMI計數(shù)總和可以粗略理解成這個細(xì)胞被捕獲到的轉(zhuǎn)錄本總量。它和nFeature_RNA是高度相關(guān)的——轉(zhuǎn)錄本總量越高自然檢出的基因種類也會更多。你在做質(zhì)控的時候通常會看到nCount_RNA低得離譜的那一批往往跟nFeature低的細(xì)胞高度重合也就是同一條低質(zhì)量線。但兩邊的信息不完全冗余有些細(xì)胞盡管nFeature看著還行但nCount極低說明這個細(xì)胞整體的轉(zhuǎn)錄活躍度已經(jīng)很低了很可能是處于應(yīng)激或者瀕死狀態(tài)。Seurat官方的質(zhì)控教程里也特別強調(diào)nFeature和nCount要組合起來看這就是為什么要畫FeatureScatter的原因——兩個指標(biāo)做散點圖能很直觀地看出哪些點偏離了正常的相關(guān)性區(qū)間。2.3 線粒體基因比例判斷細(xì)胞狀態(tài)的重要線索percent.mt指的是線粒體基因的UMI計數(shù)占總UMI的比例。這個指標(biāo)算是質(zhì)控里的黃金指標(biāo)因為它在生物學(xué)機制上特別說得通一個健康細(xì)胞擁有完整的細(xì)胞膜和細(xì)胞質(zhì)細(xì)胞質(zhì)里的mRNA整體是保留完好的而一個細(xì)胞膜破損、胞質(zhì)丟失的垂死細(xì)胞細(xì)胞質(zhì)里的mRNA大量流失但線粒體因為存在于相對穩(wěn)定的細(xì)胞器結(jié)構(gòu)中往往留在細(xì)胞碎片里于是線粒體轉(zhuǎn)錄本占總轉(zhuǎn)錄本的比例就會異常升高。用人話翻譯一下線粒體比例高說明這個細(xì)胞的漿漏了大概率處于凋亡或機械破損狀態(tài)。這個指標(biāo)在過濾時極其有效因為它是獨立于nFeature和nCount的另一個維度能把前面兩個指標(biāo)看不出來的漏漿細(xì)胞揪出來。這也是為什么Seurat官方教程里會算出一個percent.mt來幾乎沒有一個單細(xì)胞項目會跳過它。Seurat里有一行固定操作pbmc[[percent.mt]] - PercentageFeatureSet(pbmc, pattern ^MT-)。注意這行的人源數(shù)據(jù)是^MT-小鼠則是^mt-大小寫不要搞錯否則算出來的線粒體比例永遠(yuǎn)是0。這個低級錯誤我在多個同學(xué)的腳本里見過血的教訓(xùn)。3. 完整代碼從Cell Ranger輸出到一個干凈的單細(xì)胞對象理解了指標(biāo)含義之后就該上手操作了。下面這套代碼是我在實際項目中打磨過的版本可以直接復(fù)制修改使用但我不建議無腦跑因為每個步驟背后的參數(shù)含義你需要看懂。3.1 準(zhǔn)備環(huán)境與加載數(shù)據(jù)質(zhì)控這一步對R和Seurat的版本要求并不苛刻R 4.2以上、Seurat 4.x或5.x都行。我目前主要用Seurat 5但下面這套代碼在4.x上同樣能跑。需要的基礎(chǔ)依賴包括dplyr和Matrix通常裝Seurat的時候會自動帶上來。library(Seurat) library(dplyr) library(Matrix) # 讀取Cell Ranger的filtered_feature_bc_matrix目錄 # 目錄下應(yīng)該有barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz三個文件 counts - Read10X(data.dir filtered_feature_bc_matrix/) # 如果拿到的是h5文件這樣讀 # counts - Read10X_h5(filtered_feature_bc_matrix.h5)這里有個小知識點Read10X()讀進來的是一個稀疏矩陣行是基因列是細(xì)胞barcode。Cell Ranger在得到這個矩陣之前其實做了一步空滴過濾EmptyDrops把很大一部分沒有細(xì)胞但背景RNA污染產(chǎn)生的barcode去掉了。但即便經(jīng)過了這一步矩陣?yán)锸Ｏ碌腷arcode也仍然混雜著大量低質(zhì)量細(xì)胞——Cell Ranger的空滴過濾只能篩掉絕對空篩不掉破了但還有殘留轉(zhuǎn)錄本的細(xì)胞。這就是為什么我們還需要Seurat來再做一輪質(zhì)控。3.2 創(chuàng)建Seurat對象并計算質(zhì)控指標(biāo)# 創(chuàng)建Seurat對象 sce - CreateSeuratObject(counts counts, project my_project, min.cells 3, min.features 200)min.cells 3的意思是一個基因至少在3個細(xì)胞里檢測到才被保留。這是為了過濾掉那些極低表達、基本屬于背景噪聲的基因一般設(shè)在3~10之間太小吃不干凈噪聲太大會誤傷一些真實但低表達的基因。min.features 200的意思是一個細(xì)胞至少檢測到200個基因才被保留這一步做的是一次非常粗的預(yù)過濾先把那些一看就是空泡的垃圾barcode剔除掉。注意這個200只是一個保底的初始值真正的精細(xì)過濾在下一步。然后是計算線粒體比例# 計算線粒體基因比例人源用 ^MT-小鼠用 ^mt- sce[[percent.mt]] - PercentageFeatureSet(sce, pattern ^MT-)PercentageFeatureSet()的本質(zhì)就是把匹配了^MT-這些基因的總UMI計數(shù)求和再除以這個細(xì)胞所有基因的總UMI最后乘以100。所以這個指標(biāo)是一個百分比數(shù)值范圍從0到100。3.3 質(zhì)控指標(biāo)可視化拿到三個指標(biāo)之后先別急著設(shè)閾值一定要可視化。不看圖就設(shè)閾值約等于閉著眼開車。下面這幾行代碼就能幫你快速了解數(shù)據(jù)長什么樣# 3個質(zhì)控指標(biāo)的分布小提琴圖 VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3, pt.size 0.01) # nCount與percent.mt的散點圖 plot1 - FeatureScatter(sce, feature1 nCount_RNA, feature2 percent.mt) # nCount與nFeature的散點圖 plot2 - FeatureScatter(sce, feature1 nCount_RNA, feature2 nFeature_RNA) plot1 plot2小提琴圖能讓你一眼看到三個指標(biāo)的分布形態(tài)和異常尾巴。散點圖則能幫你檢查這些指標(biāo)的聯(lián)動關(guān)系是否合理。正常情況下nCount和nFeature應(yīng)該呈明顯的正相關(guān)如果你看到一大團細(xì)胞nCount很低但nFeature卻很高或者反過來那都是信號異常需要格外注意。3.4 核心過濾操作看完了分布圖才輪到核心的過濾代碼。以下是我在多數(shù)外周血和腫瘤組織樣本上驗證過的一套初始參數(shù)# 過濾 # 閾值需要結(jié)合你的樣本類型和數(shù)據(jù)分布靈活調(diào)整下面只是常見起始值 sce_filtered - subset(sce, subset nFeature_RNA 200 nFeature_RNA 5000 percent.mt 15)注意subset()函數(shù)會把不符合條件的細(xì)胞直接從對象中移除。過濾完之后建議再看一次小提琴圖和散點圖確認(rèn)過濾后的數(shù)據(jù)分布是否合理細(xì)胞數(shù)變化多少。這一步是所有質(zhì)控中最核心的一環(huán)也是最容易被忽視的一環(huán)——很多人過濾完就直接進入NormalizeData了根本不做驗證。3.5 過濾后的輸出與常用檢查# 查看過濾前后細(xì)胞數(shù)變化 ncol(sce) ncol(sce_filtered) # 查看過濾后對象的詳細(xì)信息 sce_filtered # 查看每個細(xì)胞的質(zhì)控指標(biāo)統(tǒng)計 summary(sce_filtered$nFeature_RNA) summary(sce_filtered$nCount_RNA) summary(sce_filtered$percent.mt)如果過濾比例適中一般濾掉10%~30%說明數(shù)據(jù)質(zhì)量還可以。如果一次性濾掉了一半以上就要回頭想想是組織解離過程出了問題還是閾值設(shè)得太狠還是這個樣本本身就比較臟這些都需要結(jié)合實際情況判斷而不是簡單調(diào)低閾值把細(xì)胞硬留下來。4. 閾值不是拍腦袋基于分布特征確定質(zhì)控參數(shù)上面的代碼里給的nFeature_RNA 200 nFeature_RNA 5000 percent.mt 15是一個起步參考值離精準(zhǔn)還有距離。為什么這么說因為這個參數(shù)組合是我在特定樣本類型上調(diào)出來的換一個樣本類型、換一個測序平臺它就可能完全不適用。Seurat官方教程里的經(jīng)典參數(shù)是線粒體比例小于5%那是基于PBMC數(shù)據(jù)集給出的示例很多初學(xué)者直接把它套到自己的腫瘤或組織樣本上結(jié)果濾掉了大量真實細(xì)胞。4.1 默認(rèn)參數(shù)為什么經(jīng)常翻車單細(xì)胞轉(zhuǎn)錄組的本質(zhì)是抽樣測序抽樣的效率跟建庫批次、10X試劑盒版本v2、v3、Chromium X、測序深度都有關(guān)系。同樣的細(xì)胞類型在不同批次里測出來的nFeature中位數(shù)可能相差上千。如果你拿一個中位nFeature只有800的冷凍腫瘤樣本去套用PBMC的200~2500區(qū)間那你的前半段正常細(xì)胞就會被無情誤殺。更讓人頭疼的是細(xì)胞類型差異。我處理過一個中性粒細(xì)胞比例很高的樣本中性粒細(xì)胞nFeature中位數(shù)不到500按nFeature必須大于600的常規(guī)思路過濾整個中性粒細(xì)胞群幾乎被團滅?？墒侨思抑行粤＜?xì)胞本來就是短命、轉(zhuǎn)錄本少的細(xì)胞你把它們?nèi)珓h了后續(xù)分析還怎么做這就是精準(zhǔn)兩個字的分量過濾不是要套公式而是要理解你的數(shù)據(jù)里有哪些細(xì)胞、它們的指標(biāo)形態(tài)是怎樣的。4.2 用MAD法則和分布形態(tài)來定閾值那怎么定閾值才算靠譜我自己的習(xí)慣是兩步走先看分布再用分位數(shù)輔助判斷。第一步用核密度圖看nFeature和nCount的整體形態(tài)特別注意有沒有一個矮子峰# 畫核密度圖觀察nFeature_RNA的分布 plot(density(sce$nFeature_RNA), main nFeature density, col steelblue) # 觀察percent.mt的分布 plot(density(sce$percent.mt), main percent.mt density, col firebrick)如果nFeature_RNA的密度圖在左邊低基因數(shù)區(qū)域出現(xiàn)一個獨立的小峰或者峰型很扁那說明數(shù)據(jù)里確實混了一批低質(zhì)量的細(xì)胞。而percent.mt的密度圖如果出現(xiàn)一個拖得很長的尾巴或者主峰明顯右偏那線粒體高比例的細(xì)胞一定不在少數(shù)。第二步就是更定量的做法。很多人會采用中位數(shù)±3倍絕對中位差MAD來幫助確定閾值。它的邏輯跟平均值±3倍標(biāo)準(zhǔn)差類似但MAD對異常值更不敏感更適合單細(xì)胞這種厚尾分布的數(shù)據(jù)# 用中位數(shù)±3倍MAD估計合理范圍僅供初步參考 lower_bound - median(sce$nFeature_RNA) - 3 * mad(sce$nFeature_RNA) upper_bound - median(sce$nFeature_RNA) 3 * mad(sce$nFeature_RNA) cat(nFeature_RNA 參考范圍:, lower_bound, -, upper_bound, \n)這套方法的好處是閾值完全由你自己的數(shù)據(jù)驅(qū)動而不是生搬硬套別人的經(jīng)驗。但它同樣有局限性——如果你的數(shù)據(jù)里低質(zhì)量細(xì)胞占比太高MAD本身就會被它們拉著走導(dǎo)致算出來的范圍過寬。所以我的建議是MAD算出來的范圍用來做參考最終閾值還是回到密度圖上人工確認(rèn)一下看看切的位置是否正好落在主峰和低質(zhì)量小峰之間的洼地上。4.3 不同樣本類型的初始閾值參考為了讓你有個大概的手感我把自己處理過的幾類樣本的初始閾值整理成了表格。注意這只是第一輪粗過濾的參考值不代表所有項目通用樣本類型nFeature_RNA 參考下限nFeature_RNA 參考上限percent.mt 參考上限PBMC 外周血200~4004000~500010~15腫瘤組織新鮮解離300~5006000~800015~20冷凍組織200~3004000~600020~25中性粒細(xì)胞富集樣本100~2003000~400010~15snRNA-seq核測序500~8008000~100001~5表格里的數(shù)值跨度比較大因為每批數(shù)據(jù)的質(zhì)量波動實在太大。我一直強調(diào)一個觀點表格里的數(shù)字只是用來入門的不是你交論文時的最終答案。最終答案永遠(yuǎn)要靠你對自己樣本的理解和數(shù)據(jù)分布的實際觀察來回答。這也是精準(zhǔn)質(zhì)控和應(yīng)付差事式質(zhì)控的核心差別。5. 過濾完別急著開心怎么驗證質(zhì)控效果過濾代碼跑通了細(xì)胞數(shù)也減少了并不代表質(zhì)控做完了。很多人在這時候馬上進入NormalizeData、PCA、UMAP的下一個階段結(jié)果最后發(fā)現(xiàn)分群有問題才想起回頭找質(zhì)控的問題。我更建議你在過濾之后用幾分鐘做幾項驗證把異常在早期就暴露出來。5.1 過濾前后指標(biāo)的對比最直接的驗證就是對比過濾前后的指標(biāo)分布# 過濾前后指標(biāo)分布對比 VlnPlot(sce, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3) VlnPlot(sce_filtered, features c(nFeature_RNA, nCount_RNA, percent.mt), ncol 3)過濾后你應(yīng)該看到nFeature_RNA和nCount_RNA的低值尾被切掉了整體分布更集中percent.mt的高值尾巴被削平大部分細(xì)胞集中在一個相對正常的區(qū)間小提琴圖整體變得干凈不再有明顯的離群點。如果過濾后你的nFeature主峰反而被砍掉一大半或者percent.mt的分布幾乎沒有變化那說明閾值設(shè)得有問題要么誤傷了大群正常細(xì)胞要么壓根沒把低質(zhì)量細(xì)胞濾除干凈。5.2 檢查高表達基因列表是否合理這是一個非常容易被忽略但很聰明的檢查方法。正常情況下一個健康單細(xì)胞數(shù)據(jù)集中高表達基因應(yīng)當(dāng)是各類細(xì)胞類型的功能基因比如PBMC里應(yīng)該是MALAT1、B2M、各類核糖體蛋白基因等。但如果你的數(shù)據(jù)里混入了大量垂死細(xì)胞或破裂細(xì)胞高表達基因排行榜會被線粒體基因霸榜——因為線粒體基因在破損細(xì)胞里占比太高了。我有個習(xí)慣性操作質(zhì)控后跑一下每個細(xì)胞的top基因看看排名靠前的是不是有大量MT-基因。具體可以這樣# 檢查每個細(xì)胞中最顯著高表達的基因 # 取前500個高表達基因 top_genes - rownames(sce_filtered)[ order(Matrix::rowSums(sce_filteredassays$RNA$counts), decreasing TRUE) ] head(top_genes, 30)如果你看到head出來的30個基因里有大量MT-開頭的基因那就要警惕線粒體基因過濾可能還不夠干凈尤其在有核紅細(xì)胞、心肌組織等線粒體天然豐富的樣本里。至于要不要進一步往上提高percent.mt閾值同樣還是回到分布圖和生物學(xué)實際情況來決定。5.3 雙重體到底要不要處理過濾驗證階段最容易被忽略的其實是雙重體doublets。Seurat的質(zhì)控指標(biāo)可以篩掉低質(zhì)量細(xì)胞但對雙重體的效果非常有限——因為一個雙重體的nFeature通常并不是特別離譜它可能是兩個中等基因數(shù)的細(xì)胞粘在一起生成的nFeature值恰好落在正常范圍的高值區(qū)蒙混過關(guān)。要嚴(yán)格去雙就得用專門的工具。目前最常用的是DoubletFinder# 這個步驟建議在質(zhì)控完成、數(shù)據(jù)標(biāo)準(zhǔn)化之后運行 # 簡單用法示意 # library(DoubletFinder) # sce_filtered - NormalizeData(sce_filtered) # sce_filtered - FindVariableFeatures(sce_filtered) # sce_filtered - ScaleData(sce_filtered) # sce_filtered - RunPCA(sce_filtered) # sce_filtered - FindNeighbors(sce_filtered, dims 1:20) # sce_filtered - FindClusters(sce_filtered, resolution 0.5) # sce_filtered - RunUMAP(sce_filtered, dims 1:20) # # 然后基于聚類結(jié)果預(yù)測doublet # sweep.res - paramSweep(sce_filtered, PCs 1:20, sct FALSE) # sweep.stats - summarizeSweep(sweep.res, GT FALSE) # bcmvn - find.pK(sweep.stats) # ...DoubletFinder的原理是通過人工模擬doublet看哪些細(xì)胞在特征空間上與模擬的doublet最相似從而打分預(yù)測。這個工具能篩掉相當(dāng)一部分漏網(wǎng)的doublet但它也有誤殺風(fēng)險而且是基于聚類結(jié)果預(yù)測的如果前面聚類本身就不好預(yù)測的可靠性也會受影響。所以我通常建議在分群結(jié)果出來后、確定marker基因之前結(jié)合doublet注釋結(jié)果做一次排查盡量不要在只有幾十個細(xì)胞的支持下強行去雙。5.4 與后續(xù)分析流程的正確銜接質(zhì)控是起點不是終點。過濾干凈之后后續(xù)步驟必須銜接好# 完整的下游分析銜接 sce_filtered - NormalizeData(sce_filtered, normalization.method LogNormalize, scale.factor 10000) sce_filtered - FindVariableFeatures(sce_filtered, selection.method vst, nfeatures 2000) sce_filtered - ScaleData(sce_filtered) sce_filtered - RunPCA(sce_filtered, npcs 30) sce_filtered - RunUMAP(sce_filtered, dims 1:20) sce_filtered - FindNeighbors(sce_filtered, dims 1:20) sce_filtered - FindClusters(sce_filtered, resolution 0.5) # 后續(xù)差異分析 markers - FindAllMarkers(sce_filtered, only.pos TRUE, min.pct 0.25, logfc.threshold 0.25)這里我不展開講每個步驟的細(xì)節(jié)但要注意的是質(zhì)控越靠前后面全部步驟都能受益如果你在整合多個樣本最好先分別對每個樣本做各自的質(zhì)控再進入整合流程。把不同質(zhì)量的樣本混在一起再統(tǒng)一過濾往往會出現(xiàn)好的樣本被拉到差的樣本水平這種尷尬處境。6. 實戰(zhàn)中的特殊樣本那些容易翻車的質(zhì)控場景寫代碼部分的質(zhì)控大家都能跑通真正拉開差距的是處理那些非典型樣本。我在多個項目里踩過坑下面幾個場景幾乎每個做單細(xì)胞的人都可能遇到專門拿出來說。6.1 免疫細(xì)胞樣本的線粒體比例天然偏高如果你是做免疫細(xì)胞、尤其是組織浸潤免疫細(xì)胞的一定不要迷信Seurat官方教程的percent.mt 5。5%這個閾值來自PBMC外周血的樣本背景干凈、細(xì)胞狀態(tài)比較均一。但組織里的T細(xì)胞、巨噬細(xì)胞在解離過程中經(jīng)歷了機械剪切、酶消化應(yīng)激狀態(tài)下線粒體活性會升高percent.mt就很容易沖到10%~20%甚至更高。我處理過一批腫瘤浸潤T細(xì)胞樣本percent.mt的中位數(shù)就有14%如果按5%切直接砍掉一大半T細(xì)胞后面根本沒法分析。這種情況下我的建議是把閾值放到20%甚至25%同時結(jié)合nFeature和nCount來判斷——一個T細(xì)胞如果percent.mt是20%但nFeature有1200它大概率還是一條好漢但如果percent.mt 20%同時nFeature只有150那它大概率確實快死了。6.2 腫瘤樣本低質(zhì)量與惡性信號混雜的難題腫瘤組織是單細(xì)胞項目里最復(fù)雜的一類樣本。腫瘤細(xì)胞本身存在拷貝數(shù)變異、倍性異常、增殖活躍它們的nFeature和nCount往往很高高到如果用nFeature_RNA 5000一刀切會把真正的腫瘤細(xì)胞當(dāng)成doublet去掉。另一方面腫瘤微環(huán)境里又有大量浸潤的免疫細(xì)胞、基質(zhì)細(xì)胞各類細(xì)胞之間的nFeature跨度極大。在這種樣本里做質(zhì)控我的策略是在解離后先做一個偏保守的粗過濾確保把明顯的垃圾先清掉然后在聚類之后用inferCNV等工具確認(rèn)腫瘤細(xì)胞群單獨看它們的質(zhì)控指標(biāo)分布再做一輪針對性的二次檢視。這一步很關(guān)鍵——如果你在前期就下手太狠可能把一個重要的腫瘤亞群給做沒了后面再怎么注釋都找不回來。6.3 snRNA-seq核測序千萬不要當(dāng)成普通單細(xì)胞處理單細(xì)胞核測序snRNA-seq這幾年特別火尤其在冷凍組織中優(yōu)勢明顯因為核不容易受到解離損傷的影響。但它的質(zhì)控指標(biāo)跟普通單細(xì)胞完全不一樣核里面不含完整的細(xì)胞質(zhì)線粒體基因本來就極少percent.mt通常只有1%以下同時核轉(zhuǎn)錄本的捕獲效率更高nFeature中位數(shù)往往比完整細(xì)胞更高。如果有人拿著普通單細(xì)胞的標(biāo)準(zhǔn)去過濾snRNA-seq數(shù)據(jù)大概率會把數(shù)據(jù)搞壞。snRNA-seq的質(zhì)控重點是排除帶大量細(xì)胞質(zhì)殘留的完整細(xì)胞碎片和空的細(xì)胞核而不是糾結(jié)線粒體比例。另外snRNA-seq中低nFeature并不一定代表低質(zhì)量——有些靜止細(xì)胞的核轉(zhuǎn)錄本本來就少。如果后面還要結(jié)合ATAC數(shù)據(jù)做多組學(xué)分析質(zhì)控更要謹(jǐn)慎。6.4 污染與背景RNA的影響做過組織樣本的朋友應(yīng)該都有體會組織解離之后總有一些細(xì)胞碎片和游離RNA它們會被隨機包裹進GEM微滴。這種背景RNA污染會造成一個現(xiàn)象所有細(xì)胞的nFeature和nCount都被稀釋了——具體來說垃圾細(xì)胞帶著高背景RNA會讓整體數(shù)據(jù)出現(xiàn)一個低質(zhì)量團塊。Seurat的普通質(zhì)控對這種污染有一定的清理效果但未必能完全解決。如果背景污染比較嚴(yán)重我更推薦用CellBender或者SoupX這類專門的去污染工具在質(zhì)控之前先把背景RNA信號扣除掉。這個操作能讓后續(xù)閾值判斷更精準(zhǔn)但要注意這類工具對計算資源和參數(shù)設(shè)置有一定要求不是所有項目都必需。一般建議是如果聚類后所有細(xì)胞類型都帶有同樣一組高表達基因、UMAP圖分群模糊并且高表達基因列表里有明顯的環(huán)境RNA特征比如各種血紅蛋白基因才需要考慮去污染。7. 寫在最后的實操心得文章寫到這兒代碼和理論都講得差不多了。最后分享幾點只有自己上手跑過幾批數(shù)據(jù)才會總結(jié)出來的體會。第一質(zhì)控閾值不要追求一次到位。我自己現(xiàn)在跑一個新數(shù)據(jù)集很少會第一次就定出完美的閾值。通常是先按照分布圖設(shè)一個大概范圍過濾后看看細(xì)胞數(shù)和UMAP效果再回來微調(diào)。這個跑一遍—看結(jié)果—回來調(diào)的循環(huán)在單細(xì)胞項目里極其正常不要覺得麻煩。第二質(zhì)控日志一定要記。具體包括過濾前多少細(xì)胞、過濾后多少細(xì)胞、三個指標(biāo)的閾值各是多少、當(dāng)時看了哪些圖、為什么這么定。不要高估自己的記憶力。一篇論文從投稿到返修可能隔了大半年如果到時候?qū)徃迦藛柲隳愕倪^濾標(biāo)準(zhǔn)是怎么確定的而你連跑時的參數(shù)都忘了那場面會很尷尬。我現(xiàn)在每個項目都會留一個quality_control_notes.md文件隨手記錄當(dāng)時的判斷邏輯。第三不同樣本之間不要套用同一個閾值。哪怕兩個樣本來自同一個病人的同一類組織只要批次不同、建庫時間不同就應(yīng)該各看各的分布再分別定閾值。把兩個樣本合并成一個對象后統(tǒng)一過濾是新手最容易犯的錯誤之一它會直接導(dǎo)致高質(zhì)量的樣本被你誤傷。第四質(zhì)控的上限不只是不濾錯還要考慮下一步需要什么。比如你后面想做擬時序分析那G0/G1期的細(xì)胞就要盡量保留如果后面想做轉(zhuǎn)錄因子活性分析轉(zhuǎn)錄本種類不夠豐富的細(xì)胞可能支撐不起結(jié)論。每一種下游分析對數(shù)據(jù)質(zhì)量的要求都不同。質(zhì)控參數(shù)不是越嚴(yán)格越好而是越契合你的分析目標(biāo)越好。單細(xì)胞轉(zhuǎn)錄組質(zhì)控這條路說簡單也簡單就是看三個指標(biāo)、設(shè)三個閾值說復(fù)雜也復(fù)雜因為每個數(shù)據(jù)集都有它自己的脾氣。如果你把這篇內(nèi)容看懂了再回到自己的數(shù)據(jù)上去實踐幾輪應(yīng)該能明顯感受到精準(zhǔn)質(zhì)控和套參數(shù)跑流程之間的差距。希望這些經(jīng)驗?zāi)軒湍闵僮咭稽c彎路。