PUseqClust:一種RNA-seq 數(shù)據(jù)聚類分析方法*

石險峰,劉學軍,張 禮

1(南京航空航天大學 計算機科學與技術學院,江蘇 南京 211106)

2(南京林業(yè)大學 信息科學技術學院,江蘇 南京 210037)

通訊作者:劉學軍,E-mail:xuejun.liu@nuaa.edu.cn

第二代測序技術又稱下一代高通量DNA 測序技術(NGS),與基因芯片相比,具有通量高、速度快、成本低的優(yōu)點.該技術的提出,使基因組學、基因表達學和表觀遺傳學產(chǎn)生了革命性的改變,為人們解決存在的生物問題和全面透徹地分析物種的基因組、轉錄組提供了重要工具[1].其中,基于深度測序的RNA-seq 測序技術被廣泛應用于轉錄組研究.RNA-seq 可以對某一特定物種的轉錄組和基因組序列進行全面分析,并獲得特定生理或者病理狀態(tài)下的所有信息,極大地提高了測序速度并且降低了成本.其主要思想是:通過序列對比,將讀段(read)定位到參考基因組或者轉錄組上,獲得量化的表達值[2,3].隨著測序技術的不斷發(fā)展,生物數(shù)據(jù)庫中存儲了海量的RNA-seq 讀段數(shù)據(jù),通過讀段計數(shù),可量化基因表達水平.而如何處理與分析這些基因表達值,并從中獲得有用的生物學意義,已經(jīng)成為當前熱門的研究方向.

目前,許多RNA-seq 測序實驗除了對某一特定條件下的物種進行分析,還在多條件多重復實驗下進行測序,如不同的溫度[4]、不同的組織[5]、不同的時間點[6]等條件.在自然界環(huán)境下,外界的條件會快速并且隨機改變,所以為了適應這些變化,基因在不同條件下會表現(xiàn)出不同的表達水平[7].通過分析這些不同表達模式的數(shù)據(jù),我們可以獲得客觀數(shù)據(jù)中所包含的生物意義.由于生物的性狀是多個基因共同調(diào)控的結果,而聚類分析根據(jù)不同的表達模式,將基因聚到不同的類簇.通常假設功能相關的基因會聚到同一類簇,由此可揭示基因的未知生物功能.因此,聚類分析對于處理這些多條件多重復性的RNA-seq 數(shù)據(jù)顯得尤為重要.但是當前大部分聚類分析主要應用于基因芯片數(shù)據(jù),應用到RNA-seq 數(shù)據(jù)的聚類研究相對較少.本文旨在提出一種適用于多條件多重復性的RNA-seq 數(shù)據(jù)的聚類分析框架,能夠有效處理RNA-seq 數(shù)據(jù)存在的各種噪聲和偏差,獲得更具生物學意義的聚類結果.

RNA-seq 數(shù)據(jù)的聚類分析研究工作雖已有一定進展,但仍存在不足.如文獻[5]使用了K-means 算法對4 個不同組織下的表達模式進行了簡單聚類并做后續(xù)分析,但實際上這些傳統(tǒng)算法,如分層算法、K-means 算法、自組織映射(SOM)算法,都是基于啟發(fā)式的算法,難以比較各自的優(yōu)劣,并且沒有一個特定的原則可以確定最優(yōu)的類簇個數(shù),給聚類分析增加了不少困難[8-10].一些基于模型的聚類算法應用到RNA-seq 數(shù)據(jù)上,并獲得了較好的聚類結果.這些方法中,大部分使用泊松分布對RNA-seq 數(shù)據(jù)進行建模,以減少讀段非均勻分布的影響[11-13].然而實驗結果表明,RNA-seq 數(shù)據(jù)與泊松模型相比較具有更高的差異性,即數(shù)據(jù)的方差大于數(shù)據(jù)的均值,這將導致過離散現(xiàn)象[14].所以越來越多的模型采用負二項分布模擬讀段數(shù)據(jù),以便更好地處理RNA-seq 數(shù)據(jù)的過離散現(xiàn)象.MBclust[15]采用了負二項式分布模擬讀段計數(shù)進行聚類,通過估計數(shù)據(jù)的散度系數(shù),獲得了相比K-means 等傳統(tǒng)方法更好的聚類結果.

