基于差異共表達(dá)分析的肝癌特異性基因的篩選與驗(yàn)證

岳宇巍， 王化琨

（黑龍江大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，哈爾濱150080）

0 引 言

肝癌是我國最常見的惡性腫瘤之一，主要包括兩種病理組織學(xué)類型：肝細(xì)胞癌（Hepatocellular carcinoma，HCC）和肝內(nèi)膽管細(xì)胞癌（Intrahepatic cholangiocarcinoma，ICCA），其中HCC占我國肝癌總數(shù)的83.9%～92.3%［1］。據(jù)統(tǒng)計(jì)，全球每年約有70萬例新發(fā)肝癌患者，其中大約有35萬例肝癌患者在中國［2-3］。

隨著DNA測序技術(shù)的成熟，基因表達(dá)譜數(shù)據(jù)已廣泛應(yīng)用于癌癥研究，如應(yīng)用基因的差異表達(dá)分析（differential expression analysis，DEA）方法和生物信息學(xué)工具比較腫瘤和正常組織的基因平均表達(dá)水平的差異，挖掘癌癥相關(guān)的分子標(biāo)志物［4-5］。DEA是通過識(shí)別基因的高表達(dá)或低表達(dá)篩選潛在的腫瘤標(biāo)志物，但并沒有充分利用微陣列數(shù)據(jù)，因?yàn)樗皇褂昧藖碜赃x定基因的信息，而未使用來自整個(gè)轉(zhuǎn)錄組的信息，且沒有考慮它們之間的相互作用［6］。差異共表達(dá)分析（Differential coexpression analysis，DCA）可以作為DEA的補(bǔ)充，通過比較共表達(dá)網(wǎng)絡(luò)，認(rèn)為具有強(qiáng)烈改變的連接性的基因在疾病表型中起重要作用。隨著公開轉(zhuǎn)錄組學(xué)研究的快速積累，結(jié)合多個(gè)轉(zhuǎn)錄組學(xué)研究的共表達(dá)分析，可以提供更準(zhǔn)確和穩(wěn)健的結(jié)果［7］。本文在Marjan等研究5種人類組織的差異共表達(dá)和平均表達(dá)水平的混雜效應(yīng)的基礎(chǔ)上，結(jié)合GEO數(shù)據(jù)庫篩選25個(gè)肝組織數(shù)據(jù)集，根據(jù)肝癌發(fā)生、發(fā)展的3個(gè)階段，構(gòu)建了新的健康、肝炎和肝硬化的特異性基因?qū)脖磉_(dá)分?jǐn)?shù)集合，并計(jì)算了肝癌共表達(dá)基因的特異性分?jǐn)?shù)［8］。經(jīng)驗(yàn)證得到了較好的結(jié)果，篩選出肝癌特異性共表達(dá)基因，用STRING數(shù)據(jù)庫對這些特異性共表達(dá)基因構(gòu)建蛋白-蛋白相互作用（Protein-protein interaction，PPI）網(wǎng)絡(luò)，應(yīng)用Cytoscape軟件得到了Hub基因和基因模塊，同時(shí)，利用GEPIA在線分析工具得到關(guān)鍵基因的差異表達(dá)信息及患者生存曲線。利用DAVID在線分析進(jìn)行GO和KEGG功能富集分析，篩選出與肝癌發(fā)生發(fā)展相關(guān)的關(guān)鍵基因和通路，從基因相互作用角度為肝癌的發(fā)病分子機(jī)制提供補(bǔ)充和依據(jù)。

1 方法

1.1 數(shù)據(jù)描述

從GEO數(shù)據(jù)庫（http：//www.ncbi.nlm.nih.gov/geo/）下載Affymertix人基因組U133 2.0芯片［HG-U133＿Plus＿2］同一平臺(tái)號(hào)（安捷倫GPL570平臺(tái)）的25個(gè)肝組織基因表達(dá)譜數(shù)據(jù)。首先應(yīng)用R語言affy包中的函數(shù)對數(shù)據(jù)進(jìn)行預(yù)處理，包括背景校正（Rma）、標(biāo)準(zhǔn)化（Quantiles）、PM校正（Pmonly）和匯總（Medianpolish），然后使用Gemma異常值檢測算法去除異常樣本［9］，再根據(jù)樣本信息用R語言sva包的ComBat函數(shù)移除批次效應(yīng)，接下來過濾掉平均表達(dá)值低的探針，并根據(jù)樣本信息分成4類（健康、肝炎、肝硬化和肝癌）樣本，得到33個(gè)基因表達(dá)數(shù)據(jù)矩陣，數(shù)據(jù)集及樣本分類信息見表1。設(shè)置每個(gè)類別的基因平均表達(dá)值在前80%的基因被認(rèn)為表達(dá)。為了進(jìn)一步過濾肝癌數(shù)據(jù)集的基因，結(jié)合了用R語言處理后的TCGA數(shù)據(jù)庫中腫瘤純度大于60%的340個(gè)肝癌樣本，選擇平均Counts值在前75%的基因，與GEO得到的肝癌基因取交集得到肝癌基因的研究范圍。

