田園芳，陳 偉,2*

(1. 華北理工大學生命科學學院 河北 唐山 063000；2. 成都中醫(yī)藥大學中醫(yī)藥創(chuàng)新研究院 成都 611137)

RNA 可變剪接(alternative splicing, AS)是真核生物中剪接體選擇性剪切外顯子形成不同RNA 異構體的過程[1]，是調(diào)節(jié)基因表達、產(chǎn)生蛋白質(zhì)分子多樣性的關鍵環(huán)節(jié)[2]。文獻[3-4]發(fā)現(xiàn)，90%以上的人類基因都會經(jīng)歷可變剪接?？勺兗艚硬粌H增加了生物分子的復雜性、多樣性，還與疾病的發(fā)生有關[5]。如癌基因ETS 中外顯子7b 的剪接與細胞增殖的減少有關[6]。由于剪接因子MBNL3 的調(diào)節(jié)，lncRNA PXN-AS1 的外顯子4 被保留在轉錄本中，促進了肝癌的發(fā)生[7]。丙酮酸激酶前體mRNA 在剪接過程中保留了外顯子10，產(chǎn)生的亞型PKM2 過度表達，導致了腫瘤的發(fā)生[8]?？勺兗艚赢a(chǎn)生的雌激素受體α 和β 的變體ERα46 和ERβ1 則與乳腺癌密切相關[9]。因此，對可變剪接調(diào)控機制的研究尤為重要。

RNA 可變剪接并非獨立的生物過程，而是與轉錄過程存在著時空上的偶聯(lián)[10]。除了位于外顯子和內(nèi)含子中的順式和反式元件[11-12]，可變剪接還受到組蛋白修飾、DNA 甲基化等表觀遺傳因素的調(diào)節(jié)。隨著DNA 元件百科全書(encylopedia of DNA elements, ENCODE)計劃的深入開展，組蛋白修飾參與可變剪接調(diào)控的現(xiàn)象也被逐漸發(fā)現(xiàn)，尤其是出現(xiàn)在內(nèi)含子/外顯子區(qū)的組蛋白甲基化和組蛋白乙?；揎椗c細胞系特異性可變剪接密切相關[13]。組蛋白乙?；揎棇NA 可變剪接的調(diào)控在神經(jīng)細胞粘附分子(neural cell adhesion molecule, NCAM)基因中被發(fā)現(xiàn)。NCAM 基因中盒式外顯子的切除與外顯子中高含量的H3K9ac 密切相關[14]。NCAM基因盒式外顯子中H3K9ac 改變了該區(qū)域的染色質(zhì)結構，造成轉錄過程中RNA Pol II 的移動速率加快，從而導致NCAM 基因中盒式外顯子發(fā)生了可變剪接。在小鼠胚胎干細胞分化為神經(jīng)元的過程中，組蛋白乙?；揎椷€調(diào)控了Nf1 基因可變外顯子23a 和Fas 基因可變外顯子6 的可變剪接[15]。另外，成纖維原細胞生長因子受體基因(fibroblast growth factor receptor2, FGFR2)中 外 顯 子 區(qū) 的H3K36me3 是可變剪接調(diào)控蛋白的識別標記。FGFR2 基因中存在一對互斥外顯子IIIb 和IIIc，F(xiàn)GFR2-IIIb 只在上皮細胞中表達，而FGFR2-IIIc卻在間質(zhì)細胞中表達[16]。通過分析間質(zhì)細胞和上皮細胞中FGFR2 基因的組蛋白修飾后發(fā)現(xiàn)，與上皮細胞相比，間質(zhì)細胞的FGFR2 基因外顯子中富含H3K36me3。因此，染色質(zhì)重塑復合物MRG15 通過與FGFR2 基因中的H3K36me3 相互作用，能夠招募多聚嘧啶結合蛋白(polypyrimidine tract-binding protein, PTB)與FGFR2-IIIb 外顯子側翼的內(nèi)含子剪接抑制子結合，從而使得FGFR2-IIIb 在間質(zhì)細胞中被切除。

