基因調(diào)控網(wǎng)絡中的癌癥標記物預測方法

覃桂敏，劉佳妍，殷 雨，楊璐瓊

(西安電子科技大學計算機科學與技術學院，陜西西安710071)

復雜疾病是一類不服從孟德爾遺傳定律的疾病，其發(fā)生過程涉及多個復雜的生物學過程，并由多種遺傳物質協(xié)同控制。近年來，關于復雜疾病的癌癥標記物研究層出不窮。隨著高通量技術的發(fā)展，越來越多的生物數(shù)據(jù)，譬如基因表達數(shù)據(jù)、體細胞突變數(shù)據(jù)、蛋白質互作用數(shù)據(jù)等被證實可以應用于癌癥標記物的預測中。

在生物醫(yī)學界對癌癥發(fā)病機制的認知進入到分子層面后，逐漸有其他相關領域的研究人員對癌癥標記物進行了多種多樣的分析。文獻[1]以食管癌(ESophageal CArcinoma，ESCA)為背景檢測了50例食管癌樣本和20例正常樣本，發(fā)現(xiàn)HOTAIR分子的異常表達可用于食管癌的早期診斷。除此以外，對于曾經(jīng)被證實和一些特定癌癥的特異過程相關的分子，也進行了重新探討。SLC52A3一直被認為可以促進核黃素的吸收，調(diào)節(jié)機體的代謝過程。然而文獻[2]表明，SLC52A3在轉錄過程中因為轉錄起始位點的不同，會產(chǎn)生兩個轉錄變體SLC52A3a和SLC52A3b。其中的SLC52A3a在食管癌中表達機制異常，具體表現(xiàn)為促進食管癌細胞的增殖。

胃癌(STomach ADenocarcinoma，STAD)和食管癌這些惡性腫瘤屬于復雜疾病。胃癌和食管癌早期沒有明顯癥狀，而確診后往往已經(jīng)是中晚期癌細胞擴散后，因此，致死率極高。故關于胃癌和食管癌的病理學研究一直在進行中。文獻[3]整合基因表達綜合數(shù)據(jù)庫(Gene Expression Omnibus，GEO)中的大量胃癌基因表達譜，提取出在胃癌中差異表達的433個基因，并將這些基因構建為共表達網(wǎng)絡進行網(wǎng)絡層面的分析，成功地驗證COL1A2在胃癌中明顯高表達，并且這一特性不隨患者年齡、性別以及腫瘤前后期變化。

筆者以胃癌和食管癌為生物背景提出了一種新的方法，用于預測在基因調(diào)控網(wǎng)絡中的相關致病模塊以及致病模塊中的潛在癌癥標記物。具體過程如圖1所示，整合基因表達譜、miRNA表達譜、脫氧核糖核苷酸(DeoxyriboNucleic Acid，DNA)甲基化譜，對基因和miRNA進行初步篩選，結合miRNA、轉錄因子和基因三者之間的轉錄關系、蛋白質磷酸化修飾關系分別構建胃癌和食管癌的特異性基因調(diào)控網(wǎng)絡。對兩個網(wǎng)絡進行加權共表達分析，進行層次聚類和動態(tài)切割，不斷地調(diào)整優(yōu)化參數(shù)使其呈現(xiàn)出最優(yōu)的模塊性。借助于已知的癌癥標記物在網(wǎng)絡中的連通特性，計算每個模塊的第一主成分和已知癌癥標記物的相關性最大-最小差值，使用該值排名選出特異性的模塊，分析模塊內(nèi)部的生物通路和模塊之間的關系。在特異性的模塊中，使用5種不同的距離度量方法識別潛在的癌癥標記物。結果表明，在基于距離和基于相關性的兩類方法中，基于相關性的方法識別潛在癌癥標記物時準確率更高。具體而言，在胃癌和食管癌識別中，皮爾遜相關系數(shù)的結果優(yōu)于斯皮爾曼相關系數(shù)的結果。

圖1 胃癌與食管癌的癌癥標記物研究流程圖