表1 GEO數(shù)據(jù)信息Table 1 GEO data information

1.2 差異共表達(dá)分析模型

本文的目的是找到在肝癌組織中的高共表達(dá)，而在非肝癌三個(gè)組織中低共表達(dá)的基因共表達(dá)鏈接。首先在每個(gè)數(shù)據(jù)集中，使用皮爾遜相關(guān)系數(shù)（Pearson correlation coefficient）為每個(gè)數(shù)據(jù)集建立共表達(dá)矩陣。共表達(dá)值定義為：Sij=｛corr（i，j，k） i≠j，1≤k≤33｝，其中數(shù)據(jù)集D（k）中基因i和基因j的皮爾遜相關(guān)值表示為corr（i，j，k）（這里沒有考慮負(fù)相關(guān)，把負(fù)相關(guān)定義為0）。對于肝癌的15個(gè)數(shù)據(jù)集中本文選取基因?qū)脖磉_(dá)值在前10%的基因?qū)?，并確定在n個(gè)數(shù)據(jù)集中都存在才被認(rèn)為在肝癌中高共表達(dá)并納入研究范圍，用二項(xiàng)分布作為零假設(shè)，控制錯(cuò)誤發(fā)現(xiàn)率（False discovery rate，F(xiàn)DR）為10-4，為肝癌共表達(dá)網(wǎng)絡(luò)選擇合理的密度［8］。為肝癌共表達(dá)網(wǎng)絡(luò)的每對鏈接都計(jì)算肝癌特異性分?jǐn)?shù)（Liver cancer specificity score，LCSS），首先定義兩個(gè)集合（1）和（2），集合（1）定義為在肝癌組織中基因i和基因j的共表達(dá)值，集合（2）定義為基因i和基因j分別在健康、肝炎、肝硬化數(shù)據(jù)集上的平均共表達(dá)值，即：

最后特異性分?jǐn)?shù)［8］定義為（3），即：

應(yīng)用Wilcoxon秩和檢驗(yàn)的p值比較了肝癌和其他肝組織中基因?qū)脖磉_(dá)值的秩，Wilcoxon秩和檢驗(yàn)可以檢驗(yàn)基因?qū)υ谀[瘤組織和其他三個(gè)組織的各個(gè)數(shù)據(jù)集上的共表達(dá)值是否有顯著差異，結(jié)果發(fā)現(xiàn)這兩種方法結(jié)果高度相關(guān)。由于p值越小差異越顯著，而特異性分?jǐn)?shù)越大差異越顯著，于是將p值做負(fù)對數(shù)變換，發(fā)現(xiàn)LCSS與-log10（p-值）相關(guān)性為0.88，肝癌特異性分?jǐn)?shù)可以用來表示基因?qū)υ诟伟┙M織上的共表達(dá)特異性。應(yīng)用控制錯(cuò)誤發(fā)現(xiàn)率的方法為LCSS選擇合理的閾值［8］，為15個(gè)肝癌共表達(dá)網(wǎng)絡(luò)創(chuàng)建了30個(gè)隨機(jī)子集（Random-dataset）作為LCSS的零分布，并且同樣計(jì)算隨機(jī)子集中每個(gè)鏈接的LCSS，LCSS-FDR定義為：

其中γlcssrd是隨機(jī)子集中大于閾值LCSS的個(gè)數(shù)，為肝癌組織選擇LCSS的閾值控制錯(cuò)誤發(fā)現(xiàn)率為0.01。

1.3 PPI網(wǎng)絡(luò)構(gòu)建和PPI網(wǎng)絡(luò)挖掘

STRING［10］（https：//string-db.org/cgi/input.pl/）是已知預(yù)測的蛋白質(zhì)-蛋白質(zhì)相互作用的數(shù)據(jù)庫，利用STRING數(shù)據(jù)庫構(gòu)建蛋白質(zhì)相互作用的PPI網(wǎng)絡(luò)，再應(yīng)用Cytoscape軟件進(jìn)行網(wǎng)絡(luò)的可視化，CytoHubba插件以節(jié)點(diǎn)度為篩選條件獲得Hub基因，使用MCODE插件獲得了重要的基因模塊。

1.4 關(guān)鍵基因差異表達(dá)驗(yàn)證及生存分析

