張 旭, 吳佩望, 喬 峰

(西南大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，重慶 400715)

1 引言

到目前為止，結(jié)核病仍然是全世界人類發(fā)病和因病死亡的主要原因，全世界每年有超過(guò)1000 萬(wàn)人感染結(jié)核病，160 余萬(wàn)人因結(jié)核病而死亡，更是有四分之一的人口是潛伏性結(jié)核病患者.目前中國(guó)乃至全世界在結(jié)核病防治問(wèn)題上依然存在著很大的困難，如結(jié)核病發(fā)現(xiàn)率低、病原篩查率低、確診率低、致病細(xì)菌耐藥率高等.

目前我國(guó)對(duì)結(jié)核病的診斷主要通過(guò)痰涂片和病原分離培養(yǎng).而痰涂片采用的方法靈敏度和特異度都較低，且受實(shí)驗(yàn)室條件及實(shí)驗(yàn)人員的影響較大.病原體分離培養(yǎng)則耗時(shí)較長(zhǎng)，快速生長(zhǎng)菌也至少需要3 天才有結(jié)果，時(shí)間長(zhǎng)者甚至達(dá)到數(shù)周[1].因此，利用數(shù)據(jù)分析方法探索結(jié)核病發(fā)病機(jī)制，為疾病防控提供理論依據(jù)，并促進(jìn)了結(jié)核病感染的分子診斷及治療技術(shù)的發(fā)展.

人類全血以及外周血單核細(xì)胞的基因組表達(dá)分析已被廣泛應(yīng)用于檢測(cè)活動(dòng)性肺結(jié)核患者的宿主轉(zhuǎn)錄反應(yīng)，并識(shí)別用于診斷的潛在生物標(biāo)記物.例如：針對(duì)本文所用到的數(shù)據(jù)集，Cai 等[2]應(yīng)用全基因組轉(zhuǎn)錄微陣列分析方法測(cè)定了結(jié)核病患者和對(duì)照組外周血單核細(xì)胞基因的表達(dá)量，并通過(guò)qRT-PCR 方法驗(yàn)證C1q 表達(dá)，發(fā)現(xiàn)與健康對(duì)照組和潛伏性結(jié)核感染患者相比，活動(dòng)性結(jié)核病患者外周血中C1q 表達(dá)顯著增加，證明了C1q 表達(dá)與人類結(jié)核中的活動(dòng)性疾病相關(guān)，可能是區(qū)分活動(dòng)性與潛伏性結(jié)核病以及結(jié)核性胸膜炎與非結(jié)核性胸膜炎的潛在診斷標(biāo)志物.Blankley 等[3]對(duì)活動(dòng)性結(jié)核病與健康對(duì)照組之間進(jìn)行Benjamini Hochberg 多重檢驗(yàn)校正的獨(dú)立t檢驗(yàn)鑒別了380 個(gè)差異表達(dá)的基因.Alam 等[4]利用Limma 包在健康組和活動(dòng)性結(jié)核病患者之間鑒定出了266 個(gè)差異表達(dá)基因，其中149 個(gè)上調(diào)，117 個(gè)下調(diào)；在活動(dòng)性結(jié)核病和潛在感染組之間共發(fā)現(xiàn)127 個(gè)上調(diào)基因和69 個(gè)下調(diào)基因.