1 實驗數(shù)據(jù)處理

1.1 已知癌癥標記物

已知的癌癥標記物數(shù)據(jù)來源于多個數(shù)據(jù)庫。從在線人類孟德爾遺傳數(shù)據(jù)庫(Online Mendelian Inheritance in Man，OMIM)和癌癥體細胞突變目錄數(shù)據(jù)庫(Catalogue Of Somatic Mutations In Cancer，COSMIC)中收集分別與胃癌和食管癌相關的基因，從miR2Disease、PhenomiR和人類miRNA疾病數(shù)據(jù)庫(Human Microrna Disease Database，HMDD)中收集分別與胃癌和食管癌相關的miRNA分子。將來源不同的數(shù)據(jù)進行合并后獲得的胃癌和食管癌相關基因分別有30個和17個，相關miRNA分別有381個和186個。這些基因和miRNA都將作為癌癥標記物在后續(xù)的網(wǎng)絡分析中發(fā)揮作用。

1.2 表達譜與甲基化譜

基因表達譜、miRNA表達譜和DNA甲基化譜均來源于癌癥基因圖譜(The Cancer Genome Atlas，TCGA)數(shù)據(jù)庫。對于基因表達譜和miRNA表達譜，將其涉及的樣本統(tǒng)一為相同的患者，并且要求同時擁有癌癥組織和癌旁正常組織樣本。除此以外，還需要對基因和miRNA表達譜數(shù)據(jù)進行數(shù)據(jù)清洗，對基因表達譜的清洗過程為去除在50%的樣本上表達值小于1的基因。經(jīng)過該過程，將胃癌基因表達譜的規(guī)模從20 531個基因和450個樣本降為17 059個基因和32對樣本；將食管癌基因表達譜的規(guī)模從20 531個基因和196個樣本降為17 150個基因和10對樣本。miRNA表達譜的清洗分為兩部分:對于已知的癌癥標記物miRNA,去除在50%的樣本上表達值缺失的分子;對于其余分子，則去除在90%的樣本上有缺失值的分子。篩選過的miRNA表達譜中仍然存在一部分缺失值，在本實驗中使用最小二乘插值法處理。經(jīng)過上述步驟，胃癌miRNA表達譜的規(guī)模為477個分子和32對樣本；食管癌miRNA表達譜的規(guī)模為471個分子和10對樣本。在DNA甲基化譜中，位點和分子之間是多對多的映射關系。在此，將甲基化位點對應為基因和miRNA分子，方法為均值法，即對每個分子所涉及的位點甲基化值取均值，以該值為對應分子的甲基化值。處理后的甲基化譜分為基因甲基化譜和miRNA甲基化譜兩個部分。胃癌的基因甲基化譜包含20 162個基因和137個樣本，胃癌miRNA甲基化譜包含1 504個miRNA和137個樣本。食管癌的基因甲基化譜包含20 164個基因和80個樣本，食管癌的miRNA甲基化譜包含1 507個miRNA和80個樣本。

1.3 調(diào)控數(shù)據(jù)與磷酸化數(shù)據(jù)

基因表達調(diào)控模式中的轉錄調(diào)控包含：miRNA對基因的調(diào)控、miRNA對轉錄因子的調(diào)控、轉錄因子對基因的調(diào)控以及轉錄因子對miRNA的調(diào)控。在構建基因調(diào)控網(wǎng)絡中，調(diào)控關系是重要的一環(huán)，它提供了轉錄因子、miRNA和基因之間的連接關系。從miRanda、PITA、TargetScan、基于句子的文本挖掘轉錄調(diào)控關系(Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining，TRRUST)、人類轉錄調(diào)控相互作用(Human Transcriptional Regulation Interactions，HTRI)、Mirtrans和TransmiR數(shù)據(jù)庫中收集了上述的4種調(diào)控關系。為了提升數(shù)據(jù)的可靠性和準確性，只選擇同時在兩個數(shù)據(jù)庫出現(xiàn)的miRNA-轉錄因子和miRNA-基因調(diào)控關系。在經(jīng)過該過程后，調(diào)控關系中miRNA-轉錄因子的調(diào)控關系數(shù)量是13 768，miRNA-基因的調(diào)控關系數(shù)量是123 493，轉錄因子-基因的調(diào)控關系數(shù)量是53 855，轉錄因子-miRNA的調(diào)控關系數(shù)量是7 036，轉錄因子-轉錄因子的調(diào)控關系數(shù)量是5 462。

