基于改進TF-IDF算法的基因通路富集方法

徐淑坦，冷銀輝，陳明

1.上海海洋大學信息學院，上海 201306；2.農(nóng)業(yè)農(nóng)村部漁業(yè)信息重點實驗室，上海 201306

前言

在生物醫(yī)學相關研究領域，隨著高通量測序技術的發(fā)展，組學數(shù)據(jù)的規(guī)模也呈指數(shù)級增長。從龐大的組學數(shù)據(jù)中，可以利用生物信息學技術挖掘與疾病發(fā)生機制相關的通路，對疾病的診斷和治療具有重要意義。

在過去10 多年，已經(jīng)開發(fā)了很多基因功能富集分析方法來識別各種疾病相關的通路［1-3］?；跀?shù)據(jù)來源和算法大致可以將基因功能富集分析方法分為4 大類：過代表分析（Over-Representation Analysis,ORA）、功能集打分（Functional Class Scoring, FCS）、基于通路拓撲結(jié)構（Pathway Topology,PT）和基于網(wǎng)絡拓撲結(jié)構（Network Topology,NT）的方法［3］。ORA方法首先選定一組感興趣的基因作為基因列表，然后對該基因列表與通路中的基因集做交集，找出它們共同的基因并進行計數(shù)，最后利用統(tǒng)計檢驗的方式來評估計數(shù)值是否顯著高于隨機，即待測通路在基因列表中是否顯著富集?；蚣患治觯℅ene Set Enrichment Analysis, GSEA）屬于FCS 類方法，是最具代表性的基因功能富集分析方法，該方法利用基因表達和表型數(shù)據(jù)對所有基因進行排序，然后計算基因集Kolmogorov-Smirnov（KS）統(tǒng)計量，即在基因排序列表中靠近兩端程度的得分，最后通過置換方法評估基因集的顯著性［4］。FCS 類方法將通路中的基因視作獨立個體，實際上通路中的基因通過復雜的相互作用來影響細胞的發(fā)育、分化或疾病等生物學過程［5］。之后一些基于PT 的方法開始考慮通路的信息，Liu 等［6］提出一種基于定向隨機游走的方法來推斷通路活性，即利用基因在定向通路中的結(jié)構信息來評估每個基因在通路中的重要性，然后使用重要性對基因加權進行分析。Deng 等［7］利用蛋白質(zhì)互作數(shù)據(jù)和通路的基因集構建了蛋白質(zhì)-通路交互網(wǎng)絡，從全局層面優(yōu)化富集分析。Yang等［8］提出一種基于PT 的通路富集方法，該方法根據(jù)基因節(jié)點的全局上游或下游位置和通路中的連接度數(shù)來評估節(jié)點的重要性，富集出那些擁有更多在上游或中樞節(jié)點中的基因的通路。后來的基于NT 的方法中，很多方法也借鑒了GSEA 的思想。Winterhalter 等［9］提出基于蛋白質(zhì)互作網(wǎng)絡的GSEA 方法JEPETTO，利用網(wǎng)絡的拓撲結(jié)構信息分析基因和通路之間的功能關聯(lián)。Rahmati等［10］整合來自20個通路數(shù)據(jù)庫的核心通路，構建了龐大的蛋白質(zhì)互作網(wǎng)絡，以求富集出與生物學功能相關的通路，提出pathDIP 方法。Han 等［11］提出基于人類基因功能網(wǎng)絡的GSEA 方法NGSEA，該方法衡量基因集的富集分數(shù)，不僅考慮單個基因的表達差異，還考慮功能網(wǎng)絡中它們的相鄰基因的表達差異。Yoon 等［12］提出一種新的網(wǎng)絡加權的GSEA，對基因集進行富集分析時，將重疊基因和蛋白質(zhì)互作網(wǎng)絡結(jié)合起來，有效地識別出與基因功能相關的基因集。Zito 等［13］將圖論中衡量結(jié)點的中介中心性權重引入通路富集分析，利用蛋白質(zhì)互作網(wǎng)絡，提出網(wǎng)絡節(jié)點中心性加權的GSEA。從更高層面看，通路是生物互作網(wǎng)絡的一部分，因此通路富集分析有必要綜合考慮基因在通路的局部信息和在通路數(shù)據(jù)庫的全局信息，這一點和數(shù)據(jù)挖掘領域的常用加權算法TF-IDF（Term Frequency-Inverse Document Frequency）的思想類似，該算法綜合考慮字詞在文件的局部重要性和文件庫的全局特異性兩方面來評估字詞的權重。