應(yīng)用GEPIA［11］數(shù)據(jù)庫對20個(gè)Hub基因進(jìn)行差異表達(dá)分析及在線生存分析，差異表達(dá)分析驗(yàn)證條件為。生存分析篩選條件為LIHC數(shù)據(jù)集，置信區(qū)間為95%，Hub基因表達(dá)量與預(yù)后的關(guān)系采用Log-rank檢驗(yàn)，有統(tǒng)計(jì)學(xué)意義的差異表示為Log-rank p＜0.01或p＜0.05。應(yīng)用The human protein atlas［12］（https：//www.proteinatlas.org/）得到Hub基因的肝癌預(yù)后總結(jié)。

1.5 基因模塊的富集分析

利用DAVID［13］（https：//david.ncifcrf.gov/）在線富集分析對模塊內(nèi)的基因進(jìn)行分子生物學(xué)功能（Molecular function，MF）、生物學(xué)過程（Biological process，BP）、細(xì)胞學(xué)組分（Cellular component，CC）的GO功能富集分析，KEGG通路分析，納入標(biāo)準(zhǔn)為p＜0.05。

2 結(jié) 果

2.1 利用差異共表達(dá)分析挖掘肝癌特異性共表達(dá)基因

對預(yù)處理后的肝癌數(shù)據(jù)篩選，最終確定以15個(gè)數(shù)據(jù)集（559個(gè)樣本和8 759個(gè)基因）為肝癌組織研究對象，基因?qū)χ辽僭?2個(gè)數(shù)據(jù)集中都存在才被認(rèn)為是在肝癌組織中高共表達(dá)，從而控制肝癌共表達(dá)網(wǎng)絡(luò)密度為0.007，得到了196 589個(gè)基因?qū)?。LCSS-FDR控制分?jǐn)?shù)的閾值為0.49，閾值過濾后的肝癌特異性網(wǎng)絡(luò)（Liver cancer specific network，LCSN）包含3 698個(gè)基因節(jié)點(diǎn)和12 515條邊。在LCSN中，選擇大于網(wǎng)絡(luò)平均連通度6的976個(gè)肝癌特異性基因作為構(gòu)建PPI網(wǎng)絡(luò)的對象。

2.2 肝癌PPI網(wǎng)絡(luò)的構(gòu)建和分析

有研究表明，PPI基因?qū)Φ谋磉_(dá)相關(guān)性比非PPI基因?qū)Φ谋磉_(dá)相關(guān)性更高，在部分PPI基因?qū)χ杏^察到異常高的差異共表達(dá)值，并且與高差異表達(dá)的基因相比，高差異共表達(dá)基因富含更多的肝癌基因［14］。將976個(gè)肝癌特異性共表達(dá)基因輸入STRING數(shù)據(jù)庫構(gòu)建PPI網(wǎng)絡(luò)（圖1 a），并進(jìn)一步用Cytoscape軟件得到了20個(gè)Hub基因和連接緊密的基因模塊。

2.3 Hub基因和基因模塊的篩選和驗(yàn)證

將上述方法得到的20個(gè)Hub基因構(gòu)建成基因模塊（圖1b），利用DAVID網(wǎng)站進(jìn)行KEGG富集分析，結(jié)果（表2）確定了9個(gè)關(guān)鍵基因富集到重要癌癥通路，包括癌癥、p53信號(hào)、細(xì)胞周期、Wnt、PI3k-Akt信號(hào)和病毒致癌作用通路。差異表達(dá)分析結(jié)果顯示，有6個(gè)基因存在顯著的差異表達(dá)（圖3），生存曲線分析結(jié)果（圖4）顯示，CDK4、RAC1、CHEK1、SPP1、HDAC1和UBE2D1表達(dá)的升高和ESR1表達(dá)的降低會(huì)顯著降低肝癌患者的總體生存率（Log-rank P＜0.01）。對沒有參與這些重要通路的11個(gè)Hub基因，通過已發(fā)表資料研究最終識(shí)別HDAC1、APOB、UBE2D1、ELAVL1、ATG7和MSH2為可能與肝癌發(fā)生發(fā)展相關(guān)的關(guān)鍵基因。

圖1 肝癌特異性的976個(gè)高連通度基因所構(gòu)建的PPI網(wǎng)絡(luò)和樞紐基因模塊Fig.1 PPI network constructed by 976 highly connected genes with liver cancer specific and Hub gene module

表2 蛋白-蛋白相互作用網(wǎng)絡(luò)基因富集分析Table 2 Gene enrichment analysis of PPI network

圖2 肝癌特異性對應(yīng)的蛋白質(zhì)互作網(wǎng)絡(luò)中篩選的3個(gè)高度互聯(lián)的模塊Fig.2 Three modules with the high interconnection screened from the PPI networks with liver cancer specific

2.4 肝癌模塊參與的調(diào)控作用