蛋白質磷酸化關系與細胞中的代謝調(diào)節(jié)有關。磷酸化關系中代謝酶的濃度和反應時間都可以調(diào)節(jié)細胞的代謝。在PhosphoSitePlus數(shù)據(jù)庫中收錄了大量人類激酶的已知底物。除此以外，PhosphoNetworks數(shù)據(jù)庫也收錄了全面的磷酸化關系數(shù)據(jù)集，不僅包括從蛋白質微陣列分析得到的高通量數(shù)據(jù)，還有從權威文獻中選出的磷酸化關系對。將這兩個數(shù)據(jù)庫中的磷酸化關系進行匯總，共選出9 463條磷酸化關系。

2 網(wǎng)絡分析方法

首先，分別獲取在兩種癌癥中甲基化異常的分子。甲基化異常包括高甲基化和低甲基化兩種類型。某些DNA的低甲基化會導致致癌基因轉錄激活，另一些抑癌基因的高甲基化也被證實和它的低表達相關。為了構建的網(wǎng)絡更具特異性，選擇甲基化特異的分子作為該網(wǎng)絡的節(jié)點，即分子的甲基化值大于0.8或者小于0.2。然后，使用這些分子對調(diào)控關系和磷酸化關系進行選擇，選擇的具體條件如下：如果某一條調(diào)控關系或者磷酸化關系涉及的兩個分子中存在基因或者miRNA，則它必須是甲基化異常分子；如果某一條調(diào)控關系涉及的兩個分子中存在轉錄因子，則該分子必須在表達譜中存在表達值向量。經(jīng)過選擇的調(diào)控關系和磷酸化關系視為基因調(diào)控網(wǎng)絡的連邊。

經(jīng)過上述步驟，兩個特異性網(wǎng)絡的節(jié)點和連邊都已經(jīng)確定，網(wǎng)絡屬性為無權無向網(wǎng)絡。為了加強網(wǎng)絡的緊密性和準確性，對目前的網(wǎng)絡進行加權和進一步調(diào)整優(yōu)化連邊。結合網(wǎng)絡與表達譜，為網(wǎng)絡中現(xiàn)有的每條連邊計算兩個節(jié)點的斯皮爾曼相關系數(shù)。如果p<0.05且斯皮爾曼相關系數(shù)的絕對值大于0.3，則在網(wǎng)絡中保留該連邊，并且將該斯皮爾曼相關系數(shù)作為該連邊的權值。對所有的連邊都進行此操作后，去除了聯(lián)系不緊湊的連邊。如果某個節(jié)點和它的鄰居節(jié)點之間的連邊都被去除后，則該節(jié)點作為孤立節(jié)點也將被去除。兩個網(wǎng)絡均屬于無向有權網(wǎng)絡，且連邊權值均大于0.3。

加權基因共表達網(wǎng)絡分析(Weighted Gene Co-expression Network Analysis，WGCNA)一般用于分析多個樣本基因的表達模式，可以用于在網(wǎng)絡中尋找高度相關的類簇和模塊[4]。以拓撲重疊矩陣作為加權基因共表達網(wǎng)絡的輸入，相比于鄰接矩陣，拓撲重疊矩陣加入了一階相關的關系，提升了對網(wǎng)絡的描述準確度。拓撲重疊矩陣中元素wij的計算公式為

(1)

