外显子跳跃模式中组蛋白修饰的组合模式分析

田园芳，陈 伟,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 数据来源

人类基因组(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.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 细胞系中的可变剪接受到了组蛋白修饰间相互作用的调控，相关结果还需进一步实验验证。

本文通过对IMR90 细胞系中组蛋白修饰间的相关性进行分析，发现了外显子跳跃剪接事件中组蛋白修饰间存在明显的组合模式。通过构建贝叶斯网络，分析了排除和包含外显子中组蛋白修饰间的因果关系。由此推测，IMR90 细胞系排除和包含外显子的组蛋白修饰可能通过组合方式与剪接因子相互作用，直接或间接地参与了可变剪接的调控。