顧云婧,潘以紅,朱東月,朱 平

(江南大學理學院,中國江蘇無錫214122)

卵巢癌是女性常見的癌癥之一,居于女性癌癥死因的第四位,2012年全球共診斷出新增卵巢癌患者約238 700例[1]。研究顯示,由于早期癥狀不明顯,約70%的卵巢癌患者被診斷時已是晚期,晚期患者的5年生存率低于30%[2]。目前,卵巢癌的治療手段主要是手術切除和化療,但預后情況并不理想[3],大多數(shù)晚期患者會在18個月內復發(fā)[4]。因此,迫切需要更深刻地了解卵巢癌的預后分子機制,尋找新型治療靶標,以開發(fā)更有效的治療方法來改善卵巢癌的預后效果。

已有研究顯示,基因表達失調對腫瘤的發(fā)生和發(fā)展有重要影響,通過分析基因表達模式可以發(fā)現(xiàn)具有預后價值的癌癥生物標志物[5]。目前,這些研究主要集中在對編碼基因表達水平的分析,而對mRNA可變剪接(alternative splicing,AS)的系統(tǒng)研究還較為缺乏?？勺兗艚邮菃蝹€前體RNA分子產生多種成熟mRNA的機制,mRNA的可變剪接是高等真核生物擴大蛋白質組多樣性的重要基因調節(jié)機制之一[6～7]。因而,可變剪接模式的變化可能導致蛋白質功能的變化并影響基于這些功能的重要活動。近年來,大量的基因組學和功能性研究表明可變剪接與癌癥的發(fā)生發(fā)展密不可分,可變剪接事件參與了一系列的致瘤過程,包括侵襲、增殖、血管生成、免疫和轉移[8]。據(jù)報道,含有外顯子v8-10的CD44變體與卵巢癌的轉移和預后有關[9]；剪接變體AIMP2-DX2可能是治療卵巢癌化療耐藥性的有效輔助靶點[10]。癌癥研究的新趨勢表明,可變剪接在癌癥治療方面具有臨床潛力[11～14],可以作為癌癥診斷和治療的生物標志物,并作為新型治療藥物的潛在靶標。

本文根據(jù)卵巢癌患者的全基因組可變剪接數(shù)據(jù)及臨床數(shù)據(jù),篩選出與總生存率(overall survival,OS)相關的AS事件,并融合其對應基因的蛋白質相互作用(protein-protein interaction,PPI)數(shù)據(jù)構建卵巢癌預后加權網絡；對預后加權網絡進行拓撲分析,篩選有意義的AS事件作為卵巢癌預后特征；利用預后特征建立預后預測模型。結果顯示模型可用于預測卵巢癌患者的臨床結果,顯示了可變剪接在卵巢癌中的預后價值。

1 材料與方法

1.1 數(shù)據(jù)來源

