徐 瑩，褚以忞，楊大明，李 吉，張海芹，彭海霞

上海交通大學醫(yī)學院附屬同仁醫(yī)院內窺鏡室，上海 200336

結直腸癌是世界范圍內發(fā)病率第四的惡性腫瘤，其死亡率位列第二［1］。在我國，結直腸癌位于所有惡性腫瘤發(fā)病率的第三位，死亡率排名第五［2］。雖然伴隨早期腫瘤發(fā)現(xiàn)的增多及腫瘤篩查的推廣，結直腸癌患者的5 年生存率有所提高，但是仍有25%的患者確診時為Ⅳ期，并且有25%～50%的患者確診時為早期而后發(fā)展為轉移性疾?。?-6］。轉移性結腸癌患者的預后較差。一項基于美國人群的統(tǒng)計數(shù)據(jù)顯示，其5 年中位生存率僅有12.5%［4］。因此，尋找預判轉移的生物學標志物有著重要的意義。

微衛(wèi)星指的是廣泛分布于原核和真核生物基因組中的短的串聯(lián)重復序列，約占人類基因的10%，核心序列長度為1～6 bp；微衛(wèi)星不穩(wěn)定（microsatellite instability，MSI）指基因組中短串聯(lián)重復序列次數(shù)的增加或者減少［7-8］。MSI 狀態(tài)的檢測可以通過PCR 擴增特定微衛(wèi)星標志位點、定向二代測序（next-generation sequencing，NGS）或免疫組化檢測DNA 錯配修復（mismatch repair，MMR）蛋白表達來確定，并可根據(jù)PCR 擴增出的標志位點不穩(wěn)定數(shù)或定向NGS 中的不穩(wěn)定位點累積得分將MSI 分為高度微衛(wèi)星不穩(wěn)定（microsatellite instabilityhigh，MSI-H） 和低度微衛(wèi)星不穩(wěn)定（microsatellite instability-low，MSI-L）［9-11］。MSI 為結直腸癌的特殊分子表型，在結直腸癌中占10%～20%［7-8］。既往研究［12-15］顯示，MSI 結直腸癌淋巴結轉移和遠處轉移的發(fā)生率低于微衛(wèi)星穩(wěn)定（microsatellite stable，MSS）腫瘤。但是在轉移性結直腸癌中，MSI-H 患者的預后較差［16］。并且在早期結直腸癌中，MSI-H 的比例高達20%，而晚期轉移性結直腸癌中，該比例則顯著下降（4%～5%）［17-18］，這提示在MSI-H 結直腸癌中存在一部分高轉移潛能的腫瘤。

本研究擬從轉錄組層面進一步探索影響MSI-H 結直腸癌高轉移潛能的因素。從癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數(shù)據(jù)庫中篩選有無轉移的MSI-H結直腸癌患者中的差異表達基因（differentially expressed gene，DEG），選取其中的關鍵基因構建轉移預測列線圖，以幫助結直腸癌臨床治療及隨訪策略的制定。

1 材料與方法

1.1 數(shù)據(jù)收集

根據(jù)以下標準，從TCGA 數(shù)據(jù)庫中篩選符合要求的患者：①確診結直腸癌。②分子分型為MSI-H（通過NGS 的全外顯子組測序數(shù)據(jù)得到）［19］。③患者標注有轉移信息。

1.2 DEG篩選與功能富集分析

使用R 語言edgeR 軟件包［20］對收集到的轉移組及無轉移組中DEGs 進行分析。將差異閾值設定為誤報率（false detective rate，F(xiàn)DR）＜0.05、log2∣差異倍數(shù)（fold change，F(xiàn)C） ∣＞1，其中FDR 為P值的多重校正值，log2FC為基因表達變化的數(shù)值和方向。

使用DAVID 數(shù)據(jù)庫（the Database for Annotation，Visualization and Integrated Discovery）對DEGs 進行基因本體數(shù)據(jù)庫（Gene Ontology，GO）注釋與聚類，篩選條件定為P＜0.05、基因數(shù)≥5。

