基于TCGA SpliceSeq數據庫可變剪接事件構建結腸癌患者預后風險模型

雷 鳴，郭萌月，王若穎，倪小梅，石 瓊

（云南省腫瘤醫院 昆明醫科大學第三附屬醫院 云南省癌癥中心，云南 昆明 650118）

結腸癌（colon adenocarcinoma，COAD）是常見的惡性腫瘤之一，發生率在消化道腫瘤中僅次于胃癌和食管癌[1]。腫瘤中存在大量異常mRNA剪接方式，二代高通量測序技術能有效獲取腫瘤中mRNA的表達水平和序列信息，從而發現樣本中大量的新突變、可變剪接和基因融合。可變剪接指同一個基因轉錄形成的mRNA前體通過不同的剪切和拼接方式產生不同的成熟mRNA的過程，所獲得的mRNA被稱為轉錄本，再進行翻譯可以獲得不同的蛋白質異構體[2]。可變剪接是調節基因表達和產生蛋白質多樣性的重要機制，也是真核生物基因和蛋白質數量差異較大的重要原因[3]，在細胞生長、分化、分裂、凋亡過程中具有重要作用，其異常與腫瘤細胞的凋亡、遷移、侵襲、耐藥等密切相關[4]。癌癥基因組計劃（the Cancer Genome Atlas，TCGA）SpliceSeq數據庫提供的可變剪接事件主要有7種類型：可變受體位點（alternate acceptor site，AA）、可變供體位點（alternate donor site，AD）、可變啟動子（alternate pormoter，AP）、可變終止子（alternate terminator，AT）、外顯子跳躍（exon skip，ES）、外顯子互斥（mutually exclusive exons，ME）、內含子保留（retained intron，RI）[5]。本研究基于TCGA數據庫中大樣本COAD的RNA-Seq數據及患者臨床信息，聯合TCGA SpliceSeq數據庫中的剪接事件數據，構建高效且可靠的COAD預后風險模型，并從SpliceAid 2數據庫中下載剪接因子（splicing factor，SF）數據[6]，構建其與可變剪接事件之間的調控網絡，為預測患者預后奠定基礎。

1 材料和方法

1.1 研究對象

從TCGA數據庫（https：//cancergenome.nih.gov/）中獲取COAD患者腫瘤組織（398例）及癌旁組織（39例）RNA-seq Levels 3數據和對應臨床信息，數據使用每百萬轉錄本映射的reads數（Transcripts per million reads，TPM）進行標準化處理。COAD的剪接事件從TCGA SpliceSeq數據庫（http//bioinformatics.mdanderson.org/TCGASpliceSeq/）下載，并獲得剪接百分比（percent spliced-in index value，PSI），將PSI≥75%的可變剪接事件納入到本研究中，SF數據從SpliceAid2數據庫（http：//193.206.120.249/splicing_tissue.html）中下載，共得到403個SF。

1.2 構建預后風險模型

COAD患者臨床信息排除標準：（1）患者生存時間缺失或臨床信息＜90 d；（2）缺少病理診斷的臨床信息；（3）沒有臨床信息。用單因素Cox回歸分析進行初步變量篩選，然后將得到的相關數據納入到Lasso回歸，對數據進行特征選取和降維處理，調取R語言“glmnet”包對數據進行100次10折交叉驗證，選定最小λ值為最佳λ參數值。篩選得到特征變量后，再經多因素Cox回歸分析構建出預后風險模型（Risk score）為：

公式中Expi為各可變剪接事件在樣本中的PSI值，即可能影響生存時間的有關因素；βi為多因素Cox回歸系數。根據公式計算出每例患者的預后風險值，由于基因表達量呈偏態分布，故將計算得到的風險值中位數設為最佳臨界值，將患者分為高風險組和低風險組，并以60個月為研究終點。

1.3 預后相關可變剪接事件（prognostic related alternative splicing events，PASE）功能富集分析

將單因素Cox回歸分析篩選出的PASE導入到Cytosacpe Reactome FI插件用于蛋白互作分析。由R語言“clusterProfiler”包進行基因本體（Gene Ontology，GO）富集[包括生物過程（biological processes，BP）、分子功能（molecular function，MF）、細胞構成（cellular component，CC）]和KEGG通路分析。用R語言“UpSet”顯示7種可變剪接事件組成的預后風險模型和PASE的交互作用。將R語言的“UpSet”包用于基因的可視化。

