駱澤民，方建惠，韋良宏，陳海東，易廷莊

(1. 廣西欽州市第一人民醫(yī)院消化內(nèi)科，廣西 欽州 535000；2. 廣西欽州市第一人民醫(yī)院腫瘤內(nèi)科，廣西 欽州 535000；3. 右江民族醫(yī)學院附屬醫(yī)院，廣西 百色 533000)

結(jié)直腸癌(colorectal cancer，CRC)是消化到最常見的惡性腫瘤之一，根據(jù)最新的流行病學調(diào)查的研究發(fā)現(xiàn)，每年在全球范圍內(nèi)的新發(fā)病例超過了100萬例，更重要的是超過一半的患者確診時已經(jīng)到了晚期[1]。淋巴結(jié)轉(zhuǎn)移是晚期CRC最常見的轉(zhuǎn)移方式[2]，同樣也是CRC患者預后的決定性因素之一[3]。具體說來，對于沒有淋巴結(jié)轉(zhuǎn)移的患者，術(shù)后的5年生存率甚至超過了80%，而有淋巴結(jié)轉(zhuǎn)移的患者即便沒有遠處轉(zhuǎn)移的征象，術(shù)后的5年生存率下降到了約55%[3-4]。然而，目前對晚期CRC的淋巴結(jié)轉(zhuǎn)移的關(guān)鍵機制尚不明確，這制約了進一步對CRC淋巴結(jié)轉(zhuǎn)移針對性的臨床進展。為此，本研究通過生物信息學方法對可切除的遠處轉(zhuǎn)移性結(jié)直腸癌(Ⅳ期)淋巴結(jié)狀態(tài)進行分析，并對原發(fā)病灶的基因表達譜進行運算，通過一系列的方法探索影響Ⅳ期CRC患者淋巴結(jié)轉(zhuǎn)移相關(guān)的核心基因，為進一步對轉(zhuǎn)移相關(guān)的潛在機制的探索提供線索。

1 方法