針對以上問題，本研究融合基因在通路的局部重要性和在通路數(shù)據(jù)庫的全局特異性定義基因影響力，提出一種基于改進TF-IDF 算法的基因通路富集方法（GIGSEA 方法）。首先利用基因相互作用數(shù)據(jù)來計算通路基因的影響力，然后通過基因表達數(shù)據(jù)和表型數(shù)據(jù)來計算基因與表型的相關性值，并利用基因影響力和相關性值計算富集分數(shù)，最后通過統(tǒng)計學方法計算通路的顯著性P值。在肝細胞癌（Hepatocellular Carcinoma, HCC）和結(jié)腸直腸癌（Colorectal Cancer, CRC）數(shù)據(jù)集的實驗結(jié)果表明該方法能有效識別出與疾病相關的通路，對于今后研究疾病的發(fā)生發(fā)展機制有重要的指導意義。

1 材料與方法

1.1 數(shù)據(jù)集及預處理

本研究的通路數(shù)據(jù)來自KEGG 數(shù)據(jù)庫［14］，KEGG數(shù)據(jù)庫是代謝組學和蛋白質(zhì)組學研究中常用的生化通路數(shù)據(jù)庫，從GSEA 網(wǎng)站下載通路數(shù)據(jù)集，共包括186個通路，共有5 245個基因。

本研究的基因相互作用數(shù)據(jù)從STRING 數(shù)據(jù)庫［15］中獲得，下載地址為https://www.string-db.org/cgi/download，下載的數(shù)據(jù)是最新的11 版本的人類蛋白互作數(shù)據(jù)，共11 938 498 對，從中提取包含KEGG通路基因的基因相互作用，去除重復的相互作用對，共計977 800對。

本研究的基因表達數(shù)據(jù)包括HCC 和CRC 數(shù)據(jù)集。HCC 數(shù)據(jù)集來自TCGA 數(shù)據(jù)庫，包括來自肝癌患者的腫瘤鄰近組織的50 對配對樣本和324 個腫瘤樣本，本研究只使用50對配對樣本進行測試，即包括50 個腫瘤樣本和50 個正常樣本，該數(shù)據(jù)集下載自UCSC Xena 網(wǎng)站（https://xena.ucsc.edu/）［16］。CRC 數(shù)據(jù)集來自GEO（Gene Expression Omnibus）數(shù)據(jù)庫的GSE8671 表達譜數(shù)據(jù)，包括32 個結(jié)腸癌腺瘤組織和32個配對的癌旁組織。

1.2 方法

1.2.1 基于改進TF-IDF 算法計算基因影響力在TFIDF 算法中，TF 是詞頻（Term Frequency），即詞條t 在文件中出現(xiàn)的頻率，TF 越大，意味著t 可能是文件的關鍵詞，說明t 越重要；IDF 是逆文件頻率指數(shù)（Inverse Document Frequency），如果包含詞條t 的文件越少，IDF 越大，說明t 具有很好的類別區(qū)分能力，t的權重相應地越大。TF是對詞條t在本文件的評估，IDF是對t在文件庫的評估。

類似地，本研究的基因影響力由基因在本通路的局部重要性和在通路數(shù)據(jù)庫的全局特異性來衡量。通路基因的相互作用的關系可以抽象為一個圖，本研究定義基因的重要性為在通路中與其他基因產(chǎn)生相互作用的數(shù)量，即在通路圖中的連接度。采用基因頻率GF（Gene Frequency）來表示基因重要性，GF表示為：