卵巢癌患者的可變剪接數(shù)據(jù)下載自TCGASpliceSeq數(shù)據(jù)庫(http://projects.insilico.us.com/TCGASpliceSeq/)[15]。組數(shù)據(jù)庫(The Cancer Genome Atlas,TCGA)中mRNA可變剪接模式的公開數(shù)據(jù)庫資源。本文從中選取了393個具有臨床數(shù)據(jù)的卵巢癌患者樣本,生存時間均大于30 d,且樣本的可變剪接表達估計值(percent spliced in,PSI)百分比大于75%。同時,從TCGA數(shù)據(jù)庫中額外下載了16個卵巢癌患者樣本的RNA-seq數(shù)據(jù),并利用SUPPA軟件[16]計算樣本中可變剪接事件的PSI值。

1.2 預后相關AS事件的篩選

對于臨床數(shù)據(jù),由于每個病例隨訪時間不同,發(fā)生終點事件的可能性也不相等,因此本文選用Cox單因素回歸分析篩選預后相關的AS事件(prognostic-associated alternative splicing events,PASEs)。

Cox單因素回歸模型[17]形式:

其中,h(t|X)是具有協(xié)變量X的樣本在t時的風險函數(shù),t為生存時間；h0(t)是協(xié)變量為0時的風險函數(shù)；β為模型中待估計的回歸系數(shù)。

回歸系數(shù)β＞0時,協(xié)變量的取值越大,風險函數(shù)h(t|X)的值越大,表示病人死亡的風險越大；回歸系數(shù)β=0時,表示協(xié)變量對風險函數(shù)h(t|X)沒有影響；回歸系數(shù)β＜0時,協(xié)變量的取值越大,風險函數(shù)h(t|X)的值越小,表示病人死亡的風險越小。

對于每個可變剪接事件,將患者根據(jù)該事件PSI值的中位數(shù)進行分類,分為高、低表達兩組進行單因素回歸分析,計算危險比(hazard radio,HR)及95%置信區(qū)間(confidence interval,CI),以P＜0.05為標準篩選出與卵巢癌預后相關的AS事件。

1.3 卵巢癌預后加權網絡的構建

考慮PASEs之間的調控關系,計算預后相關AS事件之間的相關性,并融合其對應基因的PPI數(shù)據(jù)建立加權網絡。對于每一個PASE,計算它與其余所有PASEs之間的皮爾遜相關性系數(shù)(Pearson correlation coefficient,PCC)。選取PCC值排名前5%,且相應P值小于0.01的AS事件對構建網絡。相關性矩陣MAS由PASEs中符合條件的PCC構成,矩陣元素M(i,j)代表第i和j個事件之間PCC的絕對值。按如下方式對相關性矩陣進行標準化,標準化后的矩陣為

其中EAS(i,i)為矩陣MAS第i行元素的和。

除了相關性關系外,本文在網絡構建中融入蛋白質相互作用數(shù)據(jù)。將預后相關的AS事件對應基因(genes of prognostic-associated alternative splicing events,PASEGs)提交至STRING在線數(shù)據(jù)庫(https://string-db.org)[18],以得分大于0.4為標準,選取PASEGs之間的PPI數(shù)據(jù)。隨后,根據(jù)PASEGs之間PPI的置信度構造矩陣MPPI,對其進行標準化:

其中EPPI(i,i)為矩陣MPPI第i行元素的和。

將AS事件對的相關性關系和其對應基因的PPI相結合,構成網絡權重:

1.4 預后加權網絡的拓撲分析及PASEGs的功能富集分析

對預后加權網絡進行網絡拓撲特征分析,以挖掘拓撲重要性強的網絡節(jié)點。網絡特征分析包含網絡節(jié)點的度(degree)、特征路徑長度(characteristic path length,CPL)、聚類系數(shù)(clustering coefficient,CC)和小世界指數(shù)(small-world index,SW)。節(jié)點的度為該節(jié)點與其他節(jié)點直接連接的邊數(shù)。網絡的CPL為網絡中所有節(jié)點對之間的最短路徑長度的平均值。節(jié)點的CC為與該節(jié)點直接連接的節(jié)點中實際存在的邊與所有可能存在的邊的比例,網絡的CC為網絡中所有單個節(jié)點的CC平均值。通過與相同規(guī)模的ER(Erdos-Rényi)隨機網絡比較可定義小世界網絡[19],即該網絡需同時滿足以下條件:

其中,CCrandom和CPLrandom分別表示1 000個與目標網絡相同規(guī)模的ER隨機網絡的CC、CPL平均值。

進一步地,網絡小世界指數(shù)[20]被定義為:

當網絡的小世界指數(shù)бSW大于1時,該網絡為小世界網絡。這意味著,該網絡的節(jié)點高度聚集,網絡的平均最短路徑長度幾乎與隨機網絡一樣低,該生物網絡具有高效的信息傳遞能力。

通過DAVID在線工具[21]對PASEGs進行GO(gene ontology)功能富集分析。GO分析包括分子功能(molecular function,MF)、生物過程(biological process,BP)和細胞組分(cellular component,CC)3個部分。同時,運用KOBAS數(shù)據(jù)庫(http://kobas.cbi.pku.edu.cn/)[22]對PASEGs進行KEGG(kyoto encyclopedia of genes and genomes)通路富集分析,P＜0.05的通路被認為是顯著富集的通路。

1.5 預后特征的篩選及預測模型的建立

卵巢癌預后加權網絡的拓撲性質顯示了AS事件在卵巢癌中的預后價值。為了預測卵巢癌患者的臨床結果,本文根據(jù)網絡節(jié)點的度篩選預后特征,并構建預后預測模型。首先選擇預后加權網絡中節(jié)點度排名前20的AS事件作為預后特征,隨后通過Cox比例風險模型(proportional hazards model)構建卵巢癌預后預測模型,預后指數(shù)(prognosis index,PI)的計算公式為:

其中,x和β表示預后特征及其回歸系數(shù),n表示預后特征的個數(shù)。

根據(jù)預后指數(shù)PI的中值將患者分為高、低風險兩組,運用Kaplan-Meier圖表估計生存函數(shù),分析兩組之間的預后差異。此外,ROC(receiver operating characteristic)曲線[23]和 AUC(area under curve)值被用來評估預后預測模型的效率,AUC值越大,預測模型的正確率越高。

2 結果及分析

2.1 預后相關的可變剪接事件

可變剪接事件一般分為7種[24],分別為外顯子跳躍(exon skip,ES)、互斥可變外顯子(mutually exclusive exon,ME)、內含子保留(retained intron,RI)、可變 3′剪接位點(alternative acceptor site,AA)、可變 5′剪接位點(alternative donor site,AD)、可變初始外顯子(alternative promoter,AP)和可變終末外顯子(alternative terminator,AT)。

利用Cox回歸對393個樣本的所有剪接事件分別進行單因素生存分析,共篩選出290個與卵巢癌OS顯著相關的AS事件。其中,AA事件16個,AD事件26個,AP事件69個,AT事件33個,ES事件116個,ME事件1個,RI事件29個。P值排名前10的AS事件見表1。

2.2 預后加權網絡的拓撲分析

從STRING數(shù)據(jù)庫中獲取了118個PASEGs的158對互作關系,與相應AS事件對的PCC結合后構建了卵巢癌預后加權網絡。該網絡共包含281個節(jié)點,2 284條邊。圖1是利用Cytoscape軟件[25]繪制的權重大于0.1的部分網絡圖。網絡拓撲分析表明,預后加權網絡的聚類系數(shù)為0.255,特征路徑長度為2.498。1 000個相同節(jié)點數(shù)和邊數(shù)的ER隨機網絡的平均聚類系數(shù)為0.058,平均特征路徑長度為2.309。與ER隨機網絡相比,預后加權網絡滿足公式(5)和公式(6),因此預后加權網絡為小世界網絡,由公式(7)可得其小世界指數(shù)為4.064。

表1 排名前10的預后相關AS事件Table 1 Top 10 prognostic-associated alternative splicing events

圖1 卵巢癌預后加權網絡的代表性結果節(jié)點的大小代表該AS事件度的大小,邊的粗細代表兩AS事件間的權重大小。Fig.1 Representative results of prognostic weighted network for ovarian cancerThe size of the node represents the degree of the AS event,and the thickness of the edge represents the weight between the two AS events.

2.3 PASEGs的功能富集分析

通過DAVID在線工具對PASEGs進行GO富集分析,結果表明,PASEGs主要富集于細胞核、核質、RNA聚合酶II啟動子轉錄調控、DNA轉錄調節(jié)、細胞凋亡、蛋白質結合、核酸結合等29個生物學過程(圖2A)。顯然,這些生物學過程與細胞中遺傳物質的轉錄、蛋白質的合成等具有密切的相關性,可對細胞的活動起到一定的調控作用,從而影響卵巢癌的發(fā)生發(fā)展。

同時,運用KOBAS數(shù)據(jù)庫對PASEGs進行KEGG通路富集分析,結果顯示,PASEGs主要富集在腸道免疫網絡、p53信號通路、凋亡途徑、核糖體相關途徑等24條通路上。這與GO富集分析結果一致,表明細胞中遺傳物質的轉錄、蛋白質的合成及細胞凋亡等功能的異常在卵巢癌發(fā)生發(fā)展的過程中具有重要作用,可能是影響卵巢癌預后的重要途徑。排名前10的顯著富集通路見圖2B。

2.4 可變剪接預后特征及預后模型

根據(jù)卵巢癌預后加權網絡中節(jié)點的度對網絡節(jié)點進行排序,選擇排名前20的AS事件作為可變剪接的預后特征(表2),并運用Cox回歸模型建立預后預測模型。預后模型的預后指數(shù)由20個預后特征的表達值和Cox回歸分析得出的回歸系數(shù)線性組合構建。根據(jù)公式(8)計算每個患者的預后指數(shù),并將患者根據(jù)預后指數(shù)的中位數(shù)分為高風險組和低風險組。圖3A的Kaplan-Meier生存曲線顯示,高風險組的患者生存率明顯低于低風險組(logrankP＜0.000 1),表明預后特征能夠顯著區(qū)分不同預后的患者。同時,為了評估模型的靈敏度和特異性,我們進行了時間依賴性ROC曲線分析,結果如圖3B所示,預后模型的AUC值為0.767。用此模型對TCGA中的16個樣本進行測試,測試集的AUC值為0.785,顯示了模型良好的預測性能。

3 討論

隨著近幾年對可變剪接的深入研究,已有文獻證明異?？勺兗艚訉τ诼殉舶┑陌l(fā)生發(fā)展具有重要影響[9～10]。本文確定了一批與卵巢癌預后相關的AS事件,并通過融合PPI數(shù)據(jù)構造了卵巢癌預后加權網絡。此外,通過可變剪接預后特征建立了預后模型,模型顯示,這些預后特征可以較好地預測卵巢癌患者的預后情況。因此,這些預后特征可能是卵巢癌患者潛在的預后因素,可作為癌癥預后的生物標志物和治療靶點。

本研究的拓撲分析表明卵巢癌預后加權網絡為小世界網絡,小世界指數(shù)為4.064。這意味著,一方面,網絡中的節(jié)點高度聚集:當一個節(jié)點連接到另外兩個節(jié)點時,后兩個節(jié)點也傾向于彼此直接連接；另一方面,網絡中的平均最短路徑長度幾乎與隨機網絡一樣低,說明預后加權網絡具有高效的信息傳遞能力。因此,本文利用網絡節(jié)點的拓撲性質篩選了可變剪接預后特征,并以此構建了卵巢癌的預后預測模型,模型顯示了良好的預測性能。

表2 用于構建預后模型的AS預后特征Table 2 Prognostic signatures of AS events involved in prognostic model

圖2 PASEGs的富集分析(A)PASEGs的GO富集分析；(B)PASEGs的KEGG富集分析中排名前10的通路。Fig.2 The enrichment analysis of PASEGs(A)GO enrichment analysis of PASEGs；(B)KEGG enrichment analysis of PASEGs(top 10).

圖3 可變剪接預后特征的Kaplan-Meier生存曲線及ROC曲線(A)預后特征的Kaplan-Meier生存曲線；(B)預后模型的ROC曲線。紅色曲線表示訓練集,紫色曲線表示測試集。Fig.3 Kaplan-Meier plot and ROC curve of alternative splicing prognostic signatures(A)Kaplan-Meier plot of prognostic signatures；(B)ROC curve of prognostic model.Red curve and purple curve indicate training set and test set,respectively.

在可變剪接預后特征所對應的基因中,AP2B1[26]、ATF5[27]、BCL2L11[28]、ING3[29]、SOD2[30]、USP36[31]等已被證實在卵巢癌細胞增殖和生長中具有重要功能,這表明拓撲學上重要的節(jié)點往往在疾病中發(fā)揮關鍵作用。例如:與化學敏感性細胞或組織相比,AP2B1在卵巢癌化學抗性細胞或組織中始終過表達[26]；ATF5在大多數(shù)上皮性卵巢癌樣本中高表達,且與晚期臨床分期和腫瘤分化差異顯著相關[27]；BCL2L11介導的AKT磷酸化水平可影響卵巢癌細胞中順鉑的耐藥性[28]。雖然現(xiàn)階段其余基因在卵巢癌中的功能尚未得到證實,但已被證明與癌癥相關。例如:DUSP28在胰腺癌中高度表達并發(fā)揮關鍵作用,靶向DUSP28可以抑制惡性胰腺癌的發(fā)展[32]；IL17RC缺失的結直腸腫瘤細胞表現(xiàn)為侵襲、生長和轉移,提示IL17RC可作為結直腸癌的預后標志物[33]。因此,這些基因可能是卵巢癌的潛在疾病基因,值得進一步研究。

在PASEGs的GO功能富集分析中,我們發(fā)現(xiàn)了幾類一致的GO注釋。例如:預后相關基因主要參與細胞核、核質等細胞組分,與RNA聚合酶II啟動子轉錄調控、DNA轉錄調節(jié)、蛋白質結合、核酸結合等生物學過程一致。顯然,這些生物學過程與細胞中遺傳物質的轉錄、蛋白質的合成等具有密切的相關性,可對細胞活動起到一定的調控作用,進而影響卵巢癌的發(fā)生發(fā)展。PASEGs還參與了細胞凋亡過程,這與KEGG通路富集分析的結果相互印證。KEGG通路富集分析表明,PASEGs主要顯著富集在腸道免疫網絡、p53信號通路、凋亡途徑、核糖體相關途徑等,提示這些通路可能是與卵巢癌AS預后相關的重要途徑。已有研究顯示,通過激活MDM2(或HDM2)和MDM4的可變剪接,可誘導p53通路發(fā)生改變[34]。此外,卵巢癌細胞的凋亡程序受到p53家族轉錄因子的調控[35],而p53的促凋亡功能是癌癥治療的有效靶標[36]。因此,在卵巢癌的發(fā)生發(fā)展機制方面深入研究這些相互作用的AS預后相關基因將是非常有意義的。

總的來講,本文在全基因組水平上對卵巢癌的可變剪接事件進行了綜合挖掘,揭示了AS事件在卵巢癌中的預后價值。特別是在預后預測中,對mRNA可變剪接的分析可能揭示了新的癌癥驅動因素,并以此提出了理想的卵巢癌預后預測特征。