然而，生命現(xiàn)象的復(fù)雜性與不確定性使得傳統(tǒng)統(tǒng)計(jì)方法的使用存在很大的局限性，經(jīng)典統(tǒng)計(jì)學(xué)方法一系列嚴(yán)格的前提假設(shè)并不能被很好的滿足.比如，t檢驗(yàn)往往需要在數(shù)據(jù)服從正態(tài)分布的假設(shè)下使用，但基因數(shù)據(jù)多數(shù)都不滿足該假設(shè)，盲目使用將會(huì)使得研究結(jié)果與實(shí)際情況相去甚遠(yuǎn).對(duì)于數(shù)據(jù)分析結(jié)果的驗(yàn)證目前也是一個(gè)棘手的問(wèn)題.前面的這幾篇文獻(xiàn)對(duì)于選出的基因有的并沒有做進(jìn)一步驗(yàn)證，有的采用了qRT-PCR 方法驗(yàn)證，而該方法依賴于實(shí)驗(yàn)，耗時(shí)長(zhǎng)、費(fèi)用高，一般只能對(duì)極少數(shù)基因進(jìn)行驗(yàn)證.因此，為了更加準(zhǔn)確地模擬及預(yù)測(cè)具體的生物學(xué)過(guò)程，迫切需要一套基于數(shù)學(xué)、統(tǒng)計(jì)學(xué)以及計(jì)算機(jī)科學(xué)等的完整有效的組合方法，從海量的數(shù)據(jù)資源入手挖掘出隱含的、有價(jià)值的信息，與微生物學(xué)領(lǐng)域的研究方法形成互補(bǔ)，這也正是本研究的主旨.

本文先通過(guò)兩種基于貝葉斯統(tǒng)計(jì)框架的方法，即線性模型及經(jīng)驗(yàn)貝葉斯方法和信息先驗(yàn)性貝葉斯檢驗(yàn)方法相結(jié)合，篩選出了結(jié)核病的潛在易感基因.再經(jīng)過(guò)樸素貝葉斯分類器驗(yàn)證了這些易感基因的準(zhǔn)確性，突出了貝葉斯方法在基因數(shù)據(jù)分析中的重要性.然后對(duì)這些基因進(jìn)行了生物信息學(xué)分析，從生物學(xué)角度分析了結(jié)核病發(fā)病的分子機(jī)制.以期為結(jié)核病的診斷、長(zhǎng)期防控工作提供參考依據(jù).

2 數(shù)據(jù)和方法

2.1 數(shù)據(jù)來(lái)源及處理

本文所使用的數(shù)據(jù)集是來(lái)自美國(guó)國(guó)立生物技術(shù)信息中心(National Center for Biotechnology Information, NCBI)的高通量芯片表達(dá)譜數(shù)據(jù)庫(kù)中的GSE54992，其中包含了9 例結(jié)核病患者樣本(TB)，12 例結(jié)核病潛伏患者樣本和6 例正常樣本(HD)，每個(gè)樣本有22283 個(gè)探針.該數(shù)據(jù)集應(yīng)用全基因組轉(zhuǎn)錄微陣列分析方法測(cè)定了結(jié)核病患者和對(duì)照組外周血單核細(xì)胞基因的表達(dá)量[2].

本文還選取了NCBI 的另一組數(shù)據(jù)，同樣來(lái)自NCBI 的高通量芯片表達(dá)譜數(shù)據(jù)庫(kù)中的GSE83456 用作差異表達(dá)基因的驗(yàn)證，其中包含了45 例結(jié)核病患者樣本(PTB)與61 例正常樣本(HC)，每個(gè)樣本有47231 個(gè)探針.該數(shù)據(jù)集應(yīng)用全基因組轉(zhuǎn)錄微陣列分析方法測(cè)定了結(jié)核病患者和對(duì)照組外周血單核細(xì)胞基因的表達(dá)量[3].

這兩個(gè)數(shù)據(jù)集均已經(jīng)過(guò)標(biāo)準(zhǔn)化，可以直接用于分析.在將探針轉(zhuǎn)化為基因符號(hào)時(shí)，我們將對(duì)應(yīng)同一基因符號(hào)的多個(gè)探針中的表達(dá)量的最大值作為該基因的表達(dá)量.

2.2 研究方法

2.2.1 線性模型及經(jīng)驗(yàn)貝葉斯方法

