外顯子跳躍模式中組蛋白修飾的組合模式分析

田園芳，陳 偉,2*

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

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

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

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

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

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

1 材料與方法

1.1 數據來源

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

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

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

1.2 轉錄本定量

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

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

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

1.3 外顯子跳躍事件分析

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

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

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

1.4 相關性計算

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

1.5 貝葉斯網絡構建

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

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

2 結 果

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

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

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

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

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

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

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

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

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

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

比較兩組網絡拓撲結構發現，排除和包含外顯子上組蛋白修飾之間存在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 直接或間接調控了剩余的26 種組蛋白修飾。8 種組蛋白修飾(H2BK120ac、 H2BK15ac、 H2BK5ac、 H3K14ac、H3K27me3、H3K4me1、H3K56ac、H3K9me1)可能直接與剪接因子相互作用，參與了可變剪接調控。在包含外顯子中，26 種組蛋白修飾構成了35 種因果關系，如圖4 所示，其中7 種組蛋白修飾(H2A.Z、H2AK5ac、H3K18ac、H3K4ac、H3K4me2、H4K5ac、H4K8ac)直接或間接調節剩余的19 種組蛋白修飾。16 種組蛋白修飾(H2BK120ac、H2BK12ac、H2BK15ac、 H2BK20ac、 H2BK5ac、 H3K14ac、H3K27ac、 H3K27me3、 H3K36me3、 H3K4me3、H3K56ac、 H3K9ac、 H3K9me1、 H3K9me3、H4K20me1、H4K91ac)可能直接與剪接因子相互作用，參與了可變剪接的調控。

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

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

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

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

3 結 束 語

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