現(xiàn)有的這些聚類方法直接采用讀段計數(shù)對基因表達水平進行聚類分析,這類方法存在一定不足.除了實驗中cDNA 文庫的制備、cDNA 中CG 堿基含量過高等因素造成了讀段在參考序列上非均勻分布外,RNA-seq 讀段非均勻分布另外一個主要原因是真核生物普遍存在選擇性剪切(AS)現(xiàn)象[16],即:一個基因的多個共享外顯子,每次選擇不同的外顯子進行組合形成多種異構體(isoform),導致由較多異構體共享的外顯子上往往讀段分布較多,而較少異構體共享的外顯子上讀段分布較少.現(xiàn)有的聚類方法中較少考慮到讀段的這種非均勻分布特性,從而降低了聚類的準確性.此外,直接采用讀段計數(shù)進行聚類分析的方法需要處理技術性、生物性等多種不確定性,如果僅采用一個模型,難以全面模擬各種不確定性.比如,MBclust 方法采用負二項分布模擬了讀段計數(shù)的過離散特性并且考慮了生物重復性影響,但是將聚類中心固定為一個數(shù)據(jù)點,忽略了聚類中心的不確定性,降低了聚類結果的可靠性.

本文主要針對現(xiàn)有方法的不足,采用兩步方法對RNA-seq 數(shù)據(jù)進行聚類分析:首先,充分考慮讀段數(shù)據(jù)固有的各種噪聲和偏差,計算出基因的表達水平及其技術性不確定性;然后,在考慮到基因表達水平生物性不確定性情況下,對所獲得的基因表達水平進行聚類分析.以往的研究結果表明:通過適當?shù)慕y(tǒng)計分析處理相關噪聲數(shù)據(jù),對于獲得有生物學意義的分析結果具有重要意義[17,18].通過概率模型,實驗產(chǎn)生的技術性不確定性能夠被融入模型中參與估計,這使得這些方法對于噪聲更加魯棒[19].在基因芯片數(shù)據(jù)分析中已證明:如果在表達水平的后續(xù)分析中考慮其不確定性,能夠獲得更具生物意義的結果[20,21].先前的工作中,我們設計了PGSeq[22]模型對RNA-seq 數(shù)據(jù)進行基因表達水平計算,該模型考慮了讀段數(shù)據(jù)中由各種噪聲引起的偏差,選擇性剪接是導致讀段非均勻分布和多源映射問題的一個重要原因.相比其他聚類算法,我們在基因表達水平估計步驟中,著重考慮了這些問題對表達水平計算準確性的影響,故對造成讀段計數(shù)非均勻分布的各種原因考慮較為完善.在PGSeq模型的基礎上,本文提出了PUseqClust(propagating uncertainty into RNA-Seq clustering)聚類框架,考慮了基因在不同條件、不同重復樣本下的表達水平的相關性,采用多維拉普拉斯方法獲得基因表達水平的不確定性,并將計算結果傳遞到混合t 分布聚類模型[23]中,實現(xiàn)對RNA-seq 數(shù)據(jù)的聚類分析.混合t 分布聚類模型使用了具有魯棒性的t 分布模型,并且考慮了生物性重復實驗數(shù)據(jù)的不確定性.本文在模擬數(shù)據(jù)集和3 個真實數(shù)據(jù)集上,驗證了PUseqClust 的聚類性能.

1 方 法

圖1 顯示了PUseqClust 方法的流程圖.輸入數(shù)據(jù)為RNA-seq 的原始讀段序列,通常以FASTA 或者FASTQ兩種格式存儲[24],經(jīng)過多步處理,獲得最終聚類分析結果.首先進行數(shù)據(jù)預處理,選擇轉錄組序列作為參考序列,利用Bowtie2[25]進行讀段定位,獲得讀段計數(shù),然后,利用PGSeq 方法模擬基因讀段的隨機產(chǎn)生過程;其次,采用多維拉普拉斯方法[26]獲得基因表達水平及其不確定性;再次,采用一元方差分析(one-way ANOVA)獲得差異基因,為聚類分析進行簡單過濾,為聚類分析過濾掉沒有明顯變化模式的無表達或穩(wěn)定表達基因;最后,將基因表達水平及其不確定性傳遞到混合t 分布聚類模型中,進行聚類分析,獲得聚類結果.

