基于TCGA數據庫分析胃癌可變剪接與腫瘤免疫的關系

顧琦晟，張米粒，曹 燦，李繼坤

上海交通大學附屬第一人民醫(yī)院普外科，上海201620

胃癌是常見的惡性腫瘤之一。根據2018 年全球癌癥統計數據，胃癌新發(fā)病例超過100 萬例，估計死亡病例783 000 例；在所有癌癥類型中，胃癌的發(fā)病率（5.7%）和病死率（8.2%）分別居第五位和第三位［1］。雖然胃癌的早期診斷和整體治療水平不斷提高，但是進展期胃癌的療效仍不盡如人意。近年來，除了傳統的手術治療、放射治療、化學治療，腫瘤免疫治療也展現出了無限潛力。以程序性死亡蛋白1（programmed death-1，PD-1）單抗為首的腫瘤免疫檢查點抑制劑對黑色素瘤、血液系統腫瘤等非實體腫瘤顯示出了極高的臨床治療活性［2-3］；但是，對進展期胃癌的治療效果較差，這與胃癌免疫逃避的潛在機制有著密切關系［4］。因此，研究胃癌的免疫逃避機制，找到增強胃癌免疫原性的方法，是現今進展期胃癌免疫治療方面亟待解決的一個問題。

可變剪接又稱為選擇性剪接，是一項精密、復雜的轉錄后調控機制，是基因的mRNA 前體轉錄所產生的RNA 外顯子以多種方式通過RNA 剪切進行重連的過程；從而，單個mRNA 可產生不同的蛋白質構體，是增加mRNA 復雜性和蛋白質多樣性的原因之一［5］。研究［6-7］發(fā)現，在腫瘤的發(fā)生、發(fā)展過程中，伴隨著基因的可變剪接異常產生，其中不乏許多已經證實的癌基因和抑癌基因的變化，并進一步推進腫瘤的發(fā)生、發(fā)展。

目前，有關胃癌可變剪接和腫瘤免疫浸潤的關系尚待深入研究。本研究旨在利用生物信息學方法，探討胃癌中的可變剪接情況，及其與胃癌表型、臨床病理學因素以及腫瘤免疫浸潤的關系，為胃癌的免疫治療提供新的思路。

1 資料與方法

1.1 數據采集

在TCGA（The Cancer Genome Atlas）數據庫中下載375 例胃癌患者及32 例配對癌旁正常組織的轉錄組、基因組數據和所有納入研究的胃癌患者的臨床信息。在SpliceSeq 數據庫中下載與TCGA 的胃癌樣本相對應的樣本剪接百分比（percent-spliced in，PSI）矩陣，納入具有PSI 值的樣本占比<75%的剪接事件，共包含48 141 個剪接事件以及415 例胃癌組織和37 例配對癌旁正常組織樣本， 以 備 后 續(xù) 分 析［8］。 在GEO （Gene Expression Omnibus）數據庫中下載包含192 例胃癌患者轉錄組芯片數據的數據集GSE15459［9］。

1.2 數據分析及統計學處理

1.2.1 胃癌樣本可變剪接數據的主成分分析 利用missMDA 包，基 于EM-PCA （Expectation Maximization Algorithm for Principal Component Analysis）算法對參數求極大似然估計（Maximum Likelihood Estimate，MLE），對缺失數據進行插補［10］。利用FactoMineR 包對數據進行PCA 降維分析，利用ggbiplot 和factoextra 包對PCA 數據進行可視化［11］。

1.2.2 利用一致性聚類方法對胃癌樣本的可變剪接進行無監(jiān)督聚類 利用ConsensusClusterPlus包的一致性聚類方法對PSI矩陣進行k-中心點聚類（k-medoid），聚類值k取2～10，根據聚類效果選擇具有較高聚類穩(wěn)定性的k值［12］。