使用在線分析工具WEB-based GEne SeT AnaLysis Toolkit，參照京都基因和基因組百科全書（Kyoto Encyclopedia of Genes and Genomes， KEGG） 和Reactome信號通路數(shù)據(jù)庫，對DEGs進行基因集富集分析（Gene Set Enrichment Analysis，GSEA），分析DEGs 涉及的信號通路。

使用STRING 在線工具，構建上調基因（轉移組較無轉移組表達升高基因）的蛋白質互作（protein-protein interaction，PPI）網絡，通過Cytoscape 軟件［21］篩選出PPI網絡中高連接度等級前10位的樞紐基因（hub基因）。

1.3 轉移預測模型構建

選取DEGs 中調整后P值（adjustedPvalue，Padj）最小，且根據(jù)既往報道與腫瘤發(fā)生發(fā)展關聯(lián)性高的前10 個基因，通過R 語言的rms 包對DEGs 構建Logistic 回歸模型，并通過R 語言的“Boot”函數(shù)［22］，使用Bootstrap 方法，將收集的63個腫瘤樣本每次隨機抽取25個，進行15次檢驗，對構建模型進行交叉驗證。使用一致性指數(shù)（concordance index，C-index）、受試者工作特征曲線（receiver operating characteristic curve，ROC 曲線）對模型預測效能進行評價。使用R語言的rms包繪制可視化列線圖。

1.4 轉移預測模型中各基因對MSI-H 結直腸癌無進展生存期影響分析

結合收集到的TCGA 數(shù)據(jù)庫中的臨床信息，使用Log-rank檢驗并將基因中位表達值作為臨界值，采用R語言survminer 包行生存分析，分析列線圖中每個基因表達水平對MSI-H 結直腸癌無進展生存期（progression-free survival，PFS）的影響。

1.5 統(tǒng)計學方法

轉移預測模型構建使用Logistic 回歸模型。模型預測效能評價使用C-index、ROC曲線。其中C-index＞0.7則證明預測模型有可靠性；ROC 曲線的曲線下面積（area under curve，AUC）用來預測準確性，當0.5＜AUC＜1 時提示優(yōu)于隨機猜測，模型有預測價值。

2 結果

2.1 MSI-H結腸癌轉移與無轉移組間DEGs篩選

基于篩選標準，一共從TCGA 數(shù)據(jù)庫中納入63 例患者。納入患者的基本信息見表1。根據(jù)標注的轉移信息將患者分為轉移組（21例）及無轉移組（42例），轉移組為存在淋巴結轉移或/和遠端轉移，無轉移組為不存在淋巴結轉移和遠端轉移。對2 組進行轉錄組分析，共獲得245個DEGs，其中轉移組較無轉移組表達升高的有204 個，轉移組較無轉移組表達降低的有41 個（圖1）。上調及下調前10位DEG的詳細信息見表2。

表1 63例納入患者基本信息Tab 1 Basic information of 63 included patients

表2 轉移組和無轉移組間上調及下調前10位DEGsTab 2 Top 10 up-regulated and down-regulated DEGs between metastasis and non-metastasis group

圖1 MSI-H結直腸癌轉移組和無轉移組間DEGs火山圖Fig 1 Volcano plot of DEGs between metastatic group and non-metastatic group of MSI-H colorectal cancer

2.2 功能注釋與富集分析

對DEGs 進行GO 注釋及富集分析，將篩選條件定為P＜0.05，基因數(shù)≥5，在上調基因中得到8 項生物過程（biological process， BP）、 12 項細胞組分（cellular component，CC）、3 項分子功能（molecular function，MF）。上調基因的BP 主要為分泌、離子穿膜轉運、神經肽信號通路等，CC 主要為細胞外部分、質膜的組成部分、細胞外空間等，MF 主要為激素活性、氯離子通道活性和生長因子活性。在下調基因中得到1 項CC，為細胞外部分（圖2）。

圖2 GO分析上調與下調基因的BP、CC與MFFig 2 Relative BP,CC and MF of up-regulated and down-regulated genes