Fig.1 The flowchart of PUseqClust method圖1 PUseqClust 方法流程圖

1.1 模擬讀段隨機產(chǎn)生過程

PGSeq 方法[22]利用泊松-伽馬分布模擬每個外顯子讀段計數(shù)的分布,獲得了較為準確的基因以及異構體表達水平.我們采用PGSeq 模型對外顯子讀段的隨機產(chǎn)生過程進行了模擬,以便獲得基因表達水平重要參數(shù)的隨機特性.設在條件c的第r個重復樣本上基因g的外顯子i讀段計數(shù)為ygicr,由于數(shù)據(jù)都是逐個處理每個基因,因此小標g在后續(xù)大部分公式中被省略.yicr表示包含外顯子i的剪接異構體的歸一化讀段計數(shù)之和,.其中,

·wcr表示在條件c的第r個重復樣本的讀段總數(shù);

·li表示外顯子i的長度;

·Mik取值為0 或1,Mik取值為1 時,表示剪接異構體k包含外顯子i;

·ticrk表示在條件c的第r個重復樣本上剪接異構體k包含的外顯子i的讀段計數(shù),并假設其服從泊松分布:ticrk～Pois(βiαcrk),其中:βi表示外顯子i的偏差特性,為解決過離散問題,模型假設βi服從Gamma分布βi～Ga(a,b),a表示形狀參數(shù),b表示尺度參數(shù);αcrk表示在條件c的第r個重復樣本上剪接異構體k的表達水平比率.

1.2 基因表達水平及不確定性計算

得到模型參數(shù)的解后,在條件c的第r重復樣本上剪接異構體k的表達水平ticrk可以寫成:

其中,NB表示負二項分布.由此推出在條件c的第r個重復樣本上剪接異構體k表達水平的期望為

對于同一個基因,由于a,b為各條件共享,則在不同條件下,這個基因異構體表達水平期望的隨機性由αcrk的隨機性決定.

考慮到各個條件重復實驗下αcrk之間的相關性,我們首先采用多維拉普拉斯方法從公式(1)近似獲得向量α={αcrk}的分布.對公式(1)進行泰勒展開,如下:

其中,

假設基因的表達水平由對應的剪接異構體表達水平之和表示,即,得到的基因表達水平服從高斯分布.均值和方差如下:

1.3 差異基因檢測

我們采用一元方差分析進行差異基因檢測.方差分析(analysis of variance,簡稱ANOVA)是數(shù)據(jù)分析中一種常見的統(tǒng)計模型,探索因變量與自變量之間關系,用于多個樣本均數(shù)差別的統(tǒng)計假設檢驗.一元方差分析即為探索一個自變量對于因變量的觀察值影響.零假設H0為基因未發(fā)生差異表達,即對于n個k維輸入樣本表達均值都相等.為了計算統(tǒng)計顯著性,若H0成立,當總偏差平方和SST 固定不變時,檢驗統(tǒng)計量取為

其中,SSE為組內(nèi)偏差和,SSA為組間偏差平方和.對于給定顯著水平a,查找F分布臨界值F函數(shù)Fa,得到p值:若p值大于a,則接受零假設;否則,拒絕零假設.

對于采樣得到的基因表達水平,我們對于每個基因在各個條件的重復樣本下分別提取100 個正樣本數(shù)據(jù),對這些數(shù)據(jù)進行一元方差分析,設置顯著水平閾值,識別顯著差異表達的基因.

1.4 聚類分析

對于顯著差異基因,我們采用先前提出的基于不確定性的混合t 分布聚類模型PUMA-CLUSTII[23]進行聚類分析.PUMA-CLUSII 采用具有魯棒性的混合學生t 分布進行聚類,由于這是一種重尾分布,能更好地適應異常值,并且模型考慮了生物性重復不確定性及技術性不確定性,并能自動優(yōu)化到最優(yōu)類簇個數(shù).