1.2.3 基因集通路富集分析 使用limma 包的競爭性基因集測試算法（Competitive Gene Set Test Accounting for Inter-gene Correlation，CAMERA），以分子特征數據庫（Molecular Signatures Database，MsigDB）的Hallmark 基因集為功能基因集，分析相關樣本集與其他樣本的差異，分析與該樣本集密切相關的通路［13-14］。錯誤發(fā)現率（false discovery rate，FDR）<0.05 表示差異有統計學意義。利用GSVA 包的基因集變異分析（Gene Set Variation Analysis，GSVA）方法，將基因在不同樣品間的表達水平矩陣轉化為基因集在樣品間的表達水平矩陣，從而分析該基因集表達水平在各個樣本間的差異［15］。

1.2.4 可變剪接相關基因的篩選 對404 個經實驗確認的剪接相關因子與胃癌可變剪接事件變化的關系進行研究［16］。根據基因的表達水平，以中位數為界，分為高表達組和低表達組，利用limma 包進行剪接事件的差異分析。以logFC>0.3 且BH 法修正后P 值<0.05 的差異可變剪接事件作為潛在受剪接因子表達水平影響的事件，篩選影響超過20%總剪接事件且相對mRNA 平均表達水平>4的基因，作為胃癌中潛在調控可變剪接的剪接因子。

1.2.5 與胃癌免疫微環(huán)境浸潤的關系分析 利用CIBERSORT 方法和LM22 免疫細胞特征矩陣，對胃癌免疫細胞浸潤豐度情況進行反卷積分析［17］。從美國國家腫瘤研究所基因數據庫（National Cancer Institute Genomic Data Commons）獲取TCGA 胃癌樣本的腫瘤免疫相關特征數據［18］。

1.3 統計學分析

采用R 3.6.2 軟件進行統計學分析。定性資料的比較采用Kruskal-Wallis H 檢驗；定量資料以±s 表示，采用ANOVA檢驗進行比較。P<0.05表示差異有統計學意義。

2 結果

2.1 胃癌可變剪接的整體情況

提取SpliceSeq 的初始矩陣，內含48 141 個剪接事件和452 例樣本。制定嚴格的過濾標準，將有效值占比>90%且PSI 標準差>0.1 的事件納入后續(xù)研究，同時排除有效值占比<80%的樣本7 例（圖1A、1B）。經過濾，獲得8 649 個剪接事件和445 例樣本，其中腫瘤樣本409 例，癌旁正常組織樣本36例。

通過觀察4 051 個基因和相關8 649 個剪接事件的大體分布情況，發(fā)現平均每個mRNA 具有超過2 種不同形式的剪接變化，其中外顯子跳躍（exon skipping，ES）、可變初始外顯子（alternate promoter，AP）和可變終末外顯子（alternate terminator，AT）是胃癌中最常見的剪接類型，5′端和3′端可變剪接、內含子保留事件相對較少，互斥外顯子事件最少（圖1C、1D）。

圖1 TCGA中胃癌可變剪切概況Fig 1 Overview of alternative splicing of gastric cancer in TCGA

2.2 胃癌的可變剪接和主成分分析及聚類分析

使用EM-PCA 算法對剪接事件進行缺失值插補后，通過主成分分析降維。利用碎石圖展示剪接事件的主成分對剪接全景的貢獻度，其中PC1（主成分1）作為最重要的主成分，解釋了剪接全景10.7%的方差值（圖2A）。通過PCA 圖發(fā)現PC1 可以明顯區(qū)別腫瘤組織與癌旁正常組織（Wilcox test，P<0.01），提示腫瘤與癌旁正常組織中存在許多有顯著差異的剪接事件（圖2B、2C）。

利用ConsensusClusterPlus 包對剪接事件矩陣進行無監(jiān)督聚類（圖3A、3B）?；谝恢滦跃仃嚭虳elta 區(qū)域結果，認為k=4時穩(wěn)定性較佳，且能夠減少過擬合效應。聚類簇C3、C1、C2+C4 在PC1 水平明顯區(qū)分，而C2 和C4在PC2 水平區(qū)分（圖3C）。C3 與癌旁正常組織的分布距離較大，提示該聚類簇包含的樣本具有較為顯著的腫瘤相關異常剪接發(fā)生。

