TCGA數據庫4種消化系統癌癥分子特征的比較①

劉 菲 袁 銳 陳世龍 王永翠

（中國科學院西北高原生物研究所青海省作物分子育種重點實驗室，西寧 810008）

消化系統癌癥（digestive cancers，DSCs）包括咽喉癌、食管癌（esophageal carcinoma，ESCA）、胃癌（stomach adenocarcinoma，STAD）等消化道癌癥和胰腺癌、膽囊癌、肝癌等消化道附屬器官癌癥，是存在于世界范圍內的健康問題［1-2］。其中，ESCA在癌癥病死率中排名第六，直結腸癌被認為是DSCs中發病率最高的一種癌癥，每年全球新增病例約100萬例，死亡55萬余人［3-5］。在過去的20年中，結腸癌（colon adenocarcinoma，COAD）、直腸癌（rectum adenocarcinoma，READ）、STAD均為全球發病率最高的10種癌癥類型之一。2018年，美國記錄了約319 160例DSCs新病例和160 820例與DSCs相關的病死病例，并且DSCs患者的五年生存率均低于20%［6-7］。隨著醫療水平的不斷提高，美國癌癥協會（ACS）、疾病控制與預防中心（CDC）、美國國家癌癥研究所（NCI）和北美中央癌癥登記處協會（NAACCR）在協作提供的美國癌癥發生和發展趨勢的年度報告中提出，2011年至2015年間，DSCs的病死率均呈下降趨勢，但其整體生存率依然很低［8］。

癌癥治療最終目的是根除癌細胞，同時盡量減少治療過程中對正常組織的潛在損傷［9］。傳統的腫瘤治療方法是依靠外界力量殺死腫瘤細胞，通常采用的手術治療方案會對正常組織產生損傷，并且預后效果不好，容易復發，即便患者在接受手術切除前后進行放射治療或者化學藥物治療，預后仍然很差［10-12］。由于常規治療方案對轉移或復發的DSCs治療效果很差，尋找其他新型治療方法必然成為癌癥治療的新趨勢［13-16］。近年來，靶向治療和免疫治療等腫瘤生物治療方法已經是相對比較前沿的癌癥治療手段［16］。依據理論，人體正常細胞和腫瘤細胞的分子特征存在很大不同，這些差異可以作為靶點，從而針對已經明確的分子靶點進行靶向治療［17］。已有大量研究通過機器學習或差異分析的方法揭示一些基因的異常表達或甲基化會影響癌癥的發展，并得到一些潛在的癌癥預后基因及其通路［18-19］。比如DAVID等［20］同時采用差異分析和機器學習的方法鑒定出多種與腫瘤發展相關的分泌途徑有關基因，SAGHALEYNI等［21］則將隨機森林的方法用于癌癥基因組圖譜（cancer genome atlas，TCGA）、基因型組織表達數據（genotype-tissue expression，GTEx）、法國肉瘤組織這3種基因表達數據庫中，用以識別新型診斷標志物，并通過實驗對其結果進行驗證。

盡管之前已有大量的工作對癌癥相關的分子數據進行了廣泛的研究，但大多數研究都集中于單一分子特征或單個癌癥類型之間的關聯上。因此，一些學者開始基于癌癥集群進行綜合分析（泛癌分析），以期發現更系統、更可靠的生物標志物，從而提高預后能力。如HOADLEY［22］等將12種癌癥類型進行多平臺綜合分析，基于基因突變、DNA拷貝數、DNA甲基化、mRNA多組學數據，根據分子特征綜合比較12種癌癥類型的異同，從而為根據分子亞型對癌癥進行分類奠定了基礎。此外，BASS等［23］基于體細胞拷貝數分析、外顯子分析、DNA甲基化分析、mRNA測序、miRNA測序等方法對295個原發性胃腺癌進行綜合分子評估并創建決策樹，定義了STAD的4種主要基因組亞型，表現出明顯突出的基因組特征，為不同STAD患者相應靶向藥物的開發提供指導。不論是在全國還是在全球發布的癌癥死亡順位分布中，人類高發病率和高病死率的惡性腫瘤中均包括READ、肝癌、STAD及ESCA 4種癌癥［24］。但目前采用泛癌分析的方法綜合分析不同DSCs的多種分子數據的研究依然很少，4種DSCs依然存在一些臨床信息及分子特征差異有待于進一步挖掘［25］。