對構造的胃癌和食管癌基因調(diào)控特異性網(wǎng)絡進行加權基因共表達網(wǎng)絡模塊分割的方法分為兩個步驟：首先對網(wǎng)絡進行層次聚類，然后對層次聚類樹進行動態(tài)切割。在進行層次聚類之前，先將特異性網(wǎng)絡的鄰接矩陣變?yōu)橥負渲丿B矩陣，降低原網(wǎng)絡中的噪音與假相關，得到相異性拓撲重疊矩陣再進行層次聚類。在本實驗中，使用 average-linkage 方法進行聚類，計算兩個分組每個數(shù)據(jù)點與其他數(shù)據(jù)點的距離，將所有距離的均值作為兩個分組之間的距離。這種方法對每個類簇中的大部分分子都有代表性，可以得到有代表性的結果。設置每個模塊最少的分子數(shù)量為30個，對兩個網(wǎng)絡分別進行動態(tài)混合切割，一方面尋找所有滿足以下條件的模塊：模塊中分子個數(shù)滿足設定的最低數(shù)量，從模塊中移除距離過遠的分支，模塊由中心分子緊密連接；另一方面將沒有歸屬的分子分配進與它最近的初始模塊。最終，胃癌被分為14個模塊，食管癌分為10個模塊。

在已經(jīng)劃分為模塊的基礎上，選擇特異性較強的模塊進行重點分析。對每個模塊用第一主成分法計算模塊的特征向量基因(Module Eigengene，ME)。每個模塊的特征向量基因是由該模塊中分子表達值向量線性組合而來，其值與整個模塊內(nèi)的分子表達高度相關。計算整理的所有已知癌癥標記物與模塊特征向量基因的斯皮爾曼相關系數(shù)，使用最大最小法計算每個模塊的絕對相關系數(shù)的最大值與最小值之差。絕對值之差越大，表明該模塊在抑制癌癥標記物上越有差異性。筆者挑選每個網(wǎng)絡中排名前三的模塊并且對模塊內(nèi)的分子進行候選癌癥標記物排序。

候選癌癥標記物是指特異性模塊中的非已知癌癥標記物，將模塊的特征向量基因與所有候選癌癥標記物的表達向量進行相似性度量。某個分子與模塊的特征向量基因相關性越大，說明該分子越有可能是該模塊的中心分子。在每個特異性模塊中每種相似性度量方法產(chǎn)生一種排序，取每種方法排名前10的分子進行分析驗證。相似性評估方法可以分為兩類：基于距離的度量以及基于相關性的度量。在此使用的基于距離的度量方法包括歐幾里得距離、曼哈頓距離和切比雪夫距離；基于相關性的度量方法包括皮爾遜相關系數(shù)和斯皮爾曼相關系數(shù)。

3 實驗結果分析

將調(diào)控關系與磷酸化關系的范圍收縮至甲基化異常的分子范圍內(nèi)，構造的胃癌特異性基因調(diào)控網(wǎng)絡的規(guī)模為905個節(jié)點和1 597條連邊，包含67個miRNA和838個基因(包括轉錄因子)，其中45個miRNA 和4個基因均是胃癌分子標記物。食管癌特異性基因調(diào)控網(wǎng)絡規(guī)模為783個節(jié)點和1 309條連邊，包含68個miRNA和715個基因，其中18個miRNA 和1個基因是已知食管癌分子標記物。

通過最大最小法選擇的特異性模塊的特性、模塊中分子的數(shù)量和這些分子中的已知癌癥標記物數(shù)量如表1所示。對兩種疾病的6個模塊進行富集分析，包括京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes，KEGG)、基因本體(Gene Oncology，GO)、Canonical Pathways以及Reactome Pathway這4種形式的分析。

表1 胃癌與食管癌特異性模塊分子數(shù)量

表2表示了富集條目和特異性模塊的關系。經(jīng)分析發(fā)現(xiàn)：在ESCAmodule3和ESCAmodule1中，包含較多的與細胞中的轉錄和調(diào)節(jié)有關的富集通路。ESCAmodule1中顯著富集的是PERK介導的未折疊蛋白反應(GO:0036499)，該過程由腫瘤壞死因子引導激活。在ESCAmodule3中還有MRNA失活(R-HSA-429947)和RNA代謝(R-HSA-8953854)顯著富集，這二者同屬于信號反應通路。還有模塊ESCAmodule1、ESCAmodule3、STADmodule1和STADmodule2，都參與到了翻譯初始化的調(diào)節(jié)過程中(GO:0006446)。