1.1 表達譜芯片下載及分析 本研究通過對GEO數(shù)據(jù)庫(https：//www.ncbi.nlm.nih.gov/geo/)進行檢索，檢索關(guān)鍵詞包括“colorectal cancer”，“CRC”，“l(fā)ymph node”，“metastasis”等，最終選擇GSE63596的表達譜芯片進行后續(xù)研究。GSE63596數(shù)據(jù)集采用了GPL17077平臺文件，納入了15例可切除的Ⅳ期CRC患者，其中7例術(shù)后病檢提示淋巴結(jié)陰性，8例為淋巴結(jié)陽性。男性10例、女性5例；年齡范圍在45～80歲之間。采用Agilent微陣列技術(shù)對所有15例患者的mRNA譜進行分析，差異分析采用LIMMA包進行運算，統(tǒng)計學過濾條件為：|Log2FC|>2，P<0.05。

1.2 淋巴結(jié)轉(zhuǎn)移風險模型 本研究通過對術(shù)后淋巴結(jié)的狀態(tài)以及原發(fā)病灶的表達譜數(shù)據(jù)進行分析，采用加權(quán)基因共表達網(wǎng)絡(luò)分析(weighted gene co-expression network analysis，WGCNA)方法篩選出與CRC淋巴結(jié)轉(zhuǎn)移相關(guān)的靶基因并構(gòu)建預測模型。首先分別計算了每個基因的標準差(standard deviation，SD)，剔除了SD最小的前25%的基因，去除了離群的基因和樣本，隨后通過主要連接關(guān)系和Pearson相關(guān)矩陣建立兩個基因之間的相關(guān)矩陣。并通過對網(wǎng)絡(luò)拓撲結(jié)構(gòu)的分析，確定軟閾值大小和相應的平均連通性。 拓撲重疊矩陣(topological overlap matrix，TOM)用于度量一個基因的網(wǎng)絡(luò)連通性，為了將表達譜相似的基因分類到基因模塊中，根據(jù)基于TOM的差異測度進行平均連鎖層次聚類。為了進一步分析模塊，計算模塊特征基因(MEs)的差異性，通過Pearson檢驗來評估MEs與淋巴結(jié)轉(zhuǎn)移的相關(guān)性，以確定相關(guān)模塊。選擇與轉(zhuǎn)移高度相關(guān)的模塊作為淋巴結(jié)轉(zhuǎn)移模塊進一步進行分析。

1.3 蛋白互作網(wǎng)絡(luò)構(gòu)建及LASSO回歸分析 選擇與淋巴結(jié)轉(zhuǎn)移風險最高模塊進行分析。利用STRING數(shù)據(jù)庫(https：//string-db.org/)進行蛋白-蛋白的互作(PPI)分析，刪除沒有互作關(guān)系的蛋白。同時利用R軟件的LASSO回歸分析包篩選轉(zhuǎn)移模塊中基因。logrank用于檢驗KM生存分析比較兩組之間的生存差異，進行了timeROC分析以比較預測模型的準確性和風險評分。LASSO回歸算法進行特征選擇，采用10倍交叉驗證，以上分析采用R軟件包glmnet。對于Kaplan-Meier曲線，P值和95%置信區(qū)間(CI)的危險比(HR)通過logrank檢驗和單變量Cox比例危險回歸得出。

2 結(jié)果

2.1 差異基因分析 差異基因分析采用R軟件LIMMA包進行篩選，基于|Log2FC|>2，P<0.05。為閾值篩選GSE63596數(shù)據(jù)集中淋巴結(jié)陽性與淋巴結(jié)陰性的原發(fā)性腫瘤的差異基因。根據(jù)相應平臺注釋信息，統(tǒng)計將探針I(yè)D轉(zhuǎn)換為gene symbol，如涉及多個探針對于同一個gene symbol時，取多個探針對應的平均值。分析之前采用normalize.quantiles函數(shù)進行數(shù)據(jù)標準化處理，以及主成因分析(principle component analysis，PCA)的預處理。針對數(shù)據(jù)預處理結(jié)果，通過箱線圖進行評估數(shù)據(jù)標準化情況(見圖1A)；數(shù)據(jù)批次效應情況通過對比批次去除前后可視化PCA圖進行評估(見圖1B、圖1C)。差異分析的結(jié)果提示有501個基因在淋巴結(jié)陽性的腫瘤中表達上調(diào)，另外有102個基因表達下調(diào)。差異表達基因(differentially expressed genes，DEGs)的具體分布采用火山圖(見圖1D)以及聚類分析的熱圖(見圖1E)進行可視化展現(xiàn)。

為探索淋巴結(jié)轉(zhuǎn)移的相關(guān)分子機制和潛在的信號通路，將DEGs進行GO功能富集以及KEGG信號通路富集分析。結(jié)果提示在淋巴結(jié)陽性的腫瘤組織中cAMP信號通路被顯著激活，而對白介素1的反應受到了顯著抑制(見圖2)。然而DEGs的富集分析雖然在一定程度上提示了淋巴結(jié)轉(zhuǎn)移的相關(guān)機制，但是這其中存在很多不相關(guān)因素的干擾。因此，還需要進一步分析并篩選出淋巴結(jié)轉(zhuǎn)移的核心功能分子。

A： 數(shù)據(jù)標準化后箱線圖。B、C： PCA結(jié)果圖；D：差異基因火山圖；E：差異基因表達熱圖

圖2 差異表達基因富集分析

2.2 基因模塊與臨床性狀關(guān)聯(lián)分析 本研究通過構(gòu)建WGCNA模型對經(jīng)過處理后的數(shù)據(jù)進行分析，確定了軟閾值(β)；β 值的一般范圍在1～20之間時，本研究中 β=20 (見圖3A)以及相應的平均連接度= 0.12 (見圖3B)。

圖3 WGCNA模型的軟閾值(β)以及平均連接度

通過R語言軟件中cutreeDynamic和moduleEigengenes函數(shù)繪制聚類圖，去除相似模塊的影響，共得到15個模塊(見圖4A)，每種模塊間的相關(guān)性(見圖4B)。需要分析的主要變量參數(shù)為淋巴結(jié)轉(zhuǎn)移的狀態(tài)，方塊中分別顯示了每個模塊與其的相關(guān)性系數(shù)及其P值(見圖4C)，結(jié)果發(fā)現(xiàn)淋巴結(jié)陽性與drakgrey模塊顯著相關(guān)(cor=0.64，P<0.01)。為了進一步分析模塊，計算模塊特征基因(MEs)的差異性，通過Pearson檢驗來評估MEs與淋巴結(jié)轉(zhuǎn)移的相關(guān)性，以確定相關(guān)模塊，結(jié)果也證實了drakgrey模塊的基因與淋巴結(jié)轉(zhuǎn)移的強相關(guān)性(r=0.53，P<0.01)，見表1、圖5A。將drakgrey模塊中的DEGs提取出來，構(gòu)建PPI網(wǎng)絡(luò)后刷選出8個Hub基因(見圖5B)。

圖4 WGCNA模型的核心模塊篩選

表1 drakgrey模塊與淋巴結(jié)轉(zhuǎn)移的相關(guān)的基因

圖5 淋巴結(jié)轉(zhuǎn)移的核心模塊與核心基因

2.3 LASSO 回歸分析 為確定上述方法篩選出的Hub基因?qū)︻A后的影響，基于TCGA數(shù)據(jù)庫，下載Ⅳ期結(jié)直腸癌的表達譜數(shù)據(jù)和臨床相關(guān)數(shù)據(jù)作為驗證級。在本研究中，使用R軟件包glmnet，整合生存時間、生存狀態(tài)和基因表達數(shù)據(jù)，利用LASSO回歸分析，以獲得最優(yōu)模型。本研究設(shè)置Lambda值為0.0407，獲得基因4個(DGAT1、GDPD3、PDZK1IP1、PLA2G10)，見圖6，構(gòu)建的模型公式為Riskscore=(0.1454)*DGAT1+(-0.0684)*GDPD3+(0.1404)*PDZK1IP1+(-0.1724)*PLA2G10。

圖6 LASSO 回歸分析篩選Lambda值以及構(gòu)建Riskscore模型

基于LASSO 回歸分析的構(gòu)建Riskscore模型評分，分為高風險組和低風險組，進一步使用“pheatmap”程序包繪制風險曲線圖和生存狀態(tài)圖?？梢妳?shù)基因在高危風險組的預后明顯短于低風險組(HR=2.63，95%CI：1.857～3.1112，P=0.036)。ROC分析用于模型的構(gòu)建過程中對結(jié)果預測的性能的一個評估。本研究中的預后模型分別繪制了1、3、5年OS的ROC曲線，結(jié)果提示該模型具有較強的預后判斷價值(見圖7)。