TCGA數據庫為癌癥研究提供了豐富的資源［26］。利用TCGA數據庫的20 000多種原發癌的分子特征和33種癌癥類型的正常樣本，結合臨床數據挖掘不同DSCs的分子特征并加以比較分析，對DSCs的預后具有重要研究價值。因此，本文通過TCGA分析DSCs（ESCA、STAD、READ、COAD）分子數據和臨床信息，從而挖掘DSCs分子特征的異同以及臨床基本特征的關系。其中，分子數據包括基因突變數據、基因拷貝數變異數據、基因表達數據、基因甲基化數據，所采用的分析方法主要有差異分析、富集分析、生存分析以及相關性分析。本文關于4種DSCs類型異同點的提出可以為深入理解DSCs的分子特征及相應分子靶向治療方案的設計提供有效依據。

1 資料與方法

1.1 資料 通過TCGAbiolinks包從TCGA公共數據 庫 官 網 接 口 獲 得DSCs（ESCA、STAD、READ、COAD）分子和樣本臨床數據，分子數據包括基因突變數據、拷貝數變異數據、基因表達數據以及甲基化數據，樣本分為癌癥樣本和正常樣本，正常樣本即癌旁組織樣本。此外，在TCGA數據庫中4種癌癥類型患者的臨床資料，包括患者的年齡、性別、生存時間、生存狀態、癌癥階段分期、發病位置、存活時間。刪除生存狀態缺失和生存時間不明確的患者數據，臨床信息數據最終包括ESCA患者信息185例、STAD患 者434例、READ患者169例、COAD患者458例。

1.2 方法

1.2.1 單一數據類型分析 對每種癌癥類型的不同分子數據分別進行數據處理。將基因突變數據采用MaxMIF算法進行篩選，挑選每種癌癥類型中30個最大的突變影響函數值（MIF）所對應的基因作為癌癥驅動基因［27］。拷貝數變異數據則通過文獻報道的顯著基因和由COSMIC數據庫公布的致癌基因進行篩選，從而分析拷貝數變異情況［16，23，28］。

甲基化數據與基因表達數據做相似處理，對采用log2標準化處理后的基因表達數據通過R Limma包做數據差異分析，為得到包含顯著的癌癥驅動基因在內的差異表達基因（differential expression gene，DEGs），對數據差異分析采用不同的標準［ESCA：|logFC|＞1.2，矯正P值（FDR）＜0.05；STAD：|logFC|＞1.2，FDR＜0.001；READ：|logFC|＞1.5，FDR＜0.001；COAD：|logFC|＞1.5，FDR＜0.001］，而后將差異分析篩選出的DEGs通過DAVID數據庫（https：//david.ncifcrf.gov）進行GO和KEGG通路富集分析，-log10P值＞2差異具有統計學意義。其中GO分析包括生物學過程（biological process，BP）、細胞組成（cell composition，CC）和分子功能（molecular function，MF）3種注釋。對完成在編碼區內篩選的甲基化位點數據通過R Limma包做數據差異分析，篩選|logFC|＞0.5，FDR＜0.001的位點基因為差異甲基化基因（differential methylation gene，DMGs），為得到其GO和KEGG通路，通過DAVID數據庫對DMGs進行GO和KEGG通路富集分析，GO分析同樣包括生物學過程、細胞組成和分子功能3種注釋。