用DAVID富集分析［15］篩選到具有高互聯(lián)和生物學(xué)意義的3個(gè)基因模塊（圖2），其中模塊1包含了120個(gè)基因，形成了210種相互作用的關(guān)系。表2列舉了模塊富集分析的主要功能，模塊1基因富集的主要功能為蛋白質(zhì)泛素化和泛素蛋白轉(zhuǎn)移酶活性，KEGG通路為泛素介導(dǎo)的蛋白質(zhì)水解。模塊2基因富集的主要功能為脂蛋白代謝過程、脂蛋白顆粒相關(guān)功能和多聚腺苷酸核糖核酸，KEGG通路為剪接體和PPAR信號(hào)通路。模塊3基因富集的主要功能為胞漿，KEGG通路為P13K-Akt和癌癥的中心碳代謝通路。

圖3 20個(gè)Hub基因的差異表達(dá)分析結(jié)果（紅色為腫瘤組，灰色為健康組）Fig.3 Differential expression analysis results of 20 Hub genes（red：tumor group，gray：healthy group）

圖4 生存分析結(jié)果及患者預(yù)后的生存曲線（紅色為腫瘤組，藍(lán)色為健康組）Fig.4 Survival analysis results and patient prognosis survival curve（red：tumor group，blue：healthy group）

3 討 論

隨著高通量測序和芯片技術(shù)的日益成熟，生成大規(guī)模、多組織的基因表達(dá)數(shù)據(jù)已經(jīng)成為現(xiàn)實(shí)。在疾病研究中，以基因表達(dá)譜為研究對象、利用生物信息學(xué)工具分析的腫瘤研究較多，在眾多基因表達(dá)數(shù)據(jù)中挖掘肝癌新型的標(biāo)志物，為肝癌的診斷與治療靶點(diǎn)選擇及預(yù)后判斷提供參考具有重要意義。

本研究以差異共表達(dá)分析方法通過肝癌發(fā)生發(fā)展過程的3類數(shù)據(jù)集（健康、肝炎和肝硬化），準(zhǔn)確識(shí)別了高差異共表達(dá)的肝癌基因?qū)?，結(jié)合生物信息學(xué)工具，對973個(gè)肝癌特異性基因進(jìn)行生物信息學(xué)分析，最終得到了3個(gè)基因模塊和20個(gè)Hub基因。在20個(gè)Hub基因中，KEGG通路富集分析結(jié)果（表2）顯示基因CCND1、CDK4、RAC1、CHEK1、RAC2、TP53、ESR1和SPP1主要參與了p53信號(hào)、細(xì)胞周期、Wnt、PI3k-Akt信號(hào)與病毒致癌作用等和肝癌發(fā)生發(fā)展相關(guān)的重要通路，這些參與重要癌癥通路的基因已經(jīng)被廣泛研究。除了這些Hub基因，已有研究資料確認(rèn)，HDAC1［16-19］、APOB［20］、UBE2D1［21］和ELAVL1［22］雖然沒有參與這些重要的癌癥通路，但是這些基因在肝癌的發(fā)生發(fā)展過程中起著重要的作用，并且這些基因并沒有顯著的差異表達(dá)（圖3），傳統(tǒng)的差異表達(dá)方法篩選不到這些基因，所以傳統(tǒng)的差異表達(dá)研究方法并沒有全部利用基因表達(dá)譜的全部信息，并且生存分析的結(jié)果顯示CDK4、RAC1、CHEK1、ESR1、SPP1、HDAC1、UBE2D1、ATG7和PRPF8的表達(dá)差異會(huì)明顯降低肝癌患者的總體生存率，這些基因在人類蛋白質(zhì)圖譜的驗(yàn)證中除基因ESR1外，均有對肝癌不利的預(yù)后。除了以上基因，MSH2和ATG7基因也可能與肝癌相關(guān)，其中MSH2有顯著差異表達(dá)（圖3），并且MSH2和ATG7的高表達(dá)也對患者的生存率有顯著影響（Log-rank p＜0.05），但是這兩個(gè)基因現(xiàn)在還沒有被納入肝癌研究的靶點(diǎn)。本文提出的差異共表達(dá)分析方法有效的識(shí)別出了肝癌的關(guān)鍵基因，可以應(yīng)用在其他類型的多數(shù)據(jù)集研究，以選擇其他復(fù)雜疾病的關(guān)鍵基因。本文結(jié)果可以作為肝癌差異表達(dá)分析研究結(jié)論的補(bǔ)充，為肝癌的診斷和治療靶點(diǎn)選擇及預(yù)后判斷提供參考，多個(gè)數(shù)據(jù)集的聯(lián)合分析使結(jié)論更加具有穩(wěn)健性。