組蛋白修飾間還存在因果關聯(lián)，多種組蛋白修飾組合在一起形成級聯(lián)，共同調(diào)控基因表達[17-18]。在可變剪接過程中，不同類型的組蛋白修飾可以通過協(xié)同或拮抗方式調(diào)控剪接復合因子的招募，從而實現(xiàn)對RNA 剪接過程的調(diào)控。文獻[19]發(fā)現(xiàn)，人胚肺成纖維細胞系(IMR90 cell line)中BIN1 基因的可變剪接就是多種組蛋白修飾(H3K36me3、H3K4me3、H2BK12ac、H4K5ac)協(xié)同作用的結果。這些工作既為研究組蛋白甲基化和組蛋白乙酰化修飾調(diào)控可變剪接提供了理論依據(jù)，又顯示出從組蛋白修飾等表觀遺傳因素中挖掘新信息是認識可變剪接調(diào)控機制的新途徑。

外顯子跳躍模式是哺乳動物最常見的可變剪接模式[20-21]。文獻[22-23]發(fā)現(xiàn)了CD4+T 細胞外顯子跳躍模式中多種組蛋白修飾在外顯子和內(nèi)含子中富集程度的差異性，并利用組蛋白修飾差異信息對包含和排除外顯子進行了識別。通過構建組蛋白修飾間的相互作用網(wǎng)絡，文獻[24]還分析了CD4+T 細胞的外顯子跳躍模式中組蛋白修飾之間的因果關系。組蛋白修飾不僅可以通過改變RNA Pol II 的延伸率或招募剪接因子參與可變剪接的調(diào)節(jié)[25]，還能通過彼此之間的相互作用調(diào)控可變剪接。通過對人胚胎干細胞系(H1 cell line)的轉錄組學和表觀遺傳組學數(shù)據(jù)進行關聯(lián)分析，文獻[26]發(fā)現(xiàn)組蛋白修飾的動態(tài)變化與細胞特異性剪接機制相關。最近，文獻[27]發(fā)現(xiàn)了組蛋白修飾在H1 和IMR90 細胞系的外顯子跳躍模式中的協(xié)同分布規(guī)律，并在MCF10a、K562 和HeLa 等細胞系中進行了驗證，遺憾的是其并未闡明組蛋白修飾間的因果關系。

鑒于此，本文以IMR90 細胞系中外顯子跳躍剪接事件為研究對象，分析了28 種組蛋白修飾在外顯子跳躍模式排除和包含外顯子上的相關性，通過構建貝葉斯網(wǎng)絡推斷了組蛋白修飾間的因果關系。

1 材料與方法

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