其中，ni,j是基因j在通路pi中的度數(shù)，分母則是通路pi中所有基因的度總和。

一個基因的特異性表現(xiàn)為基因在所有通路中出現(xiàn)的頻度，頻繁出現(xiàn)在很多通路中的基因，它們對通路的影響相對較??；僅在少數(shù)通路中出現(xiàn)的基因其特異性高，它們的差異表達對通路的影響就大。本研究定義逆通路頻率（Inverse Pathway Frequency,IPF）來表示基因特異性，IPF表示為：

其中，|P|是數(shù)據(jù)庫中的基因通路的總數(shù)；|i:gj∈pi|表示包含基因gj的通路數(shù)目，即ni,j≠0的通路數(shù)目。

某基因在特定通路內(nèi)的重要性越高，并且在通路數(shù)據(jù)庫中的特異性越高，那么該基因影響力越大。合并GF和IPF計算基因影響力GI(gi,j)，表示為：

其中，GI(gi,j)表示基因j對通路i的影響力大小。

1.2.2 GIGSEA 方法假定有N個基因、T個樣本的基因表達數(shù)據(jù)，T個樣本包括兩種表型pos 和neg，樣本數(shù)量分別為t1 和t2；給定一個通路pi，包括Q個基因。GIGSEA方法主要過程如下：

（1）利用基因相互作用數(shù)據(jù)統(tǒng)計每個基因在通路中的連接度數(shù)量，計算出GF（圖1a）。然后統(tǒng)計每個基因在全部通路中出現(xiàn)的頻數(shù)，計算出基因IPF。最后合并GF 和IPF計算GI(gi,j)。

（2）考慮到通路中有一些基因不在基因表達數(shù)據(jù)的基因集中，對兩個集合取交集，計算交集內(nèi)的基因與表型的相關性值，計算公式為：

其中，函數(shù)mean(x)、std(x)分別指表型x的基因表達值的平均值和標準差。最后按相關性值從大到小排序得到列表L=［［g1,r1］,…,［gj,rj］,…,［gm,rm］］，包括M個基因（圖1b）。

（3）計算通路的富集分數(shù)。從列表L 的第一個基因開始，當遇到一個在通路pi里面的基因（hits），則增加分數(shù)；遇到一個不在pi里面的基因（misses），則減少分數(shù)，具體公式為：

最終得到一個分數(shù)曲線（圖1c），曲線上的點到橫坐標距離的最大值即為ES0(pi)。

（4）隨機置換基因Nperm次（實驗中，Nperm取基因富集分析常用的值，1 000次）。隨機置換基因指隨機在列表L 挑選Q個基因作為通路基因。重復步驟（3），計算置換后的通路富集分數(shù)ESperm(pi)。如圖1d所示，最后統(tǒng)計|ESperm(pi)|＞|ES0(pi)|的數(shù)量Nsign，P值等于Nsign與Nperm的比值。

圖1 GIGSEA方法流程圖Figure 1 GIGSEA flowchart

為了方便比較，本研究按GSEA 方法［4］計算校正后的富集分數(shù)（Normalized Enrichment Score, NES）和錯誤發(fā)現(xiàn)率（False Discovery Rate,FDR）。

1.3 結(jié)果評價