1.2.2 不同數據類型聯合分析 將癌癥驅動基因的突變數據和基因表達數據分別結合臨床信息做生存分析，分析預后基因的預后價值。根據基因在癌癥患者中是否發生突變，將癌癥患者劃分為突變組和未突變組，采用Kaplan-Meier方法，分析癌癥驅動基因的突變與4種癌癥類型發病之間的關系，從而找到預后效果良好的癌癥驅動基因，在對4種癌癥類型的癌癥驅動基因表達數據結合臨床信息生存分析中，將劃分組別的標準改為根據基因在癌癥患者中表達水平的上下四分位數劃分，癌癥患者被劃分為高表達組和低表達組，分析癌癥驅動基因的表達與4種癌癥類型發病之間的關系，并采用Survival軟件包對結果進行圖形可視化，找到具有預后效果的癌癥驅動基因。

將突變數據分別與基因表達數據、甲基化數據做聯合分析，得到突變與基因表達、DNA甲基化之間的關系。引用皮爾遜相關系數法，利用R軟件“cor”函數計算不同癌癥類型中癌癥驅動基因的突變率，分別與癌癥驅動基因的最大表達量、癌癥驅動基因甲基化數據最大β值的相關性，并采用ggplot2軟件包對結果進行圖形可視化。

1.2.3 驗證 從TCGA官方發布的文獻中找到通過Mutsig計算得到的顯著基因［16，23，28］，將已報道的顯著基因與臨床數據結合做生存分析，根據基因在患者中是否發生突變，將癌癥患者劃分為突變組和未突變組，采用KM方法，分析顯著基因的突變與4種癌癥類型發病之間的關系，從而找到預后效果良好的顯著基因。并將已報道的顯著基因突變數據分別與基因表達數據以及甲基化數據利用R軟件“cor”函數做相關性分析。將顯著基因的一系列分析結果與癌癥驅動基因的分析結果做比較。

2 結果

2.1 單一數據類型分析結果 在對ESCA突變數據的篩選中，采用MaxMIF算法計算分別得到4種癌癥類型的30個癌癥驅動基因，并得到癌癥驅動基因在癌癥的突變數據中的大體分布（圖1）。利用t檢驗法計算不同癌癥類型中每種癌癥驅動基因的突變頻率，判斷每兩種癌癥類型之間突變頻率的差異顯著性，經計算得出，兩兩癌癥類型之間P值均＜0.001，具有極顯著性差異。另外，對各種癌癥類型的癌癥驅動基因分別做預后分析，得到預后效果顯著的癌癥驅動基因，其中，ESCA中通過突變數據得到具有預后效果癌癥驅動基因為ISG20、IL1R1，SATD癌癥驅動基因中具有預后效果的癌癥驅動基因為RAB13、SELL，READ癌癥驅動基因中預后效果良好的癌癥驅動基因為C5，COAD癌癥驅動基因中預后效果良好的癌癥驅動基因為CD226（圖2）。

圖1 癌癥驅動基因突變信息分布及樣本信息分布Fig.1 Distribution of sample informations and mutation informations based on cancer driver genes

圖2 4種癌癥類型的癌癥驅動基因預后分析Fig.2 Prognosis analysis of cancer driver genes from four types of cancer

通過癌癥驅動基因篩選拷貝數變異數據，觀察癌癥驅動基因的拷貝數變異情況。將樣本分為腫瘤組織樣本與正常組織樣本，再按基因拷貝數的正負將癌癥驅動基因分別分為拷貝數增加和減少2種，以此得到驅動基因拷貝數的4種類型：腫瘤樣本中增加、腫瘤樣本中減少、正常樣本中增加、正常樣本中減少。在4種癌癥類型中，癌癥組織樣本和正常組織樣本的個數占總樣本量的比值均接近于50%，而每種基因在癌癥組織樣本和正常組織樣本中的增加和減少的樣本分布并不均勻。針對文獻報道的顯著基因以及COSMIC數據庫公布的常見致癌基因（BRAF、CDKN2A、CTNNB1、EGFR、ERBB2、KIT、KRAS、MYC、PIK3CA、RB1、RET、SMAD4、TP53）統計拷貝數變異的情況（圖3）［16，23，28］。可以看出，在不同癌癥類型中關鍵基因的拷貝數變化在癌癥組織樣本和正常組織樣本中的分布并不均勻。