人類基因組(GRCh37 版本) cDNA 序列和基因注釋文件來自Ensembl 數(shù)據(jù)庫(https://asia.ensembl.org/index.html)。IMR90 細胞系轉錄組測序數(shù)據(jù)(RNAseq)來自GEO 數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/, GSM2400222)。使用高通量數(shù)據(jù)過濾工具Trimmomatic[28](版本0.39-2)對RNA-seq 進行處理，參數(shù)設置為“PEILLUMINACLIP: TruSeq3-PE-2.fa:2:30:10:1:true LEADING:3 TRAILING:3 SLIDING WINDOW:4:20 MINLEN:40 TOPHRED33”。

從GEO 數(shù)據(jù)庫(登錄號GSE16256)獲取IMR90細胞系的28 種組蛋白修飾數(shù)據(jù)的bed 文件，相應的GEO 樣本號如表1 所示。

表1 組蛋白修飾的GEO 樣本號

1.2 轉錄本定量

統(tǒng)計RNA-seq 比對到基因組上的短讀序列(reads)的數(shù)量，即每個轉錄本的表達量。由于測序深度、轉錄本長度等因素的影響，表達量并不能直接代表轉錄本的相對豐度。因此對這些表達量進行TPM[29]標準化：

式中，Ni表示映射到第i個轉錄本的reads 數(shù)；Li表示第i個轉錄本的長度。

使用Salmon[30](版本1.5.1)計算轉錄本的TPM值，首先在參考cDNA 序列上建立索引，參數(shù)使用默認值。生成的索引與處理后的RNA-seq 進行定量，參數(shù)設置為“-l A–validateMappings–gcBias--seqBias”。

1.3 外顯子跳躍事件分析

外顯子跳躍事件中被保留在成熟轉錄本中的外顯子為包含外顯子，被剪接的外顯子為排除外顯子。根據(jù)發(fā)生外顯子跳躍事件的基因上所有轉錄本的TPM，計算每個外顯子的包含率(percent spliced in, PSI)：

式中，TPMi表示基因的第i個包含跳躍外顯子的轉錄本的TPM 值；S1表示包含跳躍外顯子的轉錄本集合；TPMj表示基因的第j個轉錄本的TPM值；S2表示不包含跳躍外顯子的轉錄本集合。PSI的值越大，表示外顯子被包含在最終轉錄本的概率越大，PSI 的值越小，外顯子被排除的概率越大。

使用Suppa 軟件[31](版本2.3)，根據(jù)基因注釋信息生成外顯子跳躍事件，并利用TPM 計算每個跳躍外顯子的PSI，所有參數(shù)使用默認值。最終獲得36 468 個外顯子跳躍事件，其外顯子PSI 值的范圍為0～1。PSI 為1 或0 時，表示外顯子被包含或被排除在基因的所有轉錄本中。為了適當擴充數(shù)據(jù)量，定義PSI>0.85 的外顯子為包含外顯子，PSI<0.15的外顯子為排除外顯子。由于外顯子過短無法匹配到組蛋白修飾數(shù)據(jù)，因此只保留長度在150～300 bp 之間的外顯子，最后獲得5 122 個包含外顯子和4 638 個排除外顯子。

1.4 相關性計算

根據(jù)bed 文件中組蛋白修飾在基因組上的位點信息，使用BEDTools[32](版本2.30.0)計算其在排除和包含外顯子上的富集程度(reads 數(shù))。然后根據(jù)reads 數(shù)，使用R 語言(版本4.1.2)Hmisc 包中的rcorr函數(shù)分別計算排除和包含外顯子上組蛋白修飾之間的皮爾遜相關系數(shù)。由此獲得28 種組蛋白修飾間的相關性系數(shù)矩陣，用corrplot 繪制相關性熱圖，并設置統(tǒng)計顯著性閾值p=0.05。

1.5 貝葉斯網(wǎng)絡構建

貝葉斯網(wǎng)絡是一種描述變量間因果關系的統(tǒng)計推理模型，其網(wǎng)絡拓撲結構是有向無環(huán)圖(directed acyclic graph, DAG)[33]。網(wǎng)絡中包括節(jié)點和邊，節(jié)點表示隨機變量，邊表示變量之間的條件依賴關系。所有邊為單向箭頭，箭頭指向的節(jié)點為子節(jié)點，箭頭另一端的節(jié)點為父節(jié)點。在貝葉斯網(wǎng)絡中，如果兩個節(jié)點之間不存在有向邊，則說明這兩個節(jié)點彼此條件獨立；如果兩個節(jié)點間存在有向邊，則說明這兩個節(jié)點間存在因果關系。

根據(jù)1.4 節(jié)的reads 數(shù)，對排除和包含外顯子上組蛋白修飾進行離散化，“1”表示外顯子上存在組蛋白修飾(reads>1)，“0”表示外顯子上不存在組蛋白修飾(reads=0)。為了得到穩(wěn)定的組蛋白修飾相互作用網(wǎng)絡，采用10 交叉檢驗法對所得網(wǎng)絡進行驗證。首先將組蛋白修飾數(shù)據(jù)平均分成10 份，其中9 份被用作訓練集，用于構建基本的貝葉斯網(wǎng)絡，1 份被用作測試集對網(wǎng)絡的穩(wěn)定性進行驗證。使用WinMine(https://www.microsoft.com/en-us/research/project/winmine-toolkit/)構建組蛋白修飾的貝葉斯網(wǎng)絡，最終得到10 個不同的網(wǎng)絡。如果連接節(jié)點的有向邊在10 個網(wǎng)絡中均存在，則將這些節(jié)點和邊保留，用以構建最終網(wǎng)絡。利用Cytoscape[34](版本3.8.2)軟件展示所得的網(wǎng)絡拓撲結構。網(wǎng)絡中節(jié)點表示組蛋白修飾，并根據(jù)1.4 節(jié)得到的相關系數(shù)對邊賦值。

2 結 果

2.1 組蛋白修飾間的相關性分析

計算組蛋白修飾之間的相關系數(shù)，通過繪制相關性熱圖，分析組蛋白修飾間的相關性，如圖1 和圖2 所示。圖中“×”表示組蛋白修飾間的相關性不顯著，色塊大小表示相關性強弱。結果表明，IMR90細胞系中排除和包含外顯子上大部分組蛋白修飾之間存在協(xié)同(r>0)或拮抗作用(r<0)。

圖1 排除外顯子上組蛋白修飾之間的相關性

圖2 包含外顯子上的組蛋白修飾之間的相關性

具體而言，組蛋白乙?；g表現(xiàn)出了正相關性，并且具有強正相關的組蛋白修飾存在于同一組蛋白上(H2B：H2BK120ac、H2BK12ac、H2BK15ac、H2BK20ac；H3：H3K14ac、H3K18ac、H3K23ac、H3K27ac；H4：H4K5ac、H4H8ac、H4K91ac)。在CD4+T 細胞中也發(fā)現(xiàn)了組蛋白H3 和H4 上相同的組合模式，但其正相關性弱于IMR90 細胞系[24]。

組蛋白甲基化之間同時表現(xiàn)出正相關性和負相關性，如H4K20me1 和H3K79me1、H3K79me2 之間正相關，而H3K27me3 和H3K36me3、H3K79me1、H3K79me2 之間負相關。在間質(zhì)細胞中存在相同的現(xiàn)象，H3K27me3 和H3K36me3 對基因FGFR2 的外顯子Ⅲb 的保留和剪切表現(xiàn)出了拮抗作用[35]。

排除和包含外顯子上組蛋白修飾間的相關性也存在差異。如H3K36me3 在包含外顯子上與10 種組蛋白修飾(H2BK120ac、H2BK15ac、H3K14ac、H3K18ac、H3K23ac、H3K27ac、H3K4ac、H3K56ac、H4K5ac、H4K91ac)之間正相關，如圖2 所示，而在排除外顯子中則未發(fā)現(xiàn)顯著的相關性，如圖1 所示。文獻[19, 36]研究發(fā)現(xiàn)，H3K36me3 富集在包含外顯子上，推測其可能與不同的組蛋白修飾形成組合模式，參與了RNA 可變剪接的調(diào)控。

2.2 組蛋白修飾之間的因果關系分析

通過構建IMR90 細胞中組蛋白修飾間的貝葉斯網(wǎng)絡，如圖3 和4 所示，對外顯子跳躍剪接事件中組蛋白修飾間的因果關系進行推斷，發(fā)現(xiàn)貝葉斯網(wǎng)絡中涉及的組蛋白修飾不僅包括同一組蛋白中同一氨基酸不同程度的修飾，還包括不同組蛋白中同一類型的組蛋白修飾，以及不同組蛋白中不同類型的組蛋白修飾。

在網(wǎng)絡拓撲結構圖中，只有子節(jié)點的組蛋白修飾用紫色標出；只有父節(jié)點的組蛋白修飾用綠色標出；既有父節(jié)點也有子節(jié)點的組蛋白修飾用藍色標出。

比較兩組網(wǎng)絡拓撲結構發(fā)現(xiàn)，排除和包含外顯子上組蛋白修飾之間存在18 種相同的因果關系，分別為：H2A.Z→H2AK9ac→H2BK5ac、H2AK9 ac→H3K9me1、 H2A.Z→H3K27me3、 H2A.Z→H3K9me1、 H2A.Z→H3K9me3、 H2AK5ac→H2BK12ac、 H2AK5ac→H3K14ac、 H2AK5ac→H2BK15ac、 H2AK5ac→H4K91ac、 H3K18ac→H3K14ac、 H3K18ac→H2BK120ac、 H3K18ac→H3K56ac、 H3K4me2→H3K56ac、 H4K8ac→H3K79me1→H3K36me3、 H4K8ac→H3K36me3、H4K8ac→H2BK15ac。

在排除外顯子中，27 種組蛋白修飾構成了71種因果關系，如圖3 所示。H3K4me3 直接或間接調(diào)控了剩余的26 種組蛋白修飾。8 種組蛋白修飾(H2BK120ac、 H2BK15ac、 H2BK5ac、 H3K14ac、H3K27me3、H3K4me1、H3K56ac、H3K9me1)可能直接與剪接因子相互作用，參與了可變剪接調(diào)控。在包含外顯子中，26 種組蛋白修飾構成了35 種因果關系，如圖4 所示，其中7 種組蛋白修飾(H2A.Z、H2AK5ac、H3K18ac、H3K4ac、H3K4me2、H4K5ac、H4K8ac)直接或間接調(diào)節(jié)剩余的19 種組蛋白修飾。16 種組蛋白修飾(H2BK120ac、H2BK12ac、H2BK15ac、 H2BK20ac、 H2BK5ac、 H3K14ac、H3K27ac、 H3K27me3、 H3K36me3、 H3K4me3、H3K56ac、 H3K9ac、 H3K9me1、 H3K9me3、H4K20me1、H4K91ac)可能直接與剪接因子相互作用，參與了可變剪接的調(diào)控。

圖3 排除外顯子上組蛋白修飾之間的貝葉斯網(wǎng)絡

圖4 包含外顯子上組蛋白修飾之間的貝葉斯網(wǎng)網(wǎng)絡

此外，排除和包含外顯子對應的網(wǎng)絡復雜程度也存在明顯差異。排除外顯子對應的網(wǎng)絡拓撲結構復雜，組蛋白修飾組合的調(diào)控路徑長，并且部分組蛋白修飾間表現(xiàn)出了拮抗關系，如H4K8ac→H3K36me3→H3K27me3、 H3K9ac→H3K27me3，如圖3 所示。包含外顯子對應的網(wǎng)絡拓撲結構則相對簡單，組蛋白修飾組合的調(diào)控路徑較短，如圖4所示。另外IMR90 細胞中排除和包含外顯子上的H3K27ac 和H4K5ac 之間因果關系相反。

組蛋白修飾在可變剪接過程中的調(diào)控作用已被發(fā)現(xiàn)，如H3K4me3 與U2 snRNP 結合可調(diào)節(jié)剪接速率[37-38]。IMR90 細胞系包含外顯子上的調(diào)控網(wǎng)絡顯示H3K4me3 可能受到了H3K4ac 的影響，如圖4所示。而排除外顯子上H3K4me3 不僅與9 種組蛋白修飾(H3K4me2、H4K5ac、H3K4me1、H3K9ac、H3K18ac、 H2AK5ac、 H3K56ac、 H2BK15ac、H3K27ac)之間存在直接因果關系，還存在于所有的調(diào)控路徑中，如圖3 所示。在間質(zhì)細胞中，H3K36的去甲基化酶KDM2a 被募集到富含H3K27me3 的區(qū)域，保持了低H3K36me3 水平，從而促進了外顯子Ⅲb 的包含[35]。這一結果表明，H3K36me3 與H3K27me3 之間的拮抗作用調(diào)控了可變剪接。有意思的是，在IMR90 細胞系的排除外顯子中也發(fā)現(xiàn)了H3K36me3 和H3K27me3 之間的拮抗關系，如圖3 所示。由此推測，IMR90 細胞系中的可變剪接受到了組蛋白修飾間相互作用的調(diào)控，相關結果還需進一步實驗驗證。

3 結 束 語

本文通過對IMR90 細胞系中組蛋白修飾間的相關性進行分析，發(fā)現(xiàn)了外顯子跳躍剪接事件中組蛋白修飾間存在明顯的組合模式。通過構建貝葉斯網(wǎng)絡，分析了排除和包含外顯子中組蛋白修飾間的因果關系。由此推測，IMR90 細胞系排除和包含外顯子的組蛋白修飾可能通過組合方式與剪接因子相互作用，直接或間接地參與了可變剪接的調(diào)控。