分別使用KEGG 和Reactome 信號通路數(shù)據(jù)庫對DEGs 進行信號通路的GSEA，并選取上調及下調基因富集前10 位的信號通路；通過對2 個數(shù)據(jù)庫富集信號通路取交集發(fā)現(xiàn)，上調基因中神經活性物質配體-受體相互作用、代謝信號通路在2 個數(shù)據(jù)庫中都得到了富集（圖3）。

圖3 DEG信號通路的GSEAFig 3 GSEA of DEG pathway

使用STRING 工具構建了上調基因PPI 網絡（圖4A），并且通過Cytoscape 軟件篩選出其中高連接度等級前10 位的hub 基因，分別為胰高血糖素（glucagon，GCG）、生長激素抑制素（somatostatin，SST）、神經降壓素（neurotensin，NTS）、 α2-HS 糖蛋白（α2-HS glycoprotein，AHSG）、載脂蛋白B （apolipoprotein B，APOB）、嗜鉻粒蛋白B（chromogranin B，CHGB）、突觸素（synaptophysin，SYP）、胰島素樣生長因子結合蛋白3（insulin like growth factor binding protein 3，IGFBP3）、精氨酸加壓素受體2 （arginine vasopressin receptor 2，AVPR2）、 分泌粒蛋白3 （secretogranin Ⅲ，SCG3）（圖4B）。

圖4 上調基因PPI網絡及等級前10位的hub基因Fig 4 PPI network of the up-regulated genes and the top 10 hub genes

2.3 MSI-H結直腸癌基因轉移預測模型的建立

將訓練集AUC=0.975，驗證集AUC=0.920，C-index=0.832（95%CI0.798～0.866）的預測模型使用R 語言的rms包繪制可視化列線圖，即MSI-H結直腸癌基因轉移預測列線圖模型（圖5）。其中，肌動蛋白8（actin like 8，ACTL8）、鳥苷酸環(huán)化酶激活劑2B （guanylate cyclase activator 2B，GUCA2B）、 L 氧化戊二酸脫氫酶（oxoglutarate dehydrogenase L，OGDHL）及視黃醇結合蛋白4（retinol binding protein 4，RBP4）對結直腸癌轉移有較大影響。

圖5 MSI-H結直腸癌轉移風險列線圖模型Fig 5 Nomogram model of MSI-H colorectal cancer metastatic risk

2.4 轉移預測模型中各基因對PFS的影響

將列線圖中各基因表達水平對MSI-H 結直腸癌PFS的影響進行生存分析（圖6），發(fā)現(xiàn)：AC078993.1和IGLJ2的表達水平與MSI-H 結直腸癌PFS 呈明顯負相關（P=0.011，P=0.005），二者均為上調基因；其他基因的表達水平對MSI-H 結直腸癌PFS 影響無統(tǒng)計學意義（均P＞0.05）。

圖6 生存分析列線圖中10個基因對MSI-H結直腸癌PFS影響Fig 6 Survival analysis of the 10 genes in Nomogram model on PFS of MSI-H colorectal cancer

3 討論

基因組不穩(wěn)定性的產生是結直腸癌發(fā)展過程中一個重要特征，而MSI 是造成基因組不穩(wěn)定性的重要途徑之一［23］。MSI 常常反映由于DNA 錯配修復缺陷（deficient mismatch repair，dMMR）導致的DNA 復制錯誤［7-8］。在Ⅱ期和Ⅲ期結直腸癌中，MSI 分別約占20%及12%，Ⅳ期中MSI 占4%～5%［17-18］。雖然在Ⅱ期的散發(fā)性結直腸癌中，MSI 為有益于預后的標志物［24］，但是在轉移性結直腸癌中dMMR-MSI-H 預后較差。這些現(xiàn)象提示，在MSIH結直腸癌中可能存在一類高轉移潛能亞型腫瘤，而目前尚無針對影響MSI-H 結直腸癌轉移轉錄組層面的研究。本研究希望從轉錄組層面探索影響MSI-H 結直腸癌轉移、產生高轉移潛能的可能因素。