表2 富集條目和特異性模塊對照關系

模塊的富集結果表明， 經(jīng)過動態(tài)切割后生成的模塊具有明顯的功能性。STADmodule3和ESCAmodule1與免疫系統(tǒng)和腫瘤的發(fā)展有關。STADmodule1和ESCAmodule2富集的通路基本上重疊，說明這兩個模塊涉及的功能也近似。STADmodule2和ESCAmodule3參與到了細胞的分裂和增殖階段，其中STADmodule2控制調(diào)節(jié)轉錄和翻譯過程，ESCAmodule3則在核分裂方面表現(xiàn)突出。

圖2 表示了每個特異性模塊之間的關系，連線表示模塊之間的基因分子重疊或者基因與基因之間存在功能相關。STADmodule2和ESCAmodule3有較多的分子重疊，包括FTSJ2、 CITED2和LRP10。而STADmodule1只和ESCAmodule2有分子重疊，比如DNAJB4，該基因已經(jīng)被證實在轉移性黑色素瘤、乳腺癌和非小細胞肺癌中作為癌癥標記物存在(PMID:29722524，PMID:28521415，PMID:28481734)。

圖2 特異性模塊關系圖

使用5種相似性度量方法，對胃癌的3個特異性模塊和食管癌的3個特異性模塊進行候選分子排序。將每個模塊中的每種方法排名前10的分子提取出來進行對比和數(shù)據(jù)庫驗證。除已經(jīng)被COSMIC和OMIM收錄的癌癥標記物以外，一些近期公開發(fā)表的研究中提及的在癌癥中特異表達的分子作為癌癥標記物也具有一定可信度。

圖3表示在胃癌的3個特異性模塊中使用5種不同相似性度量方法的前10名分子中，可以得到文本驗證的數(shù)量。

在STADmodule1中基于距離的方法得到驗證的分子完全相同，即每種方法都有2個分子驗證成功，且這2個分子同時出現(xiàn)在3種方法中。斯皮爾曼相關系數(shù)和皮爾遜相關系數(shù)方法各有4個分子得到驗證，其中3個分子重合。斯皮爾曼相關系數(shù)方法中預測出的TPR分子，文獻[5]采用TPR探針對胃癌患者進行分析，就已經(jīng)發(fā)現(xiàn)TPR對胃癌有抑制作用。而只在皮爾遜相關系數(shù)方法中出現(xiàn)的分子是UBQLN2，UBQLN2在原癌細胞中表達水平高于癌旁組織[6]。

STADmodule2中基于距離的方法和基于相關性的方法在結果上有明顯區(qū)分。歐幾里得距離和切比雪夫距離都只有一個相同的分子得到驗證，而曼哈頓距離在前兩者的基礎上多了SOCS1分子。在關于SOCS1的抗腫瘤機制研究中表明，SOCS1可以抑制80%的胃癌細胞系增殖，這意味著SOCS1可以作為一種新型治療胃癌的方法[7]。斯皮爾曼相關系數(shù)在STADmodule2中表現(xiàn)突出，有6個分子完全得到驗證，皮爾遜相關系數(shù)也有5個分子得到驗證，與斯皮爾曼相關系數(shù)方法比較只缺少了TIMM17B。

STADmodule3中基于距離的方法表現(xiàn)欠佳，3種方法均只有1個相同的分子hsa-miR-454-3p得到驗證。hsa-miR-454-3p在胃癌細胞系中通過介導HDAC3來調(diào)節(jié)胃癌細胞的生長[8]。而斯皮爾曼相關系數(shù)結果中有3個分子得到成功驗證，皮爾遜相關系數(shù)方法在STADmodule3中優(yōu)于所有方法，有4個分子驗證成功。