本研究將顯著性P值、FDR 和|NES|的閾值分別設為0.05、0.25 和1.00，篩選出具有顯著性意義的通路，為了檢驗方法的有效性，與GSEA方法進行比較。兩種方法的顯著通路大部分是重疊的，重疊通路只是排名的略微差異，因此關注顯著通路的差集更有意義，將顯著通路差集內(nèi)的通路定義為差異通路。本研究主要從差異通路的3 個方面來說明方法的有效性，包括（1）生物學解釋驗證。通過查閱關于HCC或CRC 相關的生物學研究文獻來驗證通路與HCC或CRC 存在的某種聯(lián)系。大部分富集方法都通過這種方式來解釋，如果通路中的某些基因或者產(chǎn)物對疾病產(chǎn)生影響，那么該通路與疾病是相關的。（2）相關文獻數(shù)量。通過PubMed 生物醫(yī)學論文數(shù)據(jù)庫檢索通路和HCC 或CRC 存在聯(lián)系的文獻，利用文獻數(shù)量的多少來表示差異通路與疾病相關性的強弱。如果檢索詞全部出現(xiàn)在文獻中，那么該條文獻會出現(xiàn)在PubMed 檢索結(jié)果中，因此文獻數(shù)量從一定程度上可以反映通路與疾病的相關性。比如對于Jak Stat Signaling Pathway 與HCC 的相關性，通過關鍵詞Jak Stat、Hepatocellular Carcinoma 來進行搜索，排除掉Signaling Pathway 這些冗余的詞；而對于Asthma 和Peroxisome 通路，則采取補全關鍵詞Pathway 來檢索。（3）通路基因集對應的表達數(shù)據(jù)對表型的分類性能。每一個通路基因集對應的基因表達數(shù)據(jù)和樣本表型數(shù)據(jù)輸入到構建的支持向量機（Support Vector Machine,SVM）分類模型中，求出AUC，結(jié)合計算的P值、FDR 和||NES ，利用分類性能對兩種方法進行有效性比較［17］。AUC 是分類模型中常用的評價指標，其取值范圍為0.5～1.0，AUC 越大，分類效果越好，意味著通路基因的表達值對疾病分類越準確。比如HCC 數(shù)據(jù)集包括100 個樣本，有Normal 和Cancer 兩種表型，各50個；而通路Jak Stat Signaling Pathway中有155 個基因，則輸入SVM 模型的數(shù)據(jù)為：100 個樣本，每個樣本有155 個特征，對應155 個基因，每一個樣本的標簽對應于表型，比如Normal 為0、Cancer 為1。因此這是一個針對樣本表型的二分類問題，如果兩種表型的基因表達值存在的差異更高，則更能夠把表型區(qū)別開來，產(chǎn)生更高的AUC，意味著通路內(nèi)的基因和HCC疾病的相關性更高。

2 結(jié)果與討論

2.1 HCC數(shù)據(jù)集的結(jié)果與討論

HCC 數(shù)據(jù)集的富集結(jié)果如表1 所示，表中共有33個通路，為兩種方法顯著通路的并集。GSEA 富集出30 條通路，GIGSEA 富集出29 條通路，且富集出3個新的通路：Jak Stat 信號通路（Jak Stat Signaling Pathway）、過氧化物酶體通路（Peroxisome）、半胱氨酸和蛋氨酸代謝通路（Cysteine and Methionine Metabolism），分別排在第5、8、11 位，與GSEA 相比，它們的排名也都提高了。

表1 HCC數(shù)據(jù)集中兩種方法的富集結(jié)果Table 1 Enrichment results of two methods in HCC dataset

經(jīng)過查閱文獻發(fā)現(xiàn)，GIGSEA 的差異通路都與HCC疾病的發(fā)生機制有一定的聯(lián)系。Jak Stat信號通路的失調(diào)可能導致包括HCC 在內(nèi)的各種癌癥［18］。Tang 等［19］發(fā)現(xiàn)Jak Stat 信號通路在HCC 中維持具有腫瘤增殖能力的癌癥干細胞以及創(chuàng)建免疫抑制微環(huán)境，該通路中的STAT3 基因?qū)Π谢蚝偷鞍踪|(zhì)的調(diào)節(jié)促成腫瘤發(fā)生的概率。目前針對HCC，已經(jīng)開發(fā)出多種JAK 或STAT 小分子抑制劑和RNA 療法。過氧化物酶體通路中，過氧化物酶體包含至少50 種不同的酶［20］。Xu等［21］發(fā)現(xiàn)過氧化物酶體增殖物激活受體δ（PPARδ）是一種核轉(zhuǎn)錄因子，與腫瘤發(fā)生有關，通過PPARδ 和前列腺素信號通路之間的串擾，兩個通路共同調(diào)節(jié)人體HCC 細胞的生長。Wirtz 等［22］及Zhuang 等［23］發(fā)現(xiàn)HCC 細胞轉(zhuǎn)移與半胱氨酸和甲硫氨酸代謝通路中代謝物水平和代謝酶表達的改變有關，該通路的ENOPH1 的過表達促進細胞遷移和侵襲，而ENOPH1 下調(diào)抑制細胞遷移和侵襲，研究表明ENOPH1 可以促進HCC 進展，可以作為HCC 的生物標志物和治療靶點。