1.4 構建SF和剪接事件共表達網絡

通過Pearson檢驗分析SF表達與PASE PSI值之間的相關性。Cytoscape 3.7軟件用于構建蛋白調節網絡。核心SF（Hub）由每個節點的數值和邊界定，Hub的表達水平由腫瘤組織和癌旁組織來確定，并分析SF與非SF的PASE構成的互作網絡。

1.5 統計學方法

采用Kaplan-Meier繪制生存曲線，用Logrank秩和檢驗評估高風險組、低風險組患者的總體生存率。采用受試者工作特征（receiver operating characteristic，ROC）曲線評估預后風險模型在1、3、5年生存期的預測能力，并繪制高、低風險熱圖。采用單因素和多因素Cox回歸分析評估臨床各變量及風險評分與患者預后的相關性。SF在癌旁組織和腫瘤組織中的表達差異用Wilcoxon非配對檢驗。

2 結果

2.1 COAD中可變剪接事件和PASE統計分析

398例COAD患者中共有9 085個基因發生了35 391次可變剪接事件，每個基因平均有4次可變剪接事件發生。其中5 635個基因發生13 087次ES，3 381個基因發生7 740次AT，2 692個基因發生6 653次AP，2 124個基因發生2 917次AA，1 833個基因發生2 524次AD，1 600個基因發生2 332次RI，137個基因發生138次ME。所有可變剪接事件中，ES（36.97%）為主要類型，ME（0.38%）是發生次數最少的類型。單個基因可能發生2個或更多可變剪接事件，為確定與COAD預后相關的PASE，利用單因素Cox回歸分析對可變剪接事件進行變量篩選（P＜0.05），結果顯示有1 811個基因發生2 515個PASE。

2.2 構建蛋白互作網絡并進行GO富集、KEGG通路分析

將821個PASE（P＜0.01）所對應的649個基因導入Cytoscape Reactome FI插件中，挖掘核心調控基因，發現鏈接點數最多的核心調控基因分別是TP53、RELA、UBE21。隨后將2 515個PASE中對應的1 810個基因進行GO富集和KEGG通路的生物信息學分析，發現GO富集共計47個BP（P＜0.01），主要包括RNA剪接、細胞周期調節、mRNA加工、組蛋白修飾；14個CC（P＜0.001），主要包括細胞器、線粒體基質、染色質和中心粒；17個MF（P＜0.001），包括泛素類蛋白轉移酶活性、轉錄調控、絲氨酸/蘇氨酸激酶活性、細胞黏附。KEGG通路有10個與生存相關的具有顯著性差異（P＜0.01）的通路，其中5個主要涉及泛素介導的蛋白水解通路、剪接體通路、HIF-1信號通路、脂肪酸代謝、結直腸癌相關信號通路。見圖1。

2.3 利用PASE構建COAD預后風險模型