假設在條件c的第r個重復樣本上基因g的表達水平服從高斯分布,其中,xgcr表示真實基因表達水平,表示基因g技術性不確定性.真實基因表達水平xgcr在同一條件上假設服從高斯分布,wgc表示條件c下基因均值,ηg表示基因表達水平生物上的不確定性.對于每一個基因引入一個隱藏變量ug,則t 分布可以寫成高斯分布和Gamma 分布的卷積形式:

其中,在第k類簇上,μk和Σk分別表示均值和協(xié)方差,vk表示自由度.則均值wg的概率為

假設對于每個基因所有條件下ηg共享,且服從Gamma 分布:

此時,模型隱藏變量hg=(xg,wg,ηg,zg,ug),模型參數(shù)θ=({μk},{Σk},{vk},{αk},{βk},{πk}).

PUMA-CLUSTII 方法利用最大似然法和變分EM 算法對模型進行求解,該方法采用了MMLP 準則[27]自動確定最優(yōu)聚類數(shù).

2 數(shù)據(jù)集

由于無法獲得真實數(shù)據(jù)集中的準確基因表達模式的聚類結果,我們分別使用了模擬數(shù)據(jù)集和3 個真實數(shù)據(jù)集對所提出的聚類方法進行驗證.

2.1 模擬數(shù)據(jù)集

本文研究的聚類分析方法根據(jù)基因在多個不同實驗條件下的表達模式進行聚類,由于真實的RNA-seq 數(shù)據(jù)集缺少大規(guī)模已知聚類結果的基因表達水平數(shù)據(jù),故模擬數(shù)據(jù)模擬生成已知的特定基因表達模式,作為真實聚類標簽.我們的聚類方法對具有相似表達模式的基因進行聚類,不針對現(xiàn)實中哪一種具體的生物學波動模式,模擬數(shù)據(jù)生成方法見文獻[28]以及我們先前的工作[20,23].本文在模擬數(shù)據(jù)中人為定義了6 種不同的用數(shù)學函數(shù)描述的表達模式,且增加了一種噪聲數(shù)據(jù)來驗證算法的魯棒性.本文利用PGSeq 模型生成模擬數(shù)據(jù)集[22],對我們所提出的聚類方法進行初步評價.模擬數(shù)據(jù)集共包含7 個條件,每個條件下包含3 個重復實驗數(shù)據(jù),將各個條件下的基因差異表達水平看做已知類簇的基因表達模式,這樣的模擬數(shù)據(jù)集可以驗證聚類算法的準確性.

本文選取人類的700 個基因作為模擬數(shù)據(jù),并分成7 組,其中,最后一組作為純隨機噪聲組,用于模擬實際數(shù)據(jù)中難以聚到任一類簇中的基因表達模式.對于給定一個基因,外顯子i的偏差特性βi的分布可從MAQC 數(shù)據(jù)集[29]中的HBR 樣本實驗數(shù)據(jù)中獲得,即得到Gamma分布參數(shù)a和參數(shù)b的值.對于給定基因g,模擬數(shù)據(jù)產(chǎn)生過程如下:

(1)對于外顯子i,從伽馬分布Gamma(a,b)中采樣得到βi;

(2)變量A的值從均勻分布U(0,5)中采樣到.設C為總條件個數(shù),p為基因g所屬的數(shù)據(jù)組,然后按照如下方法獲得logαgck的值.

· 對于第1 組～第4 組:

· 對于第5 組:

· 對于第6 組:

· 對于第7 組:

由此生成具有不同模式的logαgck,然后從高斯分布N(logαgck,log(αgck)/10)重新采樣獲得含有噪聲的.在每個條件下,從高斯分布重復采樣3 次,獲得帶生物性噪聲的重復αgcrk;

(3)從泊松分布Pois(wcrliMikβiαgcrk)中采樣獲得剪接異構體k在條件c的第r重復樣本上的外顯子i的讀段計數(shù)xicrk.其中,wcr為在條件c的第r重復樣本上的讀段總數(shù);li為外顯子i的長度;Mik為外顯子i與剪接異構體k之間的關系,值為1 表示剪接異構體k包含外顯子i;