在ESCAmodule1中，3種基于距離的度量結果中均有2個分子得到驗證，分別是hsa-miR-191-3p和hsa-miR-29b-2-5p。文獻[9]通過聚合酶鏈式反應檢測在93例食管癌組織中hsa-miR-191-3p的表達情況，檢驗得知hsa-miR-191-3p的過度表達可以導致食管癌細胞的增殖和侵襲。文獻[10]也證實hsa-miR-29b-2-5p能抑制食管癌的增殖，可以作為食管癌的抑癌分子在預后治療中發(fā)揮作用。在ESCAmodule1的兩種基于相關性的度量方法中，各有5個分子得到驗證，并且這兩種方法的5個分子完全相同，但其排名順序不同。hsa-miR-191-3p出現(xiàn)在這5種排名中，表明該分子作為癌癥標記物的可能性極大。

圖3 特異性模塊驗證分子數(shù)量對比圖

在ESCAmodule2中，3種基于距離的度量方法得到的分子也完全相同，有2個分子可以得到驗證，排名不相同。HNRNPA2B1在基于相關性的兩種方法中都排在前列[11]。在皮爾遜相關系數(shù)的排名結果中還包括hsa-miR-191-5p，該分子與ESCAmodule1的5種方法中均預測出的hsa-miR-191-3p關系密切，兩個分子源于同一個miRNA前體hsa-miR-191[9]。

ESCAmodule3中基于距離的方法和基于相關性的方法驗證成功的分子數(shù)量相差甚小。歐幾里得距離和曼哈頓距離均有4個分子驗證成功，切比雪夫方法相較于前兩者并沒有成功預測出RAD51。RAD51可以通過自噬調(diào)節(jié)CHK1的穩(wěn)定性，促進食管癌細胞的生長[12]。斯皮爾曼相關系數(shù)和皮爾遜相關系數(shù)中也同樣包含RAD51，不過皮爾遜相關系數(shù)在斯皮爾曼相關系數(shù)的基礎上增加了H2AFX分子。H2AFX被證實在食管癌的細胞株中表達值顯著增加[13]。在ESCAmodule3中，RAD51在4種方法中都出現(xiàn)了，可以作為該模塊的核心候選癌癥標記物。

實驗結果顯示，在5個模塊中，基于距離的測量方法不如基于相關性的方法表現(xiàn)出彩。而在3種基于距離的測量方法中，曼哈頓方法在每個模塊中都不輸于其他兩種。在兩種基于相關性的測量方法中，皮爾遜相關系數(shù)法在5個模塊中優(yōu)于斯皮爾曼相關系數(shù)法。

4 討 論

筆者結合多種數(shù)據(jù)分別構建胃癌和食管癌的特異性網(wǎng)絡。使用轉錄因子和miRNA之間的相互調(diào)控關系、轉錄因子對基因的調(diào)控和miRNA對基因的調(diào)控，并考慮蛋白質的磷酸化修飾關系，以這些關系作為特異性網(wǎng)絡的連邊，以在胃癌中和食管癌中甲基化異常的分子作為網(wǎng)絡的節(jié)點，這樣的操作可以使網(wǎng)絡的規(guī)模精簡。甲基化在表觀遺傳學中至關重要，對于癌癥標記物研究有極大幫助[14]。筆者提出對特異性模塊內(nèi)部進行候選癌癥標記物預測，結合模塊的特征基因使用5種相似性度量方法計算候選分子的排名。比較5種相似性評估的方法，通過文獻驗證方法，評估這5種方法的準確率。實驗證明，在基于距離和基于相關性的方法中，基于相關性的方法表現(xiàn)更為出色。進一步，在兩種基于相關性的方法中，皮爾遜相關系數(shù)方法在此實驗的數(shù)據(jù)下表現(xiàn)優(yōu)于斯皮爾曼相關系數(shù)方法。在以后可以將類似的方法應用于其他的復雜疾病上，探討復雜疾病的分子機制。