數據中基因數量多而樣本量較少，因此用Lasso回歸分析進行變量篩選。為了更為精確和方便地描述可變剪接事件，本研究對COAD中每個可變剪接事件獨立編號，例如，PCNP-65964-AA，PCNP是基因名，65964是在COAD中剪接事件的唯一編號，AA是剪接事件類型。在7種可變剪接事件混合模型中，當最佳λ參數值為0.056時，共篩選出8個PASE。基于這8個PASE，用多因素Cox回歸分析進行變量篩選，最后構建出由8個PASE組成的預后風險模型（Risk score）。Risk score=（-17.504×HMGXB3-74054-RI）+（-104.165×PPP3CA-70095-ES）+（-1.947×KIAA1522-1632-AP）+（-6.171×SPINK1-73963-AT）+（-2.783×ZNF765-51718-AT）+（-13.186×PDCD2-78502-AA）+（-10.491×PTPRD-85850-ES）+（2.272×RAB3IP-23343-AP）。以風險值中位數（0.919）為最佳臨界值，將COAD患者劃分為高風險組（風險值＞0.919，173例）和低風險組（風險值≤0.919，174例）。熱圖結果顯示，隨著風險值的升高，PASE在2個組中的表達水平也明顯升高或降低。Kaplan-Meier曲線分析結果顯示，高風險組患者總體生存率較低風險組低，且2個組差異有統計學意義[風險比（hazard ratio，HR）=2.921，95%可信區間（confidence interal，CI）=1.890±4.512，P＜0.001）。用ROC曲線對預后風險模型進行預測性能的評價（COAD患者1、3、5年的生存率），結果顯示COAD患者1年生存率的曲線下面積[（area under curve，AUC）=0.860]好于3年（AUC=0.803）和5年（AUC=0.704）的。見圖2。

2.4 臨床相關參數和預后風險模型的單因素和多因素Cox回歸分析

為評估預后風險模型是否獨立于其他臨床變量，我們對其分別進行了單因素和多因素的Cox回歸分析。納入的臨床變量包括年齡、性別、息肉史、腫瘤浸潤、淋巴結轉移、遠處轉移、臨床分期。在單因素Cox回歸分析中，腫瘤浸潤、淋巴結轉移、遠處轉移、臨床分期、預后風險模型都與患者總體生存時間呈顯著負相關（P＜0.001）。經過多因素Cox回歸分析調整后，預后風險模型依然與患者總體生存時間呈顯著負相關（P＜0.001）。見表1和表2。

圖1 PASEs的生物信息學分析

圖2 PASE構建的COAD預后風險模型價值評價

2.5 PASE與其對應mRNA表達量的相關性分析

對預后風險模型的8個PASE對應的基因mRNA表達量在腫瘤組織和癌旁組織間進行比較，發現ZNF765、KIAA1522、HMGXB3、RAB3IP、PDCD2有差異（P＜0.01）。8個PASE對應的基因的生存曲線分析結果顯示，高風險組、低風險組之間沒有顯著性差異（P＞0.05）。見圖3和圖4。

2.6 生存相關SF的風險因素

從TCGA數據庫中提取到COAD相關的387個SF表達量，用單因素Cox回歸分析得到22個與生存相關的SF，其中只有RBM3、GRSF1、RBM47的HR＜1，屬于保護性因素；其他19個SF的HR＞1，屬于危險性因素。見表3。

2.7 PASE和SF共表達網絡

對2 515個PASE與22個生存相關的SF表達量進行相關性檢驗，并將結果導入到Cytoscape 3.71軟件中構建互作網絡，發現有420個PASE的PSI值與15個SF有相關性（P＜0.01），其中296個PASE與15個SF呈負相關；381個PASE與14個SF呈正相關。見圖5。

2.8 與生存相關的SF的mRNA表達量比較

我們將15個生存相關的SF在COAD癌旁組織和腫瘤組織中的表達量進行分析，發現有14個SF有顯著性差異（P＜0.05），其中C9orf78、CCDC130、CLK2、DHX38、GRSF1、ISY1、KHSRP、RBM17、RBM3、SART1和ZC3H18屬于危險因素；NRIP2、NOVA2和RBM47屬于保護性因素；NUMA1無差異（P＞0.05）。見圖6。

表1 預后風險模型單因素Cox回歸分析結果

表2 預后風險模型多因素Cox回歸分析結果

圖3 預后風險模型中8個PASE對應基因的mRNA表達量

圖4 預后風險模型中8個PASE對應基因的生存曲線分析

表3 22個生存相關SF的單因素Cox回歸分析結果

圖5 SF與生存相關剪接事件在COAD中的相關性分析

圖6 COAD癌旁組織和腫瘤組織中生存相關SF的mRNA表達量

3 討論

目前，COAD的可變剪接事件研究主要是以小樣本的測序方式進行，COAD的可變剪接事件分析還鮮有報道。本研究從TCGA SpliceSeq數據庫得到的COAD可變剪接事件和PSI值顯示，共有9 085個基因發生了35 391次可變剪接事件，ES（36.9%）為主要類型，ME（0.38%）是發生次數最少的類型。KAHLES等[7]的研究結果顯示，可變剪接事件在不同的腫瘤類型中出現的頻率各不相同，但ES都是最高的，與本研究結果相符。本研究將可變剪接事件與生存資料合并進行單因素Cox回歸分析，發現共有2 515個PASE，其所對應的基因有1 811個。“UpSet”結果顯示，部分基因至少有2種以上PASE發生，說明可變剪接事件在腫瘤發生、發展過程中起重要作用。將821個PASE（P＜0.01）所對應的649個基因導入到Cytoscape Reactome FI插件中，篩選出核心基因TP53、RELA、UBE21。TP53、RELA、UBE21是共表達網絡的核心節點，有望成為COAD治療的分子靶點。WU等[8]利用TCGA SpliceSeq數據庫下載和分析了肝癌的PASE，構建了8個風險評估模型，由ES事件組成的風險模型預測效果最好，AUC為0.898。為評估PASE在COAD中的預后診斷價值，本研究構建了預后風險模型，混合剪接事件構成的預后風險診斷模型高、低風險有顯著性差異（P＜0.001），ROC曲線的AUC為0.860（1年生存率），在COAD患者的風險分層中表現準確。為了評估預后風險模型是否獨立于其他臨床變量，本研究分別采用單因素和多因素Cox回歸分析評估臨床病理參數和預后風險模型對COAD患者預后的影響，發現預后風險模型與患者總體生存時間呈顯著負相關（P＜0.001）。對構成預后風險模型的8個PASE對應的基因的癌旁組織與腫瘤組織mRNA表達量進行分析，發現并無統計學差異（P＞0.05）。基因和其可變剪接事件可能執行不同的生物學功能，其具體的機制尚未被完全闡明，需要進行進一步探討。

在腫瘤中，SF表達異常可能會形成特定促癌剪切異構體，從而導致癌癥發生。為了闡明SF和PASE之間復雜的調控關系，本研究構建了互作網絡，有420個PASE的PSI值與15個SF有相關性（P＜0.01），其中296個PASE與15個SF存在負相關關系；381個PASE與14個SF存在正相關關系。這說明有的SF具有雙重調控功能，有的PASE同時受到不同的SF調控。在互作網絡中，SF細胞分裂周期樣激酶2與SF富含絲氨酸/精氨酸剪接因子3 2個同類型的PASE存在相反的調控關系，說明同一基因的同類型可變剪接事件可能執行不同的生物學功能。SRSF3和CLK2的mRNA表達量幾乎無相關性，而CLK2的4個可變剪接事件CLK2-8053-AP、CLK2-8054-AP、CLK2-8055-AA和CLK2-8056-AA則不屬于PASE，因此基因的mRNA表達量并不能很好地描述基因的功能，基因的生物學功能可能要依靠占優勢的可變剪接事件進行描述。SF在調控可變剪接事件中起到了重要的作用，其核酸序列的突變或者表達水平的改變都可能影響可變剪接事件[9]，但SF如何調節可變剪接事件的類型目前尚不清楚。通過上調或下調SF的表達水平都可能促進或者抑制腫瘤的進展[10]。分析COAD患者生存資料得到的22個與生存相關的SF中，大部分SF高表達時患者生存期縮短，這些SF的表達水平可能是被類似的可變剪接事件所調控。相關性分析結果顯示，15個與生存相關的SF中有14個在癌旁組織和腫瘤組織中有差異（P＜0.05），多數在腫瘤組織中表達量上調，說明在COAD發生、發展過程中與PASE相關的SF過表達，對患者預后不利。各組間核心調控因子CLK2在COAD的臨床分期、T分期、組織分級中表達有差異（P＜0.05），而在N分期和M分期中無差異（P＞0.05），說明SF的高表達多與患者預后不良有關。

可變剪接事件的發生更像是細胞為了應對外來的刺激信號而快速產生的特殊表型，用于執行特定的功能。相對于龐大的基因組，目前檢測到的剪接產物很少，剪接機制的研究還處于初級階段。本研究嘗試用多個可變剪接事件類型構建預后風險模型，得到了滿意的效果。對于數據而言，可變剪接事件可單獨成為腫瘤診斷和預后的標志物。然而，以單一類型可變剪接事件作為標志物是不夠充分的，聯合多個可變剪接事件來提高診斷和預后的敏感性，也符合臨床研究的發展趨勢。本研究還構建了SF與PASE之間的互作網絡，這為COAD的預后風險模型的構建提供了幫助。