圖7 Riskscore模型的預后預測價值

3 討論

CRC的淋巴轉(zhuǎn)移的發(fā)生機制目前仍舊存在爭議，隨著高通量技術(shù)的發(fā)展以及生物信息方法的應用，使得人們對病因以及機制的研究有了新的方法。本研究通過WGCNA以及LASSO 回歸分析最終確定了4個與淋巴轉(zhuǎn)移最為相關(guān)的Hub基因(DGAT1、GDPD3、PDZK1IP1、PLA2G10)，并基于其表達水平構(gòu)建了Riskscore評分的預測模型。

DGAT1是一種表達于內(nèi)質(zhì)網(wǎng)的多膜蛋白。它可以通過將酰基輔酶A轉(zhuǎn)化為甘油三酯來合成甘油三酯，甘油三酯在脂質(zhì)合成中起著關(guān)鍵作用[5-6]。DGAT1是MBOAT家族的成員，該家族的其他成員同樣在脂質(zhì)代謝、信號轉(zhuǎn)導和激素治療中發(fā)揮重要作用[7]。有研究認為DGAT1參與腫瘤細胞中IL-8表達的負調(diào)控[8]。IL-8可通過Wnt/β-catenin途徑介導的EMT促進腫瘤的遷移[9]。 GDPD3在維持CML干細胞靜止和TKI抵抗中起重要作用[10-11]。另外，在肌層浸潤性膀胱尿路上皮癌患者中，GDPD3的表達水平還能與新輔助化療的反應具有相關(guān)性[12]。但是尚未有足夠證據(jù)表明，GDPD3可以直接參與了腫瘤的轉(zhuǎn)移過程。PDZK1IP1是一種分子量比較小的蛋白(～17 kDa)，在多種癌癥中上調(diào)，雖然PDZK1IP1沒有酶活性或轉(zhuǎn)錄活性，但它通過調(diào)節(jié)多種細胞信號通路發(fā)揮其所描述的作用。PDZK1IP1過表達激活Notch通路，從而使該通路更為活躍[13-14]。此外，PDZK1IP1的過表達可以通過增加活性氧(ROS)減少NFκB活化和細胞自噬[15]。ROS活性的增加可以提高基因突變率和基因組的不穩(wěn)定性[16-17]。有研究報道[18]，PLA2G10的過度表達促進了SK-LMS-1細胞的生長和G1/S轉(zhuǎn)換，反之，PLA2G10的敲除減緩了細胞生長并導致G1期阻滯。而這種現(xiàn)象在CRC中同樣被觀察到，且被認為與炎癥的發(fā)生密切相關(guān)[19-20]。

在本研究中，篩選取4個基因構(gòu)建了轉(zhuǎn)移相關(guān)的預后模型，并繪制風險曲線及熱圖，希望從基因角度對預測CRC患者的預后，以期為CRC的發(fā)生、發(fā)展及轉(zhuǎn)歸的生物學過程研究提供依據(jù)。但不可否認的是，本研究結(jié)果還需進一步在臨床樣本中獲得更為可靠的驗證，以及通過體內(nèi)/外實驗驗證上述基因的臨床和生物學意義。