線性模型及經(jīng)驗(yàn)貝葉斯方法(Limma)是由Smyth 根據(jù)文獻(xiàn)[5]提出的層次模型發(fā)展而來(lái)的，針對(duì)基因組微陣列數(shù)據(jù)“大p小n”問(wèn)題的解決方法[6].Limma 整合了多種統(tǒng)計(jì)原理，它在基因表達(dá)值矩陣上運(yùn)行，矩陣中每一行代表一個(gè)基因或者探針，每一列對(duì)應(yīng)一個(gè)樣本.一方面，通過(guò)加權(quán)或廣義最小二乘法擬合多個(gè)線性模型，并以各種方式利用這些模型的靈活性.另一方面，由于基因組數(shù)據(jù)的高度平行性，它在每一個(gè)基因模型之間借用信息，允許基因之間及樣本之間存在不同程度的差異，這使得統(tǒng)計(jì)結(jié)論在樣本量較小的數(shù)據(jù)集中更為可靠.Limma 證明可以使用經(jīng)驗(yàn)Bayes 后驗(yàn)方差估計(jì)進(jìn)行精確的小樣本推斷.這種方法在小樣本的實(shí)驗(yàn)中被證明是特別有優(yōu)勢(shì)的，確保了即使在重復(fù)次數(shù)很少的情況下，也能得到可靠和穩(wěn)定的統(tǒng)計(jì)推斷.Limma 是選取差異基因的一種基礎(chǔ)的，傳統(tǒng)的方法，該方法雖然緩解了由于樣本量小而導(dǎo)致的推理結(jié)果不佳的問(wèn)題，但是在估計(jì)低方差或高方差的基因時(shí)會(huì)引入偏差.

2.2.2 信息先驗(yàn)性貝葉斯檢驗(yàn)

信息先驗(yàn)性貝葉斯檢驗(yàn)(IPBT)是由Li 等[7]于2015 提出的一種替代方法，它在某種意義上與基于貝葉斯層次模型的Limma 方法“垂直”.Limma 方法中使用了經(jīng)驗(yàn)貝葉斯方法借用同一次實(shí)驗(yàn)中不同基因的信息，而IPBT 方法則是借用在過(guò)去不同實(shí)驗(yàn)中相同的基因(或探針)，但是使用相同的技術(shù)，相同類型的芯片，在同一類型的細(xì)胞上的測(cè)量信息.該方法的關(guān)鍵思想是基于豐富的歷史數(shù)據(jù)為每個(gè)基因(或探針)指定了特異的信息性先驗(yàn)分布然后再進(jìn)行貝葉斯假設(shè)檢驗(yàn).因?yàn)椴煌幕蚓哂胁煌纳飳W(xué)功能，所以通常情況下它們的表達(dá)量會(huì)顯示出相當(dāng)多樣化的分布情況.IPBT 的模型為

因此，與Limma 在同一實(shí)驗(yàn)中借用不同基因的策略不同，IPBT 假設(shè)每個(gè)基因都有屬于自己獨(dú)特的先驗(yàn)分布，然后在一次實(shí)驗(yàn)中借用這些歷史數(shù)據(jù)中的先驗(yàn)信息來(lái)分析.經(jīng)過(guò)與仿真數(shù)據(jù)的分析對(duì)比并繪制ROC 曲線，這種方法被證明篩選結(jié)果優(yōu)于SAM、t檢驗(yàn)及Limma 等其他方法[8].

2.2.3 樸素貝葉斯分類器

樸素貝葉斯是一種構(gòu)造分類器的簡(jiǎn)單技術(shù)，是在貝葉斯算法的基礎(chǔ)上進(jìn)行相應(yīng)的簡(jiǎn)化，即假設(shè)給定目標(biāo)值時(shí)屬性之間相互條件獨(dú)立.也就是說(shuō)沒有哪個(gè)屬性變量對(duì)于決策結(jié)果來(lái)說(shuō)占有著較大或者較小的比重.雖然這個(gè)簡(jiǎn)化方式在一定程度上對(duì)貝葉斯分類算法的分類效果有稍許影響，但是在實(shí)際的應(yīng)用場(chǎng)景中，極大地簡(jiǎn)化了貝葉斯方法的復(fù)雜性[8].樸素貝葉斯主要基于條件概率模型：給定一個(gè)要分類的問(wèn)題，用向量x= (x1,x2,1,··· ,xn)代表n個(gè)特征，對(duì)于K個(gè)可能的類別中的每一類Ck，使用貝葉斯定理，將每一類的概率寫為