如表1所示，通過比較4個亞組之間患者年齡、TMN分期、病理分期、組織學分期、Lauren 分型、性別、種族等臨床特征，發(fā)現可變剪接的亞型與淋巴結轉移分期、性別無關，而與其他臨床特征有較大的關系。其中C3 型患者平均患病年齡偏大、遠處轉移可能更高、美國癌癥聯合委員會（American Joint Committee on Cancer，AJCC）分期相對較差；而C4 型患者患病年齡較小，腫瘤浸潤相對更深，多表現為彌漫性胃癌，組織學多表現為G3低分化。

圖2 通過PCA對可變剪接事件進行降維分析Fig 2 Dimensional reduction analysis of alternative splicing events by using PCA

圖3 根據可變剪接事件對胃癌樣本進行一致性聚類分析Fig 3 Consensus clustering of gastric cancer samples based on alternative splicing events

表1 基于可變剪接分類的4種亞型胃癌患者的臨床資料比較Tab1 Comparison of clinical information of patients among the 4 subtypes of gastric cancer based on alternative splicing classification

2.3 胃癌的可變剪接亞型與腫瘤標志特征和免疫微環(huán)境的關系

對4種可變剪接亞型基于包含50個主要生物活動類型的MSigDB 的標志基因集（Hallmark Gene Set）分別進行CAMERA分析，結果顯示C1型的上皮-間質轉化（epithelialmesenchymal transition，EMT）通路顯著下調；C2型的EMT通路相對上調；C3型蛋白分泌通路相對上調；C4型的細胞周期相關E2F目標基因顯著下調。該結果提示4種亞型的腫瘤生物學特性各不相同，C1型腫瘤的侵襲轉移能力相對較弱，而C4型腫瘤的增殖速度較其他類型更慢（圖4）。

如圖5所示，結合TCGA泛癌免疫全景研究結果，探索了可變剪接亞型和一些腫瘤特征的關系。結果顯示，C2 型與C4 型的基質比率、淋巴細胞比率和轉化生長因子-β（transforming growth factor-β，TGF-β）應答的特征較強；C1型的巨噬細胞調節(jié)水平較低；C3型的腫瘤浸潤性淋巴細胞（tumor infiltrating leukocyte，TIL）區(qū)域比例較低；C4 型的增殖、創(chuàng)傷修復、突變水平較低，而淋巴細胞浸潤等相對較高。

圖4 標志基因集在4個可變剪接亞型中的最顯著相關通路Fig.The most significantly related pathways of hallmark gene set in the 4 subtypes based on splicing events

基于22 種免疫細胞的簽名矩陣，利用CIBERSORT對TCGA 的轉錄組數據進行反卷積，得到胃癌組織各免疫細胞組分的比例，并通過歸一化去除腫瘤純度帶來的影響，比較各個亞型間免疫細胞組分的差異。如圖6 所示，可變剪接亞型間的腫瘤微環(huán)境組分差異明顯。C1 型的M2 型巨噬細胞豐度較低；C3 型的自然殺傷（natural killer，NK）細胞幾乎完全處于靜息狀態(tài)，調節(jié)性T 細胞豐度極低，T 細胞受體（T cell receptor，TCR）香農指數較低，而有相對高水平的活化記憶性T 細胞和M2 型巨噬細胞；C4 型活化NK 細胞、肥大細胞豐度較高，TCR 香農指數相對最高。

2.4 胃癌的可變剪接情況與剪接因子表達水平的關聯

如圖7A 所示，對胃癌的剪接因子mRNA 表達數據進行了PCA 降維，發(fā)現胃癌的4 個亞型聚集效果良好，提示剪接因子的表達水平變化是影響可變剪接事件發(fā)生的主要因素。如圖7B 所示，絕大多數剪接因子的高表達與PC1水平的降低相關，而CELF2、QKI、RBMS1等基因的表達模式和其他剪接因子相反，一定程度上起到拮抗異?？勺兗艚邮录l(fā)生的作用。