圖3 癌癥驅動基因的拷貝數變異數據分布情況Fig.3 Distribution of copy number variation data of cancer driver genes

癌癥基因表達數據和癌癥甲基化基因數據分別通過差異分析篩選，在癌癥基因表達數據中篩選到829個ESCA DEGs、394個STAD DEGs、606個READ DEGs、591個COAD DEGs。在癌癥甲基化數據中篩選到88個ESCA DMGs、36個STAD DMGs、515個READ DMGs、127個COAD DMGs。對所得到的不同癌癥類型DEGs通過DAVID進行KEGG功能富集分析，進一步了解DEGs的富集通路。KEGG富集通路顯示，與基因表達相關的差異基因顯著富集在神經活性配體與受體相互作用方面（圖4）。另外，通過上下四分位數將表達數據分為高表達組和低表達組，并以此為分組做預后分析，得到ESCA中具有預后效果的基因有29個，STAD中具有預后效果的基因有70個，READ中具有預后效果的基因有71個，COAD中具有預后效果的基因有50個。

圖4 DEGs KEGG功能富集分析Fig.4 KEGG functional enrichment analysis of DEGs

2.2 利用不同數據類型聯合分析的結果 將突變數據分別與基因表達數據做相關性分析，以不同的癌癥類型中癌癥驅動基因的突變率分別與癌癥驅動基因的最大表達量以及甲基化最大Beta值做相關性分析，其中ESCA基因突變率與基因表達量、甲基化Beta值的相關系數r為0.110、0.091，STAD基因突變率與基因表達量、與甲基化Beta值的相關系數r為-0.062、0.057，READ基因突變率與基因表達量、甲基化Beta值的相關系數r為0.200、0.110，COAD基因突變率與基因表達量、甲基化Beta值的相關系數r為0.065、0.018。相關系數的絕對值均＜0.300。

通過突變數據做預后分析從TCGA官方發布的文獻中找到的顯著基因，并將顯著基因的突變率分別與基因表達數據、甲基化數據做相關性分析，以此驗證MaxMIF算法找到的癌癥驅動基因預后率與其他分子數據相關性分析的準確性。ESCA中顯著基因32個，STAD中顯著基因23個，READ中顯著基因28個，COAD中顯著基因個數28個，具有預后效果的顯著基因所占的比例分別為9.375%、13.043%、0、7.143%，而具有預后效果的癌癥驅動基因所占的比例分別為3.333%、6.667%、3.333%、6.667%。預后基因比率差別不大，但通過癌癥驅動基因分析得到的潛在預后基因個數穩定一些。顯著基因的突變率與基因表達數據相關系數分別為0.008、0.172、-0.031、-0.130，顯著基因的突變率與甲基化數據相關系數分別為-0.232、-0.187、0.081、0.220。同樣的，顯著基因與其他分子數據的相關系數也均＜0.300。

3 討論

DSCs是常見的惡性腫瘤，由于早期無明顯不良狀況，診斷率低，發現時多瀕臨或處于中晚期，診斷費用高、治療療效差，并伴有易轉移、易復發的危險，因此發病率、病死率高，另外，我國DSCs的發病人群也日趨年輕化［28-31］。在最近的10年，隨著高通量測序技術的出現、測序成本的大幅度降低以及基因組學、轉錄組學、蛋白組學、表觀組學等技術的進步，醫療大數據得以迅速發展。大規模的研究已經繪制了基因組調控網絡，并極大地增進了對包括癌癥發生和癌癥轉移在內的多種生物學過程的分子機制的理解［32］。

已有證據表明，有些基因在多種癌癥組織中起著相似作用。例如，基因TP53通過影響NOX4表達與多種癌癥的進展建立關系，包括ESCA、甲狀腺癌、膀胱尿道癌等，并可調節胰腺癌細胞對化療藥物、信號轉導抑制劑以及營養素的敏感性［33-34］。此外，癌癥分析還發現不同癌癥之間類似的分子特征變化，例如MCM8基因在STAD、READ、肝癌、肺癌、宮頸癌、卵巢癌、淋巴瘤、頭頸部和中樞神經系統腫瘤中過度表達，ESM1基因在12種癌癥類型中過度表達并與宮頸鱗癌和腺癌、ESCA、腎透明細胞癌和腎乳頭狀細胞癌4種癌癥患者的總生存期顯著相關［35-37］。