其中Ck類變量的條件分布是

Z是一個(gè)只依賴于特征x的縮放因子，當(dāng)特征變量的值已知時(shí)是一個(gè)常數(shù).然后根據(jù)建立的樸素貝葉斯概率模型和最大后驗(yàn)概率決策準(zhǔn)則構(gòu)造樸素貝葉斯分類器[9]

樸素貝葉斯分類器具有穩(wěn)定的分類效率，在小規(guī)模數(shù)據(jù)集上的表現(xiàn)很好，能夠處理多分類任務(wù)，并且對(duì)于缺失數(shù)據(jù)不敏感，在理論上樸素貝葉斯分類器與其他分類方法相比具有最小的誤差率.

2.2.4 生物信息學(xué)分析

Gene Ontology(簡(jiǎn)稱GO)是一個(gè)國(guó)際標(biāo)準(zhǔn)化的基因功能分類體系，是一項(xiàng)重要的生物信息學(xué)計(jì)劃.而GO 富集分析是一種利用基因本體論分類系統(tǒng)解釋基因組的技術(shù)，其中基因根據(jù)其功能特征被分配給一組預(yù)定義的領(lǐng)域，如細(xì)胞組分、分子功能、生物過(guò)程[10].本文中，我們利用文獻(xiàn)[11]提供的R 語(yǔ)言包c(diǎn)lusterProfiler 中的enrichGO 函數(shù)進(jìn)行了GO 功能富集分析，并將其結(jié)果可視化.對(duì)差異表達(dá)基因進(jìn)行生物信息學(xué)分析的另一種方案是KEGG 數(shù)據(jù)庫(kù)[12]，其中我們將根據(jù)基因所涉及的通路進(jìn)行分組，并通過(guò)假設(shè)檢驗(yàn)分析這些差異表達(dá)基因是否在某些通路上存在富集.本文中我們利用文獻(xiàn)[11]提供的R 語(yǔ)言包c(diǎn)lusterProfiler 中的enrichKEGG 函數(shù)進(jìn)行了KEGG 通路分析.

3 結(jié)果與分析

3.1 篩選結(jié)果

GSE54992 數(shù)據(jù)集經(jīng)過(guò)合并探針后，我們共得到13516 個(gè)基因的表達(dá)量數(shù)據(jù).本文僅選取了其中的結(jié)核病患者和正常人兩種樣本數(shù)據(jù)進(jìn)行對(duì)比分析.經(jīng)過(guò)繪制PCA 圖，如圖1 所示，可以看出兩組樣本被明顯分開，使我們有可能通過(guò)數(shù)據(jù)對(duì)比尋找差異基因.

圖1 GES54992 兩組樣本的PCA 圖，HD 為正常樣本，TB 為患病樣本

我們先使用R 語(yǔ)言中的Limma 包處理基因表達(dá)量數(shù)據(jù)，一共得到558 個(gè)差異表達(dá)基因，其中279 個(gè)為上調(diào)表達(dá)基因，另外279 個(gè)為下調(diào)表達(dá)基因.接著我們又采用了信息先驗(yàn)性貝葉斯方法再一次篩選差異基因.基于文獻(xiàn)[7]提供的GitHub 中下載得到IPBT 包，我們一共得到821 個(gè)差異表達(dá)基因，其中上調(diào)基因293 個(gè)，下調(diào)基因528 個(gè).