圖5 可變剪接亞型與腫瘤特征的關系Fig.Association of subtypes based on splicing events and tumor characteristics

圖6 可變剪接亞型與腫瘤微環(huán)境的關系Fig.Association of subtypes based on splicing events and tumor microenvironment

圖7 通過PCA分析可變剪接亞型與剪接因子的關系Fig 7 Association of subtypes based on splicing events and splicing factors using PCA

通過差異基因分析對剪接因子進行研究，篩選log2FC>1 且FDR<0.05 的差異表達剪接因子。結果顯示C1 型與RBM25、SREK1 的上調關系密切，且受到眾多剪接因子的共同影響；C2 型與RAVER1、RBM42、ALYREF等基因上調相關；C3 型與TIA1、DDX39B、PAXBP1 等基因高表達相關；C4 型與QKI、CELF2、RBMS1 等基因的高表達相關（圖8）。

2.5 胃癌的剪接因子表達水平與腫瘤抗原提呈水平的關系

選擇圖8中陳列的所有剪接因子基因作為樣本的剪接因子簽名。此外，PSMB5、PSMB6、PSMB7、PSMB8、PSMB9、 PSMB10、 TAP1、 TAP2、 ERAP1、 ERAP2、CANX、CALR、PDIA3、TAPBP、B2M、HLA-A、HLA-B、HLA-C 共18 個基因是明確的抗原提呈細胞相關基因，采用其作為基因集分析胃癌樣本的抗原提呈水平［19］。

圖8 可變剪接亞型與主要剪接因子中位值表達水平的歸一化后熱圖Fig 8 Median normalized expression level of splicing factors in the 4 subtypes based on splicing events

利用上述剪接因子基因集和抗原提呈相關基因集，采用GSVA 分別對樣本進行評分，隨后進行相關性分析，以Pearson 相關系數表示其相關性。結果顯示（圖9）：剪接因子基因集和抗原提呈相關基因具有明顯的相關性（Pearson R=0.44，P = 0.000）。利用GEO 數據庫中包含192 例胃癌轉錄組芯片數據的隊列GSE15459 進行驗證，結果與TCGA 胃癌數據一致，提示剪接因子家族的表達對腫瘤的抗原提呈過程有重要的潛在作用。

圖9 胃癌剪接因子基因表達特征與抗原提呈基因表達特征的關系Fig 9 Scatter plots of gene expression signatures between antigen presentation and splicing factors

3 討論

可變剪接過程在正常機體中受到嚴密的調控，剪接網絡在機體中有可靠的魯棒性。而在腫瘤發(fā)生、發(fā)展的過程中，隨著腫瘤的基因組不穩(wěn)定等內在環(huán)境壓力和免疫微環(huán)境等改變帶來的外在環(huán)境壓力的產生，腫瘤細胞利用可變剪接這一重要的轉錄后調節(jié)水平變化調整自身的活性和功能，從而起到增強自身增殖信號、促進侵襲轉移、逃避免疫監(jiān)視等作用［20］。然而，目前有關可變剪接與胃癌腫瘤免疫相關的研究較為有限，值得進一步深入研究。

本研究利用剪接事件矩陣對TCGA 數據庫的胃癌可變剪接事件全景進行了分析。結果顯示，腫瘤與癌旁正常組織有顯著的剪接差異。利用ConsensusClusterPlus 包對剪接事件矩陣進行無監(jiān)督聚類，將所有樣本分為4個亞型，可以認為C3 型在剪接事件水平上的異型性最強，C1型次之，C2型和C4型與癌旁正常組織最為相近。通過分析4 個亞型的臨床特征，發(fā)現C3 型和C4 型的特征最明顯：C3 型的病理分期較差，相對容易出現遠處轉移；C4型組織學表現為低分化，多為彌漫性胃癌。

