整合生物信息學鑒定宮頸高級別鱗狀上皮內病變的關鍵基因

高銘，鄭沾福，劉興華，林泳煌，周萍，何玉清，3

(1.廣東醫科大學公共衛生學院流行病與衛生統計學教研室，廣東 東莞 523808；2.東莞市婦幼保健院兩癌篩查中心，廣東 東莞 523000；3.廣東醫科大學寮步醫院皮膚科，廣東 東莞 523400)

宮頸癌是世界第四大常見的女性惡性腫瘤，其發病率和死亡率極高。2018年有56.9萬例宮頸癌新發病例，占當年所有新發癌癥病例的3.2%[1]。目前已明確高危型人乳頭瘤病毒持續或反復感染是宮頸鱗狀上皮內瘤變(cervical intraepithelial neoplasia，CIN)以及宮頸癌發生發展的主要原因，大多數宮頸人乳頭瘤病毒感染能自行清除，但少數持續性感染會導致嚴重的不典型增生，最終導致浸潤性癌癥[2]。按照1～3等級對CIN進行分級，CIN 1級相當于宮頸低級別鱗狀上皮內病變，有發展為≥CIN 2級的風險；CIN 2級和CIN 3級則是高級別鱗狀上皮內病變(high-grade squamous intraepithelial lesions，HSIL)，≥CIN 2級有發展為宮頸癌的風險，CIN 3級是可能發展為宮頸浸潤癌的等級[3]。目前HISL的主要治療方式為盡可能保留生育功能的宮頸錐切術，但術后仍有較高的復發風險。研究表明，CIN 2～3級治療后約17%的婦女有殘留或可能發展為復發性CIN 2～3級[4-6]，還有部分HSIL發展為宮頸浸潤癌。因此，鑒定新的生物標志物以及早期篩查診斷尤為重要，可以極大地降低患者后期惡化情況，提高患者生活質量，保護患者的心理健康。

鑒定基因特異性表達模式可用于了解CIN、宮頸癌在內的多種疾病的致病機制，為診斷和治療預后評估提供新策略[7-8]。魯棒秩聚合法(robust rank aggretation，RRA)已用于基于多種疾病(如癌癥、自身免疫性疾病)的多個數據集篩選差異表達圖譜[9-11]，RRA是一種利用概率模型對排序列表進行整合的方法[12]，已有多項研究表明其是一種可靠的多組芯片數據集整合分析方法[11,13-15]。本研究旨在采用生物信息學的手段分析多數據集獲得差異表達基因并構建互作網絡，鑒定與HSIL發生發展密切相關的核心基因，進而識別新的生物標志物以及潛在的藥物治療靶點，為HSIL的研究提供新思路。

1 資料與方法

1.1數據來源 從GEO(Gene Expression Omnibus)數據庫(http://www.ncbi.nlm.nih.gov)下載HISL的基因表達數據。以“CIN”“HISL”“Gene expression”“Homo sapiens ”“Microarray”為檢索詞搜索微陣列研究。HSIL的基因表達數據集標準篩選：①實驗樣本均來源于宮頸組織；②HSIL組和對照組(正常宮頸組織)的樣本數量均≥5個；③GEO提供了原始數據或通過陣列能進行基因表達譜分析。

1.2數據集的預處理和差異性基因表達分析 從GEO數據庫下載每個陣列數據集的基因表達矩陣和相關注釋文檔，使用對應的注釋文檔將微陣列探針映射到基因符號。如果多個探針映射到同一符號，則采用平均值表示。為消除樣品之間的個體差異，將4個微陣列數據集進行分位數標準化。采用“limma”(用于微陣列數據的線性模型)R包進行分析，獲得每個微陣列中宮頸HISL組織與正常對照宮頸組織之間的差異表達基因列表，保存上調和下調的基因列表，確定差異表達基因的閾值標準為|log2FC|>1和校正后P<0.05。

1.3RRA整合分析 采用RRA識別可靠的差異基因。在進行RRA分析前，獲得每個數據集的上調和下調基因列表，這些基因列表是由HSIL組和對照組的表達倍數變化產生,然后通過“RobustRankAggregation”R包比較多個排名的基因列表。RRA分析中校正后的P值表明每個基因在最終基因列表中排名靠前的可能性。P<0.05和差異倍數變化>0.5 的基因則為差異顯著基因。

1.4信號通路富集分析 利用R包clusterProfiler[16]及ReactomePA[17]分別對差異基因進行基因本體(Gene Ontology，GO)功能注釋及京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes，KEGG)通路富集分析，以P<0.05和Q<0.05為標準篩選顯著的功能注釋及通路富集結果。