相關文獻可以反映已有研究對于通路和疾病具有相關性的支持，因此本研究對兩種方法的差異通路進行相關文獻統(tǒng)計，驗證其與HCC 的相關性。統(tǒng)計結(jié)果如表2 中文獻數(shù)量一列所示，GIGSEA 和GSEA 的差異通路的文獻數(shù)量均值分別為129 和6，GIGSEA 明顯高于GSEA，表明GIGSEA 的差異通路與HCC的相關性是更高的。

接著利用SVM 模型來測試差異通路對HCC 數(shù)據(jù)集的兩種表型的分類性能。GIGSEA 的平均P值和平均FDR 遠小于GSEA，且||NES 大于GSEA，從統(tǒng)計學上表明GIGSEA 差異通路與HCC 相關性更強。兩種方法的差異通路在SVM模型都表現(xiàn)出良好的分類效果，GIGSEA 和GSEA 的AUC 分別達到99.22%和96.99%，表明GIGSEA 差異通路的分類性能同樣優(yōu)于GSEA差異通路（表2）。

表2 HCC數(shù)據(jù)集差異通路的對比Table 2 Comparison of the differential pathways in HCC dataset

2.2 CRC數(shù)據(jù)集的結(jié)果與討論

CRC 數(shù)據(jù)集的富集結(jié)果如表3 所示，表中共有26 個通路，包括GIGSEA 和GSEA 的顯著通路。GSEA 富集出20 條通路，GIGSEA 富集出23 條通路，且富集出6個新的通路：肌動蛋白細胞骨架的調(diào)節(jié)通路（Regulation of Actin Cytoskeleton）、白細胞跨內(nèi)皮遷移通路（Leukocyte Transendothelial Migration）、補體和凝血級聯(lián)通路（Complement and Coagulation Cascades）、朊病毒病通路（Prion Diseases）、哮喘通路（Asthma）、腎素血管緊張素系統(tǒng)通路（Renin Angiotensin System），它們的排名在GIGSEA 中都有一定的提高。

表3 CRC數(shù)據(jù)集中兩種方法的富集結(jié)果Table 3 Enrichment results of two methods in CRC dataset

文獻檢索顯示除了哮喘通路，GIGSEA 的差異通路都與CRC 有一定的聯(lián)系。Kanaan 等［24］發(fā)現(xiàn)肌動蛋白細胞骨架調(diào)節(jié)通路通過細胞骨架蛋白，如Fascin-1，參與轉(zhuǎn)移性散發(fā)性CRC 的發(fā)展；與散發(fā)性CRC 相比，該通路的調(diào)控基因之間的相似遺傳多態(tài)性和突變也可能與CRC 的發(fā)育不良、癌變以及侵襲和轉(zhuǎn)移的易感性增加有關。Tremblay 等［25］發(fā)現(xiàn)E-選