基于Dunn 等的免疫編輯（immunoediting）學說，本研究中的腫瘤組織多處于免疫平衡（immune equilibrium）或免疫逃避（immune escape）階段［21］。胃癌的發(fā)生、發(fā)展與慢性炎癥密切相關，腫瘤多已能夠免疫耐受。本研究發(fā)現4 個亞型的免疫逃避機制較為不同，C1 型與C2 型的抗原提呈水平較高，但C1 型的髓系細胞和淋巴細胞浸潤程度都較低，而C2 型中M2 型巨噬細胞和調節(jié)性T 細胞等免疫抑制細胞的高豐度可能誘導了T 細胞無能，使得高豐度的CD8+T 細胞無法完成腫瘤殺傷的功能。C3 型和C4 型的抗原提呈水平較低，在C3 型中與較低的TCR多樣性相關，而C4 型與基因組穩(wěn)定、體細胞突變較少有關?？傮w上，4 個亞型潛在的免疫逃避機制分別為：C1型與低免疫浸潤水平相關；C2型與高免疫抑制環(huán)境相關；C3 型表現為低抗原提呈水平；C4 型的腫瘤免疫原性較低。

此外，剪接因子對于可變剪接的分型起到了主導作用。C2 型和C3 型腫瘤的主要剪接因子簇的表達明顯互斥。Mayeda等［22］證明了C3型剪接因子SRSF2和C2型剪接因子hnRNPA1 的拮抗關系，與本研究的結論相符。C1型的上述2 個剪接因子簇表達較為均衡，而C4 型主要以QKI、MBNL1、CELF2、RBFOX2 等 基 因 為 主。有 研究［23］發(fā)現MBNL1 和RBFOX2 的下調影響了體細胞重編程為多能干細胞過程中的剪接變化。另有研究［24］表明CELF2 通過在T 細胞激活過程中減少內含子保留事件，促進T 細胞分化和細胞因子產生。QKI 是重要的調控circRNA 生成的剪接因子［25］。有研究［26］顯示，QKI因子的下調與胃癌的預后有顯著關聯。根據剪接因子表達和抗原提呈的關系，我們推測C2 主導型剪接因子簇更容易引起異常剪接體的產生，而C3和C4主導型剪接因子簇通過無義介導的mRNA 降解等方式起到抑制異常剪接體通過翻譯發(fā)揮蛋白功能的作用。

綜上所述，本研究發(fā)現可變剪接事件的全局變化與胃癌的腫瘤生物學特性密切相關，且與腫瘤微環(huán)境的變化有潛在關聯。主要剪接因子的表達模式決定了可變剪接事件全局模式，且對于腫瘤的抗原提呈過程有著潛在的重要作用，對患者免疫治療的適用性起到了提示作用。同時，剪接因子網絡在胃癌中通過可變剪接調控適應和影響腫瘤免疫微環(huán)境的具體機制尚不明確，有待進一步研究。

參·考·文·獻

[1] Siegel RL,Miller KD,Jemal A. Cancer statistics,2019[J]. CA Cancer J Clin,2019,69(1):7-34.

[2] Mahoney KM, Freeman GJ, McDermott DF. The next immune-checkpoint inhibitors: pd-1/PD-L1 blockade in melanoma[J]. Clin Ther, 2015, 37(4):764-782.

[3] Bryan LJ, Gordon LI. Blocking tumor escape in hematologic malignancies:the anti-PD-1 strategy[J]. Blood Rev,2015,29(1):25-32.

[4] Kang YK, Boku N, Satoh T, et al. Nivolumab in patients with advanced gastric or gastro-oesophageal junction cancer refractory to, or intolerant of,at least two previous chemotherapy regimens (ONO-4538-12,ATTRACTION-2): a randomised, double-blind, placebo-controlled, phase 3 trial[J]. Lancet,2017,390(10111):2461-2471.