為了使篩選的差異表達(dá)基因結(jié)果更為準(zhǔn)確，我們對(duì)兩種方法篩選出的差異表達(dá)基因取交集，得到了319 個(gè)基因，包括114 個(gè)上調(diào)基因和205 個(gè)下調(diào)基因.為了能更直觀地看出結(jié)果，我們進(jìn)行了可視化，繪制了這319 個(gè)基因的韋恩圖和熱圖，如圖2 和圖3 所示.

圖2 IPBT 與Limma 篩選出的差異表達(dá)基因數(shù)量韋恩圖，其中319 個(gè)基因被兩種方法均識(shí)別

圖3 319 個(gè)基因及15 個(gè)樣本的表達(dá)熱圖

圖3 中橫軸是根據(jù)患病情況分組的15 個(gè)樣本，縱軸是選出的319 個(gè)基因，每個(gè)小方塊為一個(gè)基因的表達(dá)量.顏色的深淺代表該基因表達(dá)量的高低，紅色越深表明基因表達(dá)量越高，藍(lán)色越深表示基因表達(dá)量越低.從圖中我們可以明顯看出樣本被分為兩個(gè)組別，而上調(diào)與下調(diào)基因之間表達(dá)量也存在著很大的差異.表明我們通過(guò)這兩種方法篩選出的差異表達(dá)基因的效果較好.

3.2 模型驗(yàn)證

由于數(shù)據(jù)集GSE54992 的樣本本身的區(qū)分度較好(圖1)，我們選取了另一組數(shù)據(jù)GSE83456 來(lái)驗(yàn)證差異基因的分類效果.從圖4 可以看出該數(shù)據(jù)集的兩種樣本有很大部分相互交疊在一起，難以區(qū)分開.

圖4 GSE83456 兩組樣本的PCA 圖，其中HC 為正常樣本，PTB 為患病樣本

對(duì)于GSE83456 這個(gè)獨(dú)立驗(yàn)證集，我們采用兩種方式選取了特征基因建立樸素貝葉斯分類模型進(jìn)行對(duì)比，其一為IPBT 與Limma 兩種方法共同篩選得到的319 個(gè)差異表達(dá)基因，其二為隨機(jī)選取319 個(gè)基因.在70%訓(xùn)練集，30%驗(yàn)證集的條件下建立樸素貝葉斯分類模型，通過(guò)ROC 曲線我們發(fā)現(xiàn)，隨機(jī)選取特征基因用作建模的準(zhǔn)確率僅為55.6%，而通過(guò)基于貝葉斯方法篩選出來(lái)的差異表達(dá)基因用作特征基因?qū)υ摂?shù)據(jù)集進(jìn)行建模得到的分類準(zhǔn)確率超過(guò)了85%，遠(yuǎn)遠(yuǎn)高于隨機(jī)選取基因用作特征基因的表現(xiàn).因此，我們認(rèn)為篩選出的差異表達(dá)基因在結(jié)核病患病與否的預(yù)測(cè)上有較為顯著的作用.

3.3 生物信息學(xué)分析結(jié)果