(4)外顯子i在條件c的第r重復樣本上讀段計數(shù)yicr為包含外顯子i的所有剪接異構體讀段計數(shù)xicrk之和,基因g在條件c的第r重復樣本上讀段計數(shù)zgcr為所有外顯子的讀段計數(shù)yicr之和.最終生成的對數(shù)基因讀段計數(shù)歸一化后數(shù)據(jù)如圖2 所示.由圖2 可以看出:模擬數(shù)據(jù)中前6 組基因分別具有不同的表達模式,而第7 組基因沒有一致的表達模式,其代表了數(shù)據(jù)中存在的隨機噪聲,用以測試聚類方法的魯棒性.

Fig.2 The gene expression patterns of the seven clusters in the simulated data圖2 模擬數(shù)據(jù)的7 類基因表達模式

2.2 真實數(shù)據(jù)集

(1)Zhao 等人[30]采用人類激活T 細胞數(shù)據(jù)集對RNA-seq 和基因芯片在基因表達水平計算方面進行了對比研究.利用Illumina HiSeq 2000 platform 得到在0hr,2hr,4hr,6hr,24hr,72hr 這6 個時間點上的T 細胞RNA 數(shù)據(jù),其中每個時間點上包含兩個重復性實驗.本文采用一元方差分析檢測出該數(shù)據(jù)集的1 963個差異基因,并對這些基因進行聚類分析,驗證我們所提出的聚類算法的性能;

(2)?ij? 等人[31]對人類輔助T 17(Th17)細胞進行表達水平測量,利用Illumina HiSeq 2000 platform 得到了0.5hr,1hr,2hr,4hr,6hr,12hr,24hr,48hr,72hr 這9 個時間點上Th17 細胞RNA 數(shù)據(jù),其中,每個時間點上包含3 個重復性實驗.與前一個數(shù)據(jù)集處理相同,采用一元方差分析獲得了2 060 個差異基因進行后續(xù)聚類分析;

(3)GEO accession 為GSE90053 的數(shù)據(jù)集為人類多功能干細胞(pluripotent stem cell,簡稱PSC)數(shù)據(jù)集,利用Illumina HiSeq platform 2500 得到了0d,2d,8d,10d,11d,14d,21d,28d 這8 個時間點上PSC RNA 數(shù)據(jù)集,每個時間點上包含3 個重復性實驗.最終得到1 448 個差異基因進行后續(xù)聚類分析.

3 實驗結果與討論

為了驗證本文提出的方法在聚類分析方面的性能,我們使用PUseqClust 方法和其他聚類分析方法在模擬數(shù)據(jù)集和真實數(shù)據(jù)集進行聚類分析,并對比了分析結果.本文采用了兩種已有的聚類方法進行對比實驗.一種是沒有考慮表達不確定性的標準高斯混合聚類方法Mclust[32];另一種是基于負二項分布的聚類方法MBclust[15].

這兩種方法均采用RPKM 計算基因表達水平.Mclust 和MBclust 方法分別包含在R 軟件包Mclust 和MBCluster.Seq 中.

3.1 模擬數(shù)據(jù)集實驗

對于模擬數(shù)據(jù)集,我們首先在無隨機噪聲的前6 組數(shù)據(jù)上進行聚類分析,然后加入第7 組噪聲數(shù)據(jù),對聚類方法的魯棒性進行對比.在模擬數(shù)據(jù)集上,由于各個基因所屬類簇已知,所以本文在標準互信息NMI(normalized mutual information)、敏感度和特異度這3 個方面驗證各種聚類方法的性能.