擇蛋白被結(jié)腸癌細胞激活會觸發(fā)p38和ERK MAP激酶的激活，從而誘導細胞骨架重塑，導致內(nèi)皮層破裂，促進粘附的結(jié)腸癌細胞的外滲，該研究表明E-選擇蛋白介導的p38和ERK MAP激酶激活對結(jié)腸癌細胞跨內(nèi)皮遷移的調(diào)節(jié)。Matilda等［26］發(fā)現(xiàn)補體和凝血級聯(lián)通路、參與脂質(zhì)代謝的通路、急性期反應信號通路是高C 反應蛋白（C-reactive Protein, CRP）CRC 患者的主要干擾通路。細胞朊病毒蛋白（Cellular Prion Protein,PrPc）是一種細胞表面蛋白，由朊病毒通路中的PRNP 基因編碼［27］。Ong 等［28］發(fā)現(xiàn)PrPc 的過表達可能通過誘導內(nèi)皮增殖-分化開關的方式參與CRC誘導的血管生成。Chen等［29］發(fā)現(xiàn)腎素血管緊張素系統(tǒng)抑制劑的使用與CRC 風險和死亡率降低有關，該研究通過實驗證明腎素血管緊張素系統(tǒng)抑制劑使用持續(xù)時間每增加一年，CRC風險降低6%。

同樣地，通過Pubmed檢索和SVM模型分類來驗證差異通路與CRC 之間的相關性，結(jié)果如表4 所示。GIGSEA 差異通路的相關文獻數(shù)量的均值為55，而GSEA 為34；GSEA 的P值要小于GIGSEA，兩種方法的差異通路平均FDR 和||NES 都接近；GIGSEA 的平均AUC 明顯高于GSEA，分別為91.32%和85.67%，表明GIGSEA 的差異通路的分類性能優(yōu)于GSEA。綜合考慮，GIGSEA 的差異通路與CRC 的相關性比GSEA強。

表4 CRC數(shù)據(jù)集差異通路的結(jié)果對比Table 4 Comparison of the differential pathwaysin CRC dataset

2.3 基因通路富集分析網(wǎng)站

為了方便用戶使用GIGSEA進行通路富集分析，本研究基于SSM框架開發(fā)了一個在線的基因通路富集分析的網(wǎng)站，SSM 框架是Java 企業(yè)級開發(fā)領域Spring、Spring MVC 和MyBatis 框架的縮寫。本研究利用bootstrap-tablejs組件來繪制富集分析的結(jié)果展示頁面，可以對富集結(jié)果進行排序、搜索等功能（圖2）。

圖2 富集結(jié)果展示頁面Figure 2 Enrichment result display interface

本研究利用EChartsjs 組件的關系圖來進行基因通路的可視化，如圖3 所示，通路圖中的節(jié)點表示基因，節(jié)點之間的連線表示基因之間的相互作用關系，基因節(jié)點的大小與基因相關性值的絕對值成比例，基因節(jié)點的顏色與基因列表相關性值正負相關，紅色代表相關性值為正，藍色代表相關性值為負，灰色代表該基因不在基因列表中。

圖3 Echart可視化的通路局部示意圖Figure 3 Partial schematic diagram of Echart visualized pathways

3 結(jié)論

本研究提出一種基于改進TF-IDF 算法的GIGSEA 方法。首先利用通路基因相互作用數(shù)據(jù)，考慮基因在通路的局部重要性和在通路數(shù)據(jù)庫的全局特異性，計算基因的影響力；然后利用基因表達數(shù)據(jù)和表型數(shù)據(jù)計算基因與表型的相關性值；接著融合基因影響力和表型相關性值，計算通路的富集分數(shù)；最后通過置換基因的方式，考察通路是否和疾病相關。本研究利用HCC和CRC數(shù)據(jù)集來測試GIGSEA的效果。與GSEA 比較，本研究發(fā)現(xiàn)了與HCC 相關的3個新通路，以及與CRC 相關的6個新通路。除了哮喘通路，本研究都找到研究文獻來證實通路與疾病之間的相關性。利用PubMed 檢索相關文獻的結(jié)果顯示在兩個數(shù)據(jù)集中，GIGSEA 的文獻數(shù)量都遠遠多于GSEA。利用SVM 模型分類的結(jié)果顯示在兩個數(shù)據(jù)集中，GIGSEA 通路對應的表達數(shù)據(jù)的分類效果都優(yōu)于GSEA。GIGSEA 方法不僅豐富了富集分析方法，更重要的是為發(fā)現(xiàn)與疾病相關的通路提供了一種新思路。