首先從TCGA 數(shù)據(jù)庫中篩選出63 例MSI-H 結直腸癌患者的轉錄組數(shù)據(jù)。通過生物信息學方法分析了MSI-H結直腸癌有轉移和無轉移患者間的DEG，共獲得245 個DEGs。在轉移組上調的前10 位基因中，4 次跨膜結構域A12 （membrane spanning 4-domains A12，MS4A12） 是結腸特異性的鈣池調控鈣離子內流通道蛋白，在結腸癌細胞中MS4A12 蛋白表達降低會減弱癌細胞的增殖、運動性和趨化侵襲［25］。轉移組下調的前10 位基因中，F(xiàn)GL1編碼的蛋白為纖維蛋白原樣蛋白1（fibrinogen like 1，F(xiàn)GL1），屬于纖維蛋白原家族成員。FGL1在絲氨酸/蘇氨酸激酶11（serine/threonine kinase 11，LKB1）突變的肺腺癌中表達缺失會誘導上皮間質轉化和血管形成［26］。在肝細胞癌中，F(xiàn)GL1的表達缺失與肝細胞癌的低分化表型有相關性，且FGL1通過絲氨酸-蘇氨酸激酶依賴機制在肝細胞癌中起到抑癌作用［27］。肽基精氨酸脫亞胺酶3（peptidyl arginine deiminase 3，PADI3）編碼的蛋白為肽基精氨酸脫亞胺酶家族成員，PADI3 蛋白在結腸癌組織中低表達，并且PADI3 可通過熱休克蛋白90（heat shock protein，Hsp90）/ 細胞周期蛋白依賴性激酶調節(jié)亞基1（cyclin kinase subunit 1，CKS1）途徑或沉默信息調節(jié)因子2（silent information regulator2，Sirt2）/蛋白激酶B（protein kinase B，AKT）/p21 途徑發(fā)揮抑癌作用［28-29］。基于MS4A12、FGL1和PADI3既往在各種腫瘤中的功能報道，三者可能為MSI-H 結直腸癌轉移過程中較關鍵的基因，其具體作用機制值得進一步探索。

進一步通過GO分析，在BP及MF富集中發(fā)現(xiàn)MSI-H結直腸癌轉移組離子穿膜轉運、氯離子穿膜轉運及氯離子通道活性表現(xiàn)活躍。近年來，腫瘤組織中離子通道的改變越發(fā)受到重視。離子通道改變會導致細胞內離子穩(wěn)態(tài)失調，從而影響細胞體積調節(jié)、細胞膜電位變化、細胞的機械傳導及腫瘤微環(huán)境等，而這些改變都會對腫瘤細胞的增殖、凋亡、遷移、血管生成造成重要影響［30］。例如，在肝細胞癌中，氯離子選擇性通道ClC-3通過調節(jié)氯離子流調控細胞形態(tài)和體積從而促進細胞遷移［30］；在膠質瘤細胞突起中，氯離子通道ClC-3 與基質金屬蛋白酶2聚集，共同調節(jié)膠質瘤細胞的遷移和侵襲［31］等。本研究分析得到的基因富集結果提示，離子通道變化對MSI-H結直腸癌的轉移的影響值得關注。

而在CC 富集中發(fā)現(xiàn)，MSI-H 結直腸癌轉移組中細胞外部分、細胞外空間表現(xiàn)活躍。而腫瘤微環(huán)境是腫瘤細胞外的重要部分，主要包括血管、與腫瘤相關的成纖維細胞、免疫細胞和免疫抑制細胞、信號分子、骨髓源性細胞、細胞外基質等［32］。我們得到的在MSI-H 結直腸癌轉移組中細胞外成分的富集結果可能提示了腫瘤微環(huán)境對MSI-H結直腸癌轉移的發(fā)生和促進存在重要影響。