靈敏度代表所有成對基因(屬于同一類簇的基因)被劃分到同一類簇的的概率,靈敏度越高,代表算法對來自同一類簇基因的識別度越高.特異度代表所有非成對基因(不屬于同一類簇的基因)被劃分為不同類簇的概率,特異度越高,代表算法對非同一類簇基因的識別度越高[15].互信息是信息論的一種信息度量,評估一個隨機變量中所包含的另外的一種隨機變量的信息量,或者認為是一個隨機量由于已知另一個隨機變量而減少的不確定性.這里我們將其作為量化真實聚類劃分和已知聚類劃分的共有信息量的程度.我們利用文獻[33]中的方法計算互信息的值,并且將互信息歸一化至0～1 之間的標準互信息.NMI 值越接近1,表示真實聚類劃分和已知聚類劃分的相關性越高,即聚類精度越高;反之,越接近零聚類精度越低.敏感度表示來自同一類簇的成對基因被分到同一類簇的概率,特異度表示來自不同類簇的成對基因被分到不同類簇的概率.這兩個值越高,表示聚類性能越好.

圖3 顯示了模擬數(shù)據(jù)集中無噪聲組和包含噪聲組兩種情況下的NMI 值、敏感度和特異度.當無噪聲組時,對于NMI 值和敏感度,PUseqClust 方法優(yōu)于Mclust 方法,并且PUseqClust 和Mclust 都獲得了比MBclust 更準確的聚類結果;當含有噪聲組時,各種方法的性能都有所下降,但PUseqClust 仍然獲得了較為準確的聚類結果.對于特異度,圖中顯示:PUseqClust 無論在含有噪聲組或者不含有噪聲組里,都比其余兩種方法性能優(yōu)越.實驗結果證明:PUseqClust 方法采用兩步方法進行聚類分析,考慮了技術上和生物上表達水平的不確定性,對噪聲的處理能力更強,因而獲得了較準確的聚類結果.

Fig.3 The NMI,sensitivity,specificity of different methods on the simulated data圖3 模擬數(shù)據(jù)集上各種聚類方法對應NMI 值、敏感度、特異度

3.2 真實數(shù)據(jù)集實驗

對于真實數(shù)據(jù)集,因真實的類簇個數(shù)和聚類劃分未知,PUseqClust 方法采用了MMLP 準則[27]自動確定最優(yōu)類簇個數(shù),MClust 使用貝葉斯信息準則BIC(Bayesian information criterion)值[34]計算最優(yōu)類簇個數(shù),而MBclust方法的實現(xiàn)軟件沒有提供最優(yōu)類簇個數(shù)的確定方法,因此,我們只考察PUseqClust 和MClust 獲得的最優(yōu)類簇個數(shù).本文利用DAVID[35]對各種方法獲得的聚類結果進行GO(go ontology)注釋富集分析[36].GO 注釋富集分析通過尋找基因所屬的特定GO 富集功能目錄直接評估聚類結果的生物意義,其分析并找出在統(tǒng)計上顯著富集的GO 功能目錄,通過將差異基因做GO 富集分析,可以把基因按照不同的功能進行分類,達到對基因進行注釋和分類的目的.對于某個類簇的基因,GO 功能目錄越多,意味著這個類簇內(nèi)的基因在生物意義上相關程度越高.聚類分析結果應該盡可能增加相關程度相對較高的類簇數(shù)量,減少相關程度相對較低的類簇數(shù)量.富集計算結果在DAVID 上是以改進的Fisher exact test 評估,采用p值表示基因聚類結果的生物意義.GO 生物過程水平5(biological process level 5)設置閾值計數(shù)水平為5、p值為0.05 時,含有GO 富集功能目錄的類簇被認為是GO富集類簇,因此,我們設定閾值相關計數(shù)水平為5、p值為0.05,然后獲得所有的GO 富集功能目錄個數(shù),以便評價整個聚類方法的性能.

在T 細胞數(shù)據(jù)集的分析中,PUseqClust 和Mclust 方法獲得最優(yōu)類簇個數(shù)分別為17,18.Th17 細胞數(shù)據(jù)集PUseqClust 和Mclust 方法獲得最優(yōu)聚類個數(shù)分別為15,8.PSC 數(shù)據(jù)集中分別為20,21.

圖4～圖6 分別顯示了T 細胞數(shù)據(jù)集、Th17 細胞數(shù)據(jù)集和PSC 數(shù)據(jù)集上各種方法在GO 富集功能目錄指定數(shù)量范圍內(nèi)的類簇個數(shù),得到聚類結果在各個層次的類簇數(shù)量.

