基于差異共表達分析的肝癌特異性基因的篩選與驗證

岳宇巍， 王化琨

（黑龍江大學 數學科學學院，哈爾濱150080）

0 引 言

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

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

1 方法

1.1 數據描述

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

表1 GEO數據信息Table 1 GEO data information

1.2 差異共表達分析模型

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

最后特異性分數［8］定義為（3），即：

應用Wilcoxon秩和檢驗的p值比較了肝癌和其他肝組織中基因對共表達值的秩，Wilcoxon秩和檢驗可以檢驗基因對在腫瘤組織和其他三個組織的各個數據集上的共表達值是否有顯著差異，結果發現這兩種方法結果高度相關。由于p值越小差異越顯著，而特異性分數越大差異越顯著，于是將p值做負對數變換，發現LCSS與-log10（p-值）相關性為0.88，肝癌特異性分數可以用來表示基因對在肝癌組織上的共表達特異性。應用控制錯誤發現率的方法為LCSS選擇合理的閾值［8］，為15個肝癌共表達網絡創建了30個隨機子集（Random-dataset）作為LCSS的零分布，并且同樣計算隨機子集中每個鏈接的LCSS，LCSS-FDR定義為：

其中γlcssrd是隨機子集中大于閾值LCSS的個數，為肝癌組織選擇LCSS的閾值控制錯誤發現率為0.01。

1.3 PPI網絡構建和PPI網絡挖掘

STRING［10］（https：//string-db.org/cgi/input.pl/）是已知預測的蛋白質-蛋白質相互作用的數據庫，利用STRING數據庫構建蛋白質相互作用的PPI網絡，再應用Cytoscape軟件進行網絡的可視化，CytoHubba插件以節點度為篩選條件獲得Hub基因，使用MCODE插件獲得了重要的基因模塊。

1.4 關鍵基因差異表達驗證及生存分析

應用GEPIA［11］數據庫對20個Hub基因進行差異表達分析及在線生存分析，差異表達分析驗證條件為。生存分析篩選條件為LIHC數據集，置信區間為95%，Hub基因表達量與預后的關系采用Log-rank檢驗，有統計學意義的差異表示為Log-rank p＜0.01或p＜0.05。應用The human protein atlas［12］（https：//www.proteinatlas.org/）得到Hub基因的肝癌預后總結。

1.5 基因模塊的富集分析

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

2 結 果

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

對預處理后的肝癌數據篩選，最終確定以15個數據集（559個樣本和8 759個基因）為肝癌組織研究對象，基因對至少在12個數據集中都存在才被認為是在肝癌組織中高共表達，從而控制肝癌共表達網絡密度為0.007，得到了196 589個基因對。LCSS-FDR控制分數的閾值為0.49，閾值過濾后的肝癌特異性網絡（Liver cancer specific network，LCSN）包含3 698個基因節點和12 515條邊。在LCSN中，選擇大于網絡平均連通度6的976個肝癌特異性基因作為構建PPI網絡的對象。

2.2 肝癌PPI網絡的構建和分析

有研究表明，PPI基因對的表達相關性比非PPI基因對的表達相關性更高，在部分PPI基因對中觀察到異常高的差異共表達值，并且與高差異表達的基因相比，高差異共表達基因富含更多的肝癌基因［14］。將976個肝癌特異性共表達基因輸入STRING數據庫構建PPI網絡（圖1 a），并進一步用Cytoscape軟件得到了20個Hub基因和連接緊密的基因模塊。

2.3 Hub基因和基因模塊的篩選和驗證

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

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

表2 蛋白-蛋白相互作用網絡基因富集分析Table 2 Gene enrichment analysis of PPI network

圖2 肝癌特異性對應的蛋白質互作網絡中篩選的3個高度互聯的模塊Fig.2 Three modules with the high interconnection screened from the PPI networks with liver cancer specific

2.4 肝癌模塊參與的調控作用

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

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

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

3 討 論

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

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