掃碼聯(lián)系客服
單細胞RNA測序數(shù)據(jù)基因表達不同轉(zhuǎn)化方法比較
研究背景
單細胞RNA測序從細胞層面研究腫瘤異質(zhì)性，獲得腫瘤組織內(nèi)部單個細胞的基因結構和基因表達狀態(tài)，同時捕獲細胞之間通訊和發(fā)育軌跡，相比RNAseq，對基因表達定量更加精準。目前主流技術是基于barcode的單細胞測序，將單細胞懸液樣品和帶有barcode的水凝膠珠子，通過微流體芯片，包裹在一個油滴之中。在油滴中進行逆轉(zhuǎn)錄后，每個單細胞的cDNA文庫，就帶上了獨一無二的barcode。由于每個細胞中RNA含量極低，對每個轉(zhuǎn)錄本加入UMI，可有效降低PCR擴增錯誤導致的計數(shù)誤差。基因x細胞的表達計數(shù)矩陣是單細胞RNA測序數(shù)據(jù)分析基礎，無論是聚類分析還是特定細胞亞群表達基因分析，或者細胞通訊分析，都需要輸入這個計數(shù)表格數(shù)據(jù)，且基于統(tǒng)計推斷的單細胞分析工具都有個假設前提是基因在不同細胞之間表達計數(shù)具有相似方差。
研究結論
來自德國海德堡大學的科研人員在模擬和真實單細胞基準數(shù)據(jù)基礎上，比較了4類對基因x細胞的表達計數(shù)矩陣進行轉(zhuǎn)化以獲得不同細胞之間動態(tài)相似方差的統(tǒng)計方法：delta、殘差、直接表達推測和因子分析。結果發(fā)現(xiàn)：一種簡單的計數(shù)增量然后log轉(zhuǎn)化加PCA分析，跟有著較好理論基礎的殘差方法、表達推測和因子分析等有著同樣甚至更加優(yōu)異表現(xiàn)。
研究結果
基礎理論
轉(zhuǎn)錄組測序數(shù)據(jù)中高表達基因往往具有較大的計數(shù)方差，不同實驗批次、測序數(shù)據(jù)質(zhì)量、分析流程差異都可能帶來較大干擾，導致分析結果偏離真實情況。目前RNA表達數(shù)據(jù)的標準統(tǒng)計方法對方差一致的數(shù)據(jù)表現(xiàn)最好，對于方差不齊的RNA數(shù)據(jù)，一般采用抽樣分布建模來預估分布參數(shù)或者數(shù)據(jù)轉(zhuǎn)化以獲得方差齊性進行預處理。建模參數(shù)推斷往往耗時長且需要大量計算資源，更有效的方法是對計數(shù)進行轉(zhuǎn)化獲得方差齊性。對于單細胞RNA測序數(shù)據(jù)，UMI標記轉(zhuǎn)錄本計數(shù)比較符合負二項分布，各種轉(zhuǎn)換方法均在此基礎上進行。
轉(zhuǎn)化方法
delta
通過簡單的非線性轉(zhuǎn)化，如g(y) = log (y + y0)，其中y代表基因在單個細胞中的原始計數(shù)，y0為固定值，使得方差在動態(tài)范圍內(nèi)更加接近，其中y0的選擇對結果準確性很重要。單細胞RNA測序數(shù)據(jù)中每個細胞UMI總數(shù)受采樣效率和細胞大小影響，在進行方差穩(wěn)定轉(zhuǎn)化前，需要對每個細胞確定一個放縮因子s，常規(guī)方法是讓所有細胞的放縮因子之和接近于1，所以s=單個細胞UMI總數(shù)/L，其中L表示所有細胞的UMI總數(shù)之和的平均值，然后再進行方差穩(wěn)定轉(zhuǎn)化，即g(y) = log(y/s + y0)。Seurat中L設置為固定值10,000，相當于y0=0.5，α=0.5，更接近于真實數(shù)據(jù)中觀察到的過度離散狀態(tài)(overdispersion)，CPM計數(shù)本質(zhì)是L=106，相當于y0=0.005，α=50，比典型單細胞表達數(shù)據(jù)的離散度大了兩個數(shù)量級。對于模擬數(shù)據(jù)，該類別分別基于acosh轉(zhuǎn)化、y0=1、y0=1/4(α)、CPM計數(shù)、log轉(zhuǎn)化后選擇表達高變化基因（HVG）或者zscore進行額外處理等方法進行評估。
殘差
前面提到，最接近UMI標記單細胞表達計數(shù)的統(tǒng)計模型是負二項分布，因此可以采用廣義線性模型（GLM）對實際表達數(shù)據(jù)進行回歸分析，以期達到Pearson殘差最小化，而負二項分布參數(shù)均值通過每個基因擬合對應的斜率和截距以及單個基因放縮因子的log變換得到。對于模擬數(shù)據(jù)，該類別分別基于剪切和未剪切的Pearson殘差、隨機分位數(shù)殘差、選擇表達高變化基因（HVG）或者zscore進行l(wèi)og轉(zhuǎn)化后處理等方法進行評估。
基因表達推測
該方法的基本思想是以觀察到的基因計數(shù)為前提，直接推斷該基因表達狀態(tài)，貝葉斯模型是這個問題的典型應用。代表方法有基于全貝葉斯模型的Sanity，它另外包括兩種細分方法，一種通過計算基因表達log轉(zhuǎn)化的均值和方差來衡量細胞之間的距離，找出距離最近的k個細胞，另一種Santiy MAP是返回后驗最大值。Dino擬合混合負二項分布，返回隨機樣本后驗值，Normalisr返回每個計數(shù)的最小均方誤差估計。對于模擬數(shù)據(jù)，在上述提到的四種方法中進行了評估。
因子分析
該方法不對計數(shù)進行轉(zhuǎn)化，而是將多維單細胞數(shù)據(jù)降到低維空間，在低維空間中，同樣符合負二項分布的UMI計數(shù)進行因子分析。包括GLM PCA和NewWave，后者是前者的優(yōu)化版本，對以上兩種方法分別進行模擬數(shù)據(jù)評估。
方法比較
對包裹了來自相同RNA的等分試樣的液滴均勻溶液，分別比較原始計數(shù)、delta、剪切Pearson殘差、Sanity MAP等轉(zhuǎn)化方法后PCA前兩個主成分分布、均值-方差關系分布，圖1a中每個點代表一個液滴。從分析結果可以看出，對于delta方法，放縮因子仍代表強方差變異，不同級別放縮因子液滴明顯分群，這也是delta方法的致命缺點，將具有較大計數(shù)和較大放縮因子的液滴和較小計數(shù)和較小放縮因子的液滴分到了一起，而這并不符合負二項分布的均值-關系分布。相比之下，GLM殘差和Sanity MAP都能夠?qū)σ旱芜M行很好的混合。
對取自造血細胞的10x數(shù)據(jù)2597個基因表達分析結果顯示，delta方法還有一個缺陷是對低表達基因(表達均值<0.1)轉(zhuǎn)化后方差幾乎為0，不能體現(xiàn)低表達基因?qū)Y果的影響。相比之下，經(jīng)過GLM殘差轉(zhuǎn)化的方法顯示低表達基因?qū)档娜踝兓?，除非極低表達基因，而Sanity MAP轉(zhuǎn)化后表達是多樣的，顯示這種方法不直接與穩(wěn)定方差相關。
SFTPC基因在肺上皮細胞中呈現(xiàn)雙峰表達模式，挑選該基因的單細胞表達分布來衡量3類轉(zhuǎn)化方法差異。結果顯示：delta和Sanity MAP方法均能很好的體現(xiàn)SFTPC的不同表達水平細胞亞群水平，而Pearson殘差對于高表達細胞亞群的方差相比低表達細胞亞群的方差不能更大程度的縮小，所以高表達亞群細胞會被更加細分，而低表達細胞亞群會更加聚集。這會影響此類在不同細胞中差異表達基因的可視化，還有標記基因檢測和細胞聚類和分類。
圖1.不同轉(zhuǎn)化方法對比
基準數(shù)據(jù)測評
方法評估取決于分析目的，如果為了比較細胞類型特異性標記基因檢測，可通過圖1c看細胞計數(shù)分布，而本文評估目標定位于基于低維數(shù)學結構（聚類、軌跡、分支和組合）理解細胞類型和狀態(tài)的多樣性，k-NN圖體現(xiàn)了基本功能結構，通過計算k-NN圖之間的重疊度以評估結果的一致性，值越大說明重疊度越高，準確性越好。
本文主要通過以下方法構建真實數(shù)據(jù)集：
從GEO數(shù)據(jù)庫中下載10個10X轉(zhuǎn)錄組數(shù)據(jù)集，每份數(shù)據(jù)一分為二，理論上兩份數(shù)據(jù)之間k-NN圖重疊度較好，一致性較高。
基于已發(fā)表的4種數(shù)據(jù)模擬框架生成模擬數(shù)據(jù)，且本文對其中一種模擬方法進行了修正，對模擬數(shù)據(jù)不同轉(zhuǎn)化方法生成k-NN圖與參考k-NN圖比較重疊度。
對來自mcSCRB和Smart-seq3的5個深度測序數(shù)據(jù)集，對其進行下采樣到10x深度，分別比較深度測序k-NN圖與下采樣數(shù)據(jù)k-NN圖的重疊度。
對以上評估數(shù)據(jù)集分別進行PCA分析、k-NN圖繪制、k-NN與參考k-NN重疊度計算，涉及包含在四大類中的22種轉(zhuǎn)化方法，包括過度離散值固定為0、0.05和基因特異性評估，4-8個PCA亞群參數(shù)設置，k=10、50、100的最佳鄰接細胞重疊計算。對k=50的最佳鄰接圖形重疊度計算結果顯示：在三種數(shù)據(jù)集中，delta方法都得到了較高的k-NN圖重疊，特別在一分為二的數(shù)據(jù)集上體現(xiàn)比較明顯，而在模擬框架數(shù)據(jù)集上，不同轉(zhuǎn)化方法差異不明顯，但針對不同的模擬方法還是有很大的差異。下采樣數(shù)據(jù)集代表了真實數(shù)據(jù)情況，更接近于10x轉(zhuǎn)錄組數(shù)據(jù)，結果顯示不同轉(zhuǎn)化方法k-NN重疊度差異不大。
之前文章研究發(fā)現(xiàn)，在random walk模擬方法產(chǎn)生的數(shù)據(jù)集上，Sanity是識別細胞k-NN的最佳方法，本文研究人員在先將多維單細胞數(shù)據(jù)降到低維，然后利用delta和殘差轉(zhuǎn)化方法分析細胞k-NN時取得同樣優(yōu)異表現(xiàn)。通過PCA亞群數(shù)量和k-NN重疊度分布曲線可以看出，數(shù)據(jù)空間維度大小對性能評估影響很大，隨著空間維度增大（>1000），k-NN重疊度下降，特別是在兩種模擬方法數(shù)據(jù)集上體現(xiàn)比較明顯，在維度較低(<100)的一致性評估數(shù)據(jù)集和下采樣數(shù)據(jù)集中也有所體現(xiàn)。
圖2.基準數(shù)據(jù)測試結果
對10,064個輔助T細胞的10x轉(zhuǎn)錄組測序數(shù)據(jù)進行方差穩(wěn)定轉(zhuǎn)化和k-NN計算，不同方法耗費的CPU時間差距比較大，基于表達推測的方法整體比較耗時，其中Sanity distance時間最長，由于要計算細胞之間的歐式距離，GLM PCA和NewWave耗時次之，delta方法最快，保留了數(shù)據(jù)的稀疏性，除了Pearson殘差，其它殘差方法速度跟delta方法一樣快。隨著細胞個數(shù)增加，所有方法的消耗時間同步增加，大部分方法的增長曲線斜率為1，而Sanity和GLM PCA的增長斜率為1.5。
圖3.計算資源消耗統(tǒng)計
所有轉(zhuǎn)化方法綜合對比
在細胞亞群分類方面，log(y/s+1)轉(zhuǎn)化始終優(yōu)于其它方法，而且當PCA降維到合適的目標維度后，log(y/s+1)轉(zhuǎn)化在所有三個基準中都比更復雜的基于表達推斷的轉(zhuǎn)換表現(xiàn)得更好。額外的處理步驟，如重新縮放移位對數(shù)的輸出，選擇hvg或使用zscore使所有基因的方差相等等并沒有太大改善識別細胞k-NN結果。對于皮Pearson殘差和acosh變換，設置α > 0是有用的，但與使用固定值如0.05相比從輸入數(shù)據(jù)估計α值，性能沒有明顯的提高。隨著每個細胞測序深度增加，所有方法通常都有更好的k-NN重疊度，顯示出更好的分析性能，主要是由于隨著測序深度增加，隨機采樣帶來的噪音就會減小。
圖4.所有轉(zhuǎn)化方法評估結果
總結
本文從三類基準數(shù)據(jù)，包括10個單細胞測序數(shù)據(jù)一分為二的數(shù)據(jù)集、4種不同模擬方法產(chǎn)生的數(shù)據(jù)集以及從高深度單細胞RNA測序數(shù)據(jù)下采樣得到深度10的數(shù)據(jù)集，評估來自于4大類（delta、殘差、表達推測、因子分析）共22種細分計數(shù)轉(zhuǎn)化方法對細胞亞群結構的分析準確度，研究發(fā)現(xiàn)采用簡單的log(y/s+1)轉(zhuǎn)化，加上數(shù)據(jù)降維到低維度空間就可以跟有著較好理論基礎的復雜模型一樣的優(yōu)異表現(xiàn)。