對DEGs信號通路的富集分析發(fā)現(xiàn)，代謝通路在上調基因中得到富集。既往研究［33］也提示細胞代謝失調為腫瘤重要標志之一，這說明在MSI-H 結直腸癌轉移過程中代謝通路改變也有著重要的意義。而既往報道神經系統(tǒng)異常，包括神經遞質、神經營養(yǎng)因子和其受體的異常，在結腸癌肝轉移中有重要作用［34］。我們在上調基因中也富集到神經活性物質配體-受體相互作用通路，說明該方向在MSI-H 結直腸癌的轉移研究中也值得關注。通過PPI網絡及Cytoscape 軟件篩選出的hub 基因中，NTS 被報道其蛋白的表達水平與結腸癌的預后呈負相關，并且NTS可以促進多種結腸癌細胞的生長［35］；結腸癌組織中SYP蛋白表達水平與腫瘤預后呈負相關［36］；IGFBP3的基因水平與結腸癌的預后呈負相關［37］，且IGFBP3 可以促進肺腺癌的腦轉移［38］。這些都提示發(fā)生轉移的MSI-H 結直腸癌相較于無轉移腫瘤，在轉錄組水平體現(xiàn)出更強的轉移潛能。

通過對篩選出的DEGs進行列線圖模型構建，得到了預測MSI-H 結直腸癌轉移的基因列線圖。該列線圖有較好的區(qū)分度和一定的預測效能。既往有研究［39］從轉錄組層面分析結直腸癌的基因標簽，也有研究［40］提出一組可以用于預測結腸癌轉移及復發(fā)的基因標簽。但是運用這類基因標簽進行評分，其方式較復雜，需要專業(yè)繁瑣的統(tǒng)計學方法計算。本研究主要針對結直腸癌中的特定分子亞型——MSI-H 結直腸癌，特異性地構建了預測MSIH 結直腸癌轉移風險的基因表達預測模型。該模型較直觀，計算MSI-H 結直腸癌發(fā)生轉移的風險預測值較方便。隨著基因測序技術的不斷發(fā)展及在臨床應用中的推廣，對于結腸癌術后行腫瘤組織轉錄組測序的患者，可應用此類預測模型評估其發(fā)生轉移的風險，以協(xié)助術后治療及隨訪方案的制定。將列線圖中各基因表達水平對MSIH 結直腸癌PFS 影響情況進行生存分析， 得到AC078993.1和IGLJ2與PFS呈負相關，且兩者均為MSI-H轉移性結腸癌中上調基因。結合兩者的表達水平及對PFS的影響，提示其可能在MSI-H 結直腸癌發(fā)生轉移的過程中較重要。但其在腫瘤中均未有相關報道，在之后進一步的研究中值得重點關注。

本列線圖模型僅使用了單個數(shù)據(jù)庫中的數(shù)據(jù)進行內部驗證，使結果有一定的局限性。而由于MSI 結直腸癌在總結直腸癌中占10%～20%［7-8］，且MSI 結直腸癌轉移發(fā)生率較MSS 結直腸癌相對低等特點［12-15］，在各數(shù)據(jù)庫及數(shù)據(jù)集中可得到的發(fā)生轉移的MSI 結直腸癌病例數(shù)相對較少，這使得之后的研究需進一步收集臨床病例和其他數(shù)據(jù)庫、數(shù)據(jù)集中的相關數(shù)據(jù)，以加強驗證及優(yōu)化該列線圖模型。并且本研究使用的數(shù)據(jù)無法區(qū)分MSI 結直腸癌中林奇綜合征患者，之后的研究也需結合其他數(shù)據(jù)庫信息，區(qū)分出林奇綜合征及錯配修復蛋白甲基化導致的MSI-H結直腸，進行更細致的探討。

綜上，本研究從轉錄組層面探索了影響MSI-H 結直腸癌產生高轉移潛能的因素。通過對TCGA 數(shù)據(jù)庫中相關數(shù)據(jù)的生物信息學分析，獲得對MSI-H 結直腸癌發(fā)生轉移有影響的潛在關鍵基因；并且構建了有一定預測效能的MSI-H 結直腸癌轉移基因預測模型。由于本研究局限于單一數(shù)據(jù)庫的數(shù)據(jù)，未來還需通過擴大數(shù)據(jù)的收集范圍、收集臨床樣本及生物學實驗等多個方面驗證及完善本研究的結果。