本研究通過檢索TCGA數據庫發現，每種癌癥類型都具有特有的潛在預后基因，如ESCA中的ISG20、IL1R1，STAD中的RAB13、SELL，READ中的C5、GPC6、EDNRB，以及COAD中的CD226。ISG20基因是由GONGORA等［38］作為早幼粒細胞白血病核體（PML-NB）相關蛋白首先引入的。ISG20的上調表達可以促進腫瘤轉移和血管生成，并與肝細胞癌患者的無疾病復發存活期較短呈顯著相關，另外，ISG20可以促進局部腫瘤免疫，可能成為人類神經膠質瘤的潛在治療靶點［39-40］。由IL1R1編碼的具有中繼信號作用的IL-1Rα與IL-1結合，后IL-1Rα與IL-1受體輔助蛋白形成復合物，誘導了多種通路的激活，并且這些通路都是腸道腫瘤發生的關鍵［41-43］。RAB13通過多種機制影響細胞遷移，在集體遷移潛在部位也存在豐富的RAB13 RNA，從而為針對集體癌細胞入侵提供有效的打擊機會［44］。SELL基因翻譯產生L-選擇蛋白，是一種細胞黏附分子，在調節免疫細胞運輸中起重要作用，并且已經被發現其表達在幾種人類癌癥中被上調，如SELL上調表達與炎癥性微環境以及乳腺癌的良好預后有關［45］。補體系統是先天免疫系統的重要組成部分，在適應性免疫的激活和調節中起著關鍵作用，而C5是末端補體途徑中補體級聯的一個組分，補體激活導致C5裂解成小C5a分子和大C5b分子，C5a分子具有免疫刺激作用，而C5b分子是膜攻擊復合物（MAC）的初始成分［46-47］。CD226是在各種免疫細胞上表達的激活受體，充當共刺激物，增強T細胞或自然殺傷細胞的活化，并在先天性和適應性免疫調節網絡中顯示出重要意義［48-49］。可見不同癌癥類型中的預后效果良好的基因大多可能與腫瘤免疫方向有關。

對DEGs進行KEGG富集通路分析，得到在神經活性配體與受體相互作用方面顯著富集的結果，表明DEGs可能通過神經活性配體-受體相互作用途徑在4種DSCs中發揮作用。先前有研究表明，該途徑似乎在同為DSCs的肝癌進展中起關鍵作用，在肝癌的每個階段都能觀察到神經活性的配體-受體相互作用［50-51］。另一項研究說明，神經活性配體-受體相互作用途徑參與了腫瘤微環境和細胞間的通訊，其失調可能促進了癌癥的發展［52］。

綜上所述，本研究報告了ESCA、STAD、READ、COAD 4種DSCs在多組學層面的綜合分析。4種癌癥類型具有相似的分子特征，包括：4種癌癥類型的突變頻率與基因表達、甲基化水平無明顯相關性；MaxMIF軟件獲得的癌癥驅動基因中均僅有少數幾個基因具有良好的預后特征，如ESCA癌癥驅動基因中的ISG20、IL1R1，STAD癌癥驅動基因中的RAB13、SELL，READ癌癥驅動基因中的C5，以及COAD癌癥驅動基因中的CD226；拷貝數變化在癌癥組織樣本和正常組織樣本中的分布并不均勻；DEGs中均僅有少部分基因具有良好的預后特征，且KEGG通路都富集在“神經活性的配體-受體相互作用”這條通路上。同時，獲得4種癌癥類型不同的分子特征，即4種癌癥類型的癌癥驅動基因、預后良好的DEGs、DMGs均不同。總之，以上結果的提出有助于了解整個DSCs的異同，并有望為深入理解DSCs的分子特征及相應的分子靶向治療方案的設計提供有效的依據。