從圖4～圖6 中可以看出:相關程度相對較低的層次,即GO 富集功能目錄個數(shù)小于5 的范圍,在最優(yōu)類簇聚類時,T 細胞數(shù)據(jù)集、Th17 細胞數(shù)據(jù)集和PSC 數(shù)據(jù)集中,PUseqClust 方法中類簇數(shù)量為(4,2,4),Mclust 和MBclust中類簇數(shù)量分別為(6,5,8)和(6,2,8);在非最優(yōu)類簇聚類時,T 細胞數(shù)據(jù)集和Th17 細胞數(shù)據(jù)集中,PUseqClust 方法中類簇數(shù)量為(2,0,4),Mclust 和MBclust 中類簇數(shù)量分別為(2,1,8)和(4,1,7).由此可知,無論在最優(yōu)類簇聚類還是非最優(yōu)類簇聚類時,PUseqClust 方法都擁有最少的相關程度相對較低的類簇個數(shù).

從圖4～圖6 同樣可以看出:在高層次水平上,即GO 富集功能目錄個數(shù)大于50 的范圍內(nèi),PUseqClust 方法也能獲得不少于Mclust 和MBclust 的類簇個數(shù).所以我們認為:PUseqClust 方法與Mclust 和MBclust 相比較,聚類的性能更具競爭優(yōu)勢.

圖7 分別顯示了3 個真實數(shù)據(jù)集各個方法所有類簇GO 富集功能目錄個數(shù)的總和,用以檢測不同聚類算法獲得的聚類結果的生物學意義.從圖中可以明顯看出:PUseqClust 方法獲得了比Mclust 和MBclust 顯著多的GO富集功能目錄,因而獲得了最具生物學意義的聚類結果.

Fig.4 The cluster numbers under the specified catalogue numbers of GO enrichment function on the T cell data圖4 T 細胞數(shù)據(jù)集各種方法在指定GO 富集功能目錄數(shù)量范圍內(nèi)的類簇個數(shù)

Fig.5 The cluster numbers under the specified catalogue numbers of GO enrichment function on the Th17 cell data圖5 Th17 細胞數(shù)據(jù)集各種方法在指定GO 富集功能目錄數(shù)量范圍內(nèi)的類簇個數(shù)

Fig.6 The cluster numbers under the specified catalogue numbers of GO enrichment function on the PSC data圖6 PSC 數(shù)據(jù)集各種方法在指定GO 富集功能目錄數(shù)量范圍內(nèi)的類簇個數(shù)

Fig.7 The total catalogue numbers of GO enrichment function圖7 GO 富集功能目錄個數(shù)總和

4 總結

本文針對RNA-Seq 讀段數(shù)據(jù),在PGSeq 模型的基礎上,提出了PUseqClust 聚類框架.該聚類方法利用PGSeq降低讀段非均勻分布的影響,利用拉普拉斯反復考慮基因表達水平在不同條件下的不同重復樣本之間的相關性得到基因的對數(shù)表達水平及其不確定性,增加了模型的健壯性,并將不確定性傳遞到混合t 分布聚類模型進行聚類分析,加強了對噪聲的魯棒性.

本文采用模擬數(shù)據(jù)集和真實數(shù)據(jù)集驗證所提方法的性能,并與同樣基于模型的聚類方法Mclust 和MBclust方法進行了對比.實驗結果表明:本文方法在模擬數(shù)據(jù)集上顯示了更為準確的聚類性能,在真實數(shù)據(jù)集上獲得了更具生物意義的聚類結果.

本文的工作表明,由于RNA-seq 讀段數(shù)據(jù)本身具有多種技術性噪聲,而且RNA-Seq 實驗大都采用了生物性重復實驗,導致數(shù)據(jù)中存在一定的生物噪聲,這些對后續(xù)的數(shù)據(jù)處理形成了嚴峻挑戰(zhàn);同時,聚類分析本身也要考慮聚類過程中聚類中心的不確定度.因此,僅采用一個模型往往不能很好地對多種噪聲和不確定性進行模擬.我們通過采用多步分析,充分考慮各個步驟中的各種噪聲和不確定性,獲得了較好的分析結果.