[5] Stamm S, Ben-Ari S, Rafalska I, et al. Function of alternative splicing [J].Gene,2005,344:1-20.

[6] Biamonti G, Catillo M, Pignataro D, et al. The alternative splicing side of cancer[J]. Semin Cell Dev Biol,2014,32:30-36.

[7] Kozlovski I, Siegfried Z, Amar-Schwartz A, et al. The role of RNA alternative splicing in regulating cancer metabolism[J]. Hum Genet, 2017,136(9):1113-1127.

[8] Ryan M, Wong WC, Brown R, et al. TCGASpliceSeq a compendium of alternative mRNA splicing in cancer[J]. Nucleic Acids Res, 2016, 44(d1):D1018-D1022.

[9] Ooi CH, Ivanova T, Wu J, et al. Oncogenic pathway combinations predict clinical prognosis in gastric cancer[J]. PLoS Genet,2009,5(10):e1000676.

[10] Josse J, Husson F. missMDA: a package for handling missing values in multivariate data analysis[J]. J Stat Softw,2016,70(1):1-31.

[11] Lê S, Josse J, Husson F. FactoMineR: an R package for multivariate analysis[J]. J Stat Softw,2008,25(1):1-18.

[12] Wilkerson MD,Hayes DN. ConsensusClusterPlus:a class discovery tool with confidence assessments and item tracking[J]. Bioinformatics, 2010, 26(12):1572-1573.

[13] Wu D, Smyth GK. Camera: a competitive gene set test accounting for intergene correlation[J]. Nucleic Acids Res,2012,40(17):e133.

[14] Liberzon A, Birger C, Thorvaldsdóttir H, et al. The molecular signatures database hallmark gene set collection [J]. Cell systems, 2015, 1(6): 417-425.

[15] H?nzelmann S,Castelo R,Guinney J. GSVA:gene set variation analysis for microarray and RNA-seq data[J]. BMC Bioinformatics,2013,14(1):7.

[16] Kahles A, Lehmann KV, Toussaint NC, et al. Comprehensive analysis of alternative splicing across tumors from 8,705 patients [J]. Cancer Cell,2018,34(2):211-224. e216.

[17] Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles[J]. NatMethods,2015,12(5):453-457.

[18] Thorsson V,Gibbs DL,Brown SD,et al. The immune landscape of cancer[J].Immunity,2018,48(4):812-830.e14.

[19] Leone P, Shin EC, Perosa F, et al. MHC class I antigen processing and presenting machinery: organization, function, and defects in tumor cells[J].J Natl Cancer Inst,2013,105(16):1172-1187.

[20] Oltean S, Bates DO. Hallmarks of alternative splicing in cancer[J].Oncogene,2014,33(46):5311-5318.

[21] Dunn GP, Bruce AT, Ikeda H, et al. Cancer immunoediting: from immunosurveillance to tumor escape[J]. Nat Immunol, 2002, 3(11): 991-998.

[22] Mayeda A, Krainer AR. Regulation of alternative pre-mRNA splicing by hnRNP A1 and splicing factor SF2[J]. Cell,1992,68(2):365-375.

[23] Venables JP,Lapasset L,Gadea G,et al. MBNL1 and RBFOX2 cooperate to establish a splicing programme involved in pluripotent stem cell differentiation[J]. Nat Commun,2013,4:2480.

[24] Mallory MJ,Allon SJ,Qiu J,et al. Induced transcription and stability of CELF2 mRNA drives widespread alternative splicing during T-cell signaling[J]. PNAS,2015,112(17):E2139-E2148.

[25] Conn SJ, Pillman KA, Toubia J, et al. The RNA binding protein quaking regulates formation of circRNAs[J]. Cell,2015,160(6):1125-1134.

[26] Bian Y, Wang L, Lu H, et al. Downregulation of tumor suppressor QKI in gastric cancer and its implication in cancer prognosis [J]. Biochem Biophys Res Commun,2012,422(1):187-193.