我們先將差異表達(dá)基因分為上調(diào)基因與下調(diào)基因，再分別進(jìn)行了GO 功能富集分析，并根據(jù)p值大小分別選出上調(diào)和下調(diào)的各前十名，共20 項(xiàng)功能.經(jīng)過(guò)文獻(xiàn)查閱，我們發(fā)現(xiàn)選出的功能類別中有10 項(xiàng)功能與結(jié)核病發(fā)病相關(guān)，其中包括細(xì)胞遷移的正調(diào)節(jié)[13]、白細(xì)胞遷移[13]、髓樣白細(xì)胞遷移[13]、細(xì)胞趨化性[14]、白細(xì)胞趨化性[14]、血管系統(tǒng)發(fā)育的調(diào)節(jié)[15]、血管生成[15]、運(yùn)動(dòng)的正調(diào)節(jié)[16]、細(xì)胞成分運(yùn)動(dòng)的正調(diào)節(jié)以及細(xì)胞運(yùn)動(dòng)的正調(diào)節(jié)[16].表明我們所篩選出來(lái)的差異表達(dá)基因中大多數(shù)基因被富集在與結(jié)核病發(fā)病相關(guān)的基因功能中.隨后我們對(duì)選出的319 個(gè)差異表達(dá)基因進(jìn)行KEGG 通路分析.由于無(wú)論是上調(diào)或是下調(diào)的基因均可能產(chǎn)生同一生物效應(yīng)，或在同一通路中出現(xiàn)，因此，我們使用所有的差異表達(dá)基因一起進(jìn)行KEGG 通路富集分析，通過(guò)查閱文獻(xiàn)發(fā)現(xiàn)，p值列于前十的通路中，細(xì)胞因子-受體相互作用[17]、趨化因子信號(hào)通路[18]、抗生素的生物合成[19]與結(jié)核病發(fā)病相關(guān)，結(jié)核病與百日咳[20]、軍團(tuán)病[21]、阿米巴病[22]同屬呼吸道疾病并且存在并發(fā)癥的情況，更有文章表明用腫瘤壞死因子抑制劑治療類風(fēng)濕性關(guān)節(jié)炎可能會(huì)導(dǎo)致患活動(dòng)性結(jié)核病的風(fēng)險(xiǎn)顯著增加[23].此外，在所有p值小于0.05 的通路中，我們發(fā)現(xiàn)有9 個(gè)基因被富集在編號(hào)為hsa05152 的名稱為“結(jié)核病”的通路中，這9 個(gè)基因分別為：IL1A、CAMK2A、IL1B、SPHK1、NFKB1、IL6、CR1、IL10 和JAK2.

4 結(jié)論與討論

本文通過(guò)Limma 和IPBT 兩種基于貝葉斯統(tǒng)計(jì)的方法從數(shù)據(jù)集GSE54992 中篩選出319 個(gè)差異表達(dá)基因；再基于獨(dú)立驗(yàn)證集GSE83456 對(duì)篩選出的319 個(gè)差異表達(dá)基因利用樸素貝葉斯分類器進(jìn)行驗(yàn)證，使原本難以分開的兩組樣本得到了較準(zhǔn)確地分類，證明了篩選的差異表達(dá)基因的準(zhǔn)確性，同時(shí)突出了貝葉斯方法在基因數(shù)據(jù)分析中的重要性；最后對(duì)這些差異表達(dá)基因進(jìn)行GO 功能富集分析和KEGG 通路分析，找到了多個(gè)與結(jié)核病相關(guān)的功能與通路，使篩選出的差異基因得到了進(jìn)一步驗(yàn)證.該研究對(duì)結(jié)核病生物標(biāo)志物的篩選以及結(jié)核病的診斷與防控提供了新的思路.

本文巧妙地組合了三種貝葉斯方法，在傳統(tǒng)的經(jīng)驗(yàn)貝葉斯基礎(chǔ)上結(jié)合了信息先驗(yàn)性貝葉斯檢驗(yàn)來(lái)篩選差異基因，使得結(jié)果更加可靠；對(duì)選出的差異基因又利用樸素貝葉斯分類器進(jìn)行驗(yàn)證，通過(guò)基因指導(dǎo)結(jié)核病樣本分類的高準(zhǔn)確率說(shuō)明了選出基因的重要性.每個(gè)步驟都是基于貝葉斯框架，有利于減少系統(tǒng)誤差.此外，IPBT 方法還有很大的提升空間，例如，目前IPBT 提供的歷史數(shù)據(jù)僅有來(lái)自GPL96 平臺(tái)的數(shù)十種細(xì)胞的微陣列數(shù)據(jù)，今后可以更多地收集相關(guān)歷史數(shù)據(jù)以獲得更準(zhǔn)確的差異表達(dá)分析結(jié)果.