1.5蛋白互作網絡(protein protein interaction network，PPI network)構建與分析 使用STRING數據庫(http://www.string-db.org,version:11.0)，置信度標準為>0.4構建PPI網絡。采用Cytoscape(version:3.71)軟件對PPI網絡可視化，使用MCODE插件鑒定整個網絡中的核心模塊，并用Cytohubba插件識別整個網絡的關鍵基因(Hub基因)。在Cytoscape中，每個節點代表一個基因或蛋白質，節點之間的邊緣代表分子之間的相互作用。

2 結 果

2.1基因芯片信息 根據先前確定的納入標準，篩選出了4個符合本研究的芯片數據，包括GSE63514、GSE7803、GSE26278、GSE138080，每個數據集的具體信息見表1。本研究共納入分析161個樣品，其中包含56個正常對照和105個HSIL樣本。

2.2差異性基因表達分析 根據截止標準(log2FC>1，校正后的P<0.05)，按照表1的分組進行差異表達分析，使用“limma”R軟件包篩選出差異表達基因。結果顯示：GSE63514數據集獲得了1 956個差異表達基因(1 315個上調基因，641個下調基因)；GSE27678數據集獲得了1 022 個差異表達基因(567個上調基因，455個下調基因)；GSE7803數據集獲得了747個差異表達基因(372個上調基因，375 個下調基因)；GSE138080數據集獲得了1 023個差異表達基因(344個上調基因，679個下調基因)。4個微陣列的火山圖見圖1。

2.3RRA綜合分析的結果 使用“Robust Rank Aggeretation”R包分析排序好的4個微陣列的差異表達基因，確定了175個重要的差異表達基因(上調77個，下調98個)，其中前10個上調和前11個下調基因的熱圖見圖2。

2.4信號通路富集分析 通過clusterProfiler[16]包對175個差異表達基因進行GO富集分析，結果顯示主要富集在細胞分化、免疫反應、炎癥、代謝等功能類別上，見圖3。KEGG富集結果顯示，差異表達基因主要富集在補體和凝血級聯、花生四烯酸代謝、炎癥、前列腺癌、膀胱癌、細胞周期、白細胞介素-17、腫瘤壞死因子、AMP活化蛋白激酶信號通路上，見圖4。

GO：基因本體

KEGG：京都基因與基因組百科全書

2.5PPI網絡構建以及Hub基因鑒定 對175個基因進行PPI網絡構建，包含132個節點，421條邊，見圖5a。通過Cytoscape軟件進行可視化，使用MOCE插件提取1個關鍵模塊，得到10個關鍵基因：基質金屬蛋白酶9(matrix metalloproteinase-9，MMP-9)、金屬蛋白酶2組織抑制物(tissue inhibitor of metalloproteinase 2，TIMP2)、CXC趨化因子配體1[chemokine(C-X-C motif)ligand 1，CXCL1]、胞嘧啶脫苷酶(cytidine deaminase，CDA)、親珠蛋白(haptoglobin，HP)、乳鐵蛋白(lactotransferrin，LTF)、嗅覺素-4(olfactomedin 4，OLFM4)、富含半胱氨酸的分泌蛋白3(cysteine-rich secretory protein 3，CRISP3)、胱抑素B(cystatin B，CSTB)、胰蛋白酶3(protease,serine 3，PRSS3)，見圖5b。應用Cytohubba插件中MCC和DMNC兩種算法提取連接最多的前10個基因并取其交集，得到4個核心基因：LTF、CSTB、CRISP3、CDA，見圖5c。其中CXCL1、LTF在RRA法中TOP10的差異表達基因中。綜合分析，最終確定5個核心基因：細胞周期蛋白依賴性激酶抑制劑2A(cyclin-dependent kinase inhibitor 2A，CDKN2A)、富半胱氨酸的C端1(cysteine rich C-terminal 1，CRCT1)、CXCL1、LTF、CDA。

5a：Cytoscape可視化PPI蛋白互作網絡(其中紅色為上調基因，綠色為下調基因)；5b：提取的基因模塊；5c：Cytohubba兩種算法提取前10個基因的韋恩圖；PPI：蛋白互作網絡

3 討 論

高危型人乳頭瘤病毒持續感染是導致宮頸低級別鱗狀上皮內病變、HSIL和浸潤性宮頸癌發生的主要原因。HISL如果不能被早期發現和治療可能繼續發展為浸潤性宮頸癌。隨著宮頸癌篩查的推廣以及宮頸癌疫苗的接種，目前宮頸癌的發病率有所下降，但HSIL治療后的復發仍是臨床較難解決的問題。因此，對HSIL惡變前的診治顯得尤為重要。本研究的目的是篩選和鑒定與HISL密切相關的基因生物標志物，為臨床早期篩查和治療預后提供科學的理論依據。

PPI分析中節點之間的連接能被可視化，以鑒定HSIL中差異表達基因編碼的蛋白質之間的相互作用，位于中心節點的基因被認為是可能起關鍵作用的Hub基因，Hub基因通常被認為是參與整個調控網絡并有生物學功能的關鍵基因。本研究通過RRA聯合網絡拓撲異構結構對多個芯片進行綜合分析，成功鑒定了5個核心基因，其中CDKN2A、CXCL1兩個基因已有研究發現[18-21]與HSIL有直接或間接的關聯，CDKN2A、CXCL1可能通過細胞周期、p53、白細胞介素-17和腫瘤壞死因子等癌癥相關信號來靶向調控HSIL的發展進程。CDKN2A也稱為p16INK4A，是一種細胞周期蛋白依賴性激酶(cyclin-dependent kinase，CDK)抑制劑，可阻斷CDK4和CDK6介導的轉錄共加壓子的磷酸化，抑制E2F依賴的轉錄和細胞周期進程[22]。多項研究證實，幾乎在所有HSIL病例中均可以檢測到CDKN2A，其表達水平與宮頸進展的病變程度呈正相關[19]，具有較高的檢測特異性和敏感性，可作為HSIL的早期篩查和預后標志物[20-21]。CXCL1是一類具有促細胞分裂和促血管生成活性的趨化因子,常與炎癥、感染、細胞免疫相關，在多種癌癥中表達上調[23-24]，可介導腫瘤細胞的生長、增殖、存活、新血管生成和轉移。研究發現，CXCL1在不同級別的宮頸病變(CIN 1、2、3級)中差異表達，且其表達水平與病變級別呈正相關，可作為檢測HSIL的生物標志物，聯合其他趨化因子可增強預測的特異性和敏感性[18]。本研究結果與前期文獻報道[22-24]的結果一致，證明研究結果的可靠性。

CDA是嘧啶挽救途徑中的一種酶，可參與游離嘧啶的循環利用來維持細胞嘧啶池的穩定[25]，研究發現CDA在大約60%的癌細胞和組織(包括子宮頸癌)中表達下調[26]，本研究與上述研究結果一致。多項研究發現，細胞(包括HeLa細胞)中CDA缺乏可導致DNA損傷和細胞遺傳不穩定，此過程與腫瘤的發生發展密切相關[27-29]。結合現有研究，推測CDA的表達失調可能調控宮頸上皮內病變甚至惡變，但其分子機制還需進一步驗證。LTF是一種鐵結合糖蛋白，可參與體內鐵的轉運和代謝、調節機制免疫，具有抗炎、抗微生物以及抗腫瘤等生物學活性[30-31]。有研究發現，LTF在人乳頭瘤病毒感染的早期發揮抑制作用[32]，在炎癥、惡變前表達上調，但在惡變后以及癌癥中表達下調[33]；另有研究發現來源于中性粒細胞的LTF可以促進炎癥細胞因子和趨化因子(白細胞介素-6、CC趨化因子配體20等)的表達[34]。LTF可能通過炎癥通路參與HSIL的發生發展。CRCT1由表皮分化復合物編碼，在食管鱗狀細胞癌中表達下調，CRCT1過表達可促進食管鱗狀細胞癌細胞凋亡并上調凋亡相關蛋白的表達[35]。目前LTF、CRCT與HSIL發生發展有直接關聯的相關研究較少，因此仍需要進行深入研究。本研究的信號通路功能富集分析結果顯示，篩選出來的差異性表達基因主要富集在炎癥和癌癥相關的信號通路上，其中補體和凝血級聯、花生四烯酸代謝途徑是最顯著富集的3個炎癥相關通路。其中花生四烯酸代謝途徑主要用于炎癥介質的合成，可介導各種炎癥細胞因子(如單核細胞趨化蛋白1、腫瘤壞死因子、白細胞介素和干擾素)的產生與發展[36]。GO功能富集分析顯示，這些基因主要參與中性粒細胞相關的細胞分化、免疫等生物學過程，目前公認的慢性中性粒細胞炎癥涉及許多類型的上皮癌的起始階段，可誘發惡變。越來越多的證據顯示，炎癥與癌癥的發生發展密切相關[37-38]，表明本研究鑒定的關鍵基因可能與HSIL的惡性進展有關，尤其是這些炎癥細胞因子可能參與癌癥發展的后期階段，特別是轉移和擴散。

近年來，生物信息學已廣泛用于疾病相關的新標志物和潛在治療靶點的挖掘。相較于其他HSIL研究中使用單一芯片的方式[7]，本研究采用RRA進行了多芯片聯合分析，增大了樣本量，不僅合理地處理了4套芯片分析中實驗平臺不同的問題，還根據不相關輸入的零假設來識別基因，對結果誤差和背景信號噪聲具有較好的魯棒性，進而使得本研究獲得的疾病相關的基因更加可靠，這為HSIL的發展機制研究提供了新思路。

綜上所述，本研究鑒定了5個(CDA、LTF、CRCT1、CDKN2A、CXCL1)與HSIL相關的核心基因，可能為HSIL的篩查、診斷和預后提供可靠的分子生物標志物。其中LTF、CRCT1兩個核心基因是首次提出可能成為新的HSIL生物標志物，聯合已有的生物標志物可能有助于早期篩查診斷和預后監測。但本研究也存在一些不足，如篩選的4個芯片中，GSE63514和GSE7803樣本來自美國，GSE138080數據樣本來自荷蘭，GSE27678數據樣本來自英國，原始文獻并未給出樣本來源的具體種族。考慮到不同種群之間存在遺傳背景的差異，故本研究在推廣應用于其他人群可能會有一定限制，需要收集更多人群樣本進行驗證。