【Paper reading】可变剪接预测&ENCODEC数据集

可变剪接预测&ENCODEC数据集

1. Epigenome-based splicing prediction using a recurrent neural network

2. An integrative ENCODE resource for cancer genomics

---

前言

Encode是一项公共联合研究项目，目标是定义人类基因组中所有发挥“功能”的DNA元件，并研究这些元件与蛋白和RNA之间的相互作用如何调控细胞的转录活动。Encode是迄今最详细的人类基因组分析数据，是对人类生命科学的又一重大贡献。

---

以下两篇论文都是基于Encode项目，都为耶鲁大学Mark Gerstein团队做的工作。

Epigenome-based splicing prediction using a recurrent neural network

  这篇论文使用递归神经网络，基于表观基因组进行剪接预测。文章发表在PLOS COMPUTATIONAL BIOLOGY，这个期刊是计算生信方向的老杂志，口碑很不错的。论文的通讯作者是耶鲁大学的Mark Gerstein。文章的code已经放在的github上。

1.Introduction

  首先来介绍一下选择性剪接：真核生物合成蛋白质的转录过程以特定的单链DNA片段作为模板，RNA聚合酶作为催化剂，合成前mRNA，前mRNA经进一步加工后转为成熟mRNA。mRNA又在翻译过程中作为合成蛋白质的模板。剪接（splicing），是一种基因重组现象，是真核生物的前mRNA变成mRNA的过程之一。选择性剪接是剪接过程可以透过对同一个基因转录的相同前mRNA使用不同的剪接选择，产生不同的mRNA异构物。前mRNA的可变剪接在转录组多样化起着不可或缺的作用。能够影响约95%的人类蛋白质编码。

  这篇文章提出了一种新的计算方法来表征表观遗传修饰在共转录剪接过程中的作用。实验采用了类似于共转录剪接的物理特征的递归神经网络（RNN）架构。将方法应用于ENCODE中的49种细胞类型以研究剪接调控，发现不仅空间，转录过程中染色质上下文时间方向的不对称同样影响转录过程中的剪接调控。

2.Methods

  我们来看一下文章提出的模型。图1是文章的模型结构。A图是生物背景，B图是模型输入，包括DNA上下文信息和染色质标记信息。C图是RNN模型，经过训练，模型能够输出外显子是否被包含。模型有三个优势。1.该模型可以从表观遗传环境的时间序列中学习，类似于在共转录剪接过程中，当Pol II沿着引导DNA链向前移动时，表观遗传环境逐渐变化。2.RNN模型使我们能够基于相邻的DNA序列和表观遗传修饰来预测外显子的包含。3.模型具有可解释性：模型的物理相似性使我们能够解释训练后的模型权重参数，并探索指导DNA元素与周围表观遗传修饰之间的时空联系。

  提出的模型叫ESPRNN模型。使用基因组序列和表观基因组环境作为输入，模型由两个循环层组成。来自3’和5‘剪接位点的输入在第一个循环层中分别进行处理，然后在下一个循环层中合并。softmax分类器用于确定外显子的包含。我们的模型使用五个核心组蛋白修饰轨迹，在整个细胞类型中，基于LSTM的模型的平均F1得分（精度和召回率的谐和平均值）为0.8472（基于GRU的模型为0.8757）。（B）来自六种不同ENCODE细胞类型的Precision-recall曲线。（C）通过留一法分析确定剪接过程中各元素重要性的顺序，是筛选出的有助于剪接调控的表观遗传特征。（D）LSTM模型与其他模型的比较，发现LSTM模型在分类准确度方面优于其他模型。实验发现H3K36me3是预测剪接的最重要变量。

  H3K36me3是一种组织蛋白。先来简单介绍一下组织蛋白。组织蛋白分为两大类，五小类。组织蛋白修饰的命名是：

1. 先以组织蛋白名称开始，如H3；
2. 单一字母的氨基酸简称，如K代表赖氨酸，及在蛋白质的位置；
3. 修饰的种类，Me即甲基化、P即磷酸化、Ac即乙酰化及Ub即泛素化。举例来说，H3K4Me就代表组织蛋白H3从N端开始起计第4个赖氨酸的甲基化。

  那什么是3’和5’端呢：内含子总是在两端具有两个不同的核苷酸。 在5’端，DNA核苷酸是GT [（pre-mRNA）中的GU]； 在3’端他们是AG。 这些核苷酸是剪接位点的一部分。

  总体实验发现H3K36me3是预测剪接的最重要变量。有趣的是，剪接的最高预测因子之一是H3K79me2，该基因先前已显示与基因体上的H3K36me3相关。H3K9me3是一种组蛋白修饰，它可以募集衔接子蛋白（如HP1）以促进剪接因子，也被列为最佳预测因子。当我们分析模型时，发现循环单元中各个门的训练权重也表明，开放染色质（DHS），H3K27ac，K3K36me3和H3K4me1的权重比其他表观遗传特征的权重更高（如预期的那样）。同时还注意到，H3K27me3和K3K9me3在输入门处负加权，这表明这些功能对外显子包涵有负面影响。

3.Results

  文章还做了一些其他实验。实验通过检查外显子-内含子边界周围多个组蛋白修饰和DNA甲基化的富集，研究选择性剪接的表观遗传背景。
研究者在ENCODE中汇总了49种细胞类型的多个组蛋白修饰，观察到它们的富集程度是与剪切位点距离的函数。发现3‘受体和5’供体在剪切位点的100 碱基对内最有富集的趋势。组织蛋白H3K36me3和H3K27me3的强富集模式出现在外显子边界附近从3’受体开始，在进入外显子约100 bp处出现了峰富集。在5’供体处，它离外显子更近，大约50 bp。
可选的剪接调控是一个复杂的过程，需要多个剪接因子和酶的精确协调。研究表明，RNA结合蛋白（RBP）有助于转录过程中的剪接调控。图C为跨外显子-内含子边界的RBP富集。由图C可以发现剪接因子的空间结合偏好可能显示出略有不同的模式。特别是，RNA结合蛋白hnRNP A1和SRSF1在3‘SS以外的内含子区域富集，而SF3B4和hnRNP C在外显子区域富集。

  实验还测试了组蛋白修饰是否对替代外显子的包含和表达有任何影响。观察图A和B，外显子-内含子边界处H3K36me3的富集与外显子表达呈正相关，而H3K27me3呈负相关。图C可以看出H3K36me3和H3K27me3具有统计学显着富集和耗竭。这些组蛋白甲基化与外显子表达的对比趋势和相关性表明，剪接编码可以在表观基因组环境中直接或间接编码。
为了探索表观遗传调控与组织特异性剪接之间的潜在联系，实验还研究了49个ENCODE生物样品中的剪接模式。基于所有编码外显子的剪接百分比（PSI）值的相似性，使用层次聚类将生物样本分为五个类别，如图D。来自同一起源细胞的组织类型之间的剪接模式高度相关，基于表观遗传标记再现相似的聚类结果血统细胞类型形成簇C2，而脑和神经细胞簇在簇C4中。 此外，我们观察到癌细胞系在群集C3中聚集在一起。

总结

  最后来总结一下这篇文章，预测模型表明，剪切位点的表观基因组特征在确定剪接结果中起着重要作用。此外，积极的结果表明我们的模型可以扩展以从基因组和表观基因组的背景预测完整的转录组组成。未来希望可以通过使用ENCODE项目中的全部可用表观基因组数据来增加更多的深层隐藏层并增加训练样本的数量来进一步改进模型。我们的方法确实存在一些局限性，因为在经过训练的模型中可视化和评估权重矩阵的多维上下文仍然具有挑战性。我们可以应用降维技术来探究各种表观基因组信号之间关系的潜在表示。未来的研究应集中于比较正常样本和癌症样本的剪接模型，以期阐明剪接调控的表观基因组情况的差异。

An integrative ENCODE resource for cancer genomics

  这篇论文介绍了用于癌症基因组学的集成式ENCODE资源，论文今年发表在nature communications，是个非常新的杂志，基本上每期都会有一些生信类的文章，期刊影响因子在不断上涨，很有潜力。文章的通讯作者还是耶鲁大学的Mark Gerstein，和上一篇文章的是同一个团队做的。Supplementary Data的网址。

1.Introduction

  ENCODE包含数千个功能基因组学数据集，数百种细胞类型，为基因组解释提供了通用注释。然而，对于特定的应用，使用定制的注释可能更为有效。研究者通过利用多种数据丰富的ENCODE细胞类型上的高级检测方法开发了一种自定义注释。
文章的主要工作包括五个方面。

1. 首先开发了一种自定义注释，注释的关键是转录因子和RNA结合蛋白（TF和RBP）的综合性和实验性衍生网络。癌症作为一种系统失调的疾病，是这种基于网络注释的理想应用。
2. 具体来说，对于与癌症相关的细胞类型，研究者将regulator置于层次结构网络中，并测量其在肿瘤发生过程中的网络变化（重连线）。
3. 研究者还对TF-RBP的串扰进行了广泛的调查，强调了以前未表征的RBP-SUB1如何驱动异常的肿瘤表达并增强著名的致癌转录子–MYC的作用。
4. 此外，研究了肿瘤到正常细胞状态的轨迹，许多normal-to-tumor转变都趋向于茎样状态，而致癌基因的敲低则显示出相反的趋势。
5. 最后，文章将资源组织到一个统一的工作流程中，除regulator，对key elements和variants,进行优先级排序，展示排序的应用，并通过靶向实验验证。

  五个主要工作分别用五张图来展示，我后面会分别介绍这五张图。

  首先来介绍文章提出的Encodec资源。ENCODEC是针对ENCODE专门的癌症基因组学数据。首先，使用ENCODE数据，针对每种癌症，尝试找到可用的最佳肿瘤-正常配对。然后，构建了派生资源。包括（1）全面的网络，使我们能够看到网络重新布线和regulator层级的全局变化；（2）带注释的细胞类型目录，使我们能够相对于正常细胞和干细胞进行致癌性改变；（3）紧凑的非编码注释和扩展的基因定义，可能会增加统计能力来解释基因组变异（种系和体细胞）和基因表达变化。

2.Overview of the ENCODEC resources

  图1是Encodec资源的综述。图中有两个关键维度：跨细胞类型的宽度和跨测定的深度。横坐标表示229种细胞类型，纵坐标表示Encode数据不同类型的测定方法。表格构建了具有两个关键特征的深层集成注释：（1）紧凑地定义非编码元素以更精确地定位功能位点；（2）这些不连续的调控区域是与基因链接形成扩展基因定义。注释在两个方面都很紧凑：它包含的总元素较少（因为跨许多测定的深度整合消除了许多潜在的假阳性），每个元素的长度往往较短，但在功能上相关的核苷酸却更加丰富。可以通过减少多个测试负担或通过删除每个元素中的非功能性核苷酸来更清晰地定义核心区域，从而提高统计能力。
  插图1表示使用特定于细胞类型的ENCODE分析法中的注释来构建扩展的基因定义：编码和非编码元素根据它们的相互作用和相关功能而链接。插图2关联转录因子（TFs）和RNA结合蛋白（RBPs）形成层次网络，描述regulator的调节潜力。插图3 通过比较肿瘤和正常ENCODE样本中的调节网络，开发了可能与正常到肿瘤过渡情况下发生的调节变化有关的重新布线网络。

3.Levaraging ENCODE networks to prioritize regulators

  我们来具体看图2，regulator网络的层次结构。构建多模态TF-RBP网络后，我们将其系统地分为一个层次结构（图 2a，b）。a图为TF，b图为RBP，分层标准是目标表达的相关性，联合TF-RBP形成regulator层次网络。网络中较高层的元素倾向于调节较低层的元素。在层次结构中，我们发现顶层TF和RBP更显着地驱动差异表达（p值<2.2e-16，单面Wilcoxon检验）。联合TF-RBP网络还可以研究TF和RBP之间的交叉调节。有趣的是，实验发现与上层和中层相比，下层的TF-RBP交互作用更少。此外，我们注意到著名的癌基因MYC是位于层次结构顶层的主要TF之一。MYC不仅直接调节其他TF的表达，而且还靶向许多RBP。
  网络还可以在肿瘤样品中进行基因表达分析。研究者使用了基于回归的方法来系统搜索在不同癌症中最能驱动肿瘤正常差异表达的TF和RBP，对于每位患者，我们测试了调节剂的活性与其靶标的肿瘤至正常表达变化之间的相关程度。然后，我们计算了每种癌症类型中具有这些关系的患者百分比，并在图2c中显示了TF和RBP的总体趋势，红色和蓝色分别表示上调和下调。从图中可以看出MYC的靶基因在许多癌症类型中均显著上调。图 2d为MYC调节活性升高与乳腺癌的疾病特异性生存率（DSS）降低有关。使用敲除法验证MYC的调节作用。与预测一致，在MCF-7中MYC敲低后，MYC靶标的表达显着降低。
   研究者以类似于TF网络的方式分析了RBP网络，找到了与每种癌症相关的regulator。例如，RBP SUB1的ENCODE eCLIP配置文件的结合峰富集在基因的3‘UTR区域，SUB1的预测靶标在许多癌症类型中均显着上调（图 2f，左）。作为RBP，SUB1以前与癌症无关，因此研究者试图研究其作用。敲除HepG2细胞中的SUB1可以显着下调其靶标，并且SUB1靶标的衰减率低于非靶标的衰减率（图 2f，右）。此外，我们发现SUB1靶标的上调可能导致某些癌症类型患者的生存率下降。
  然后，研究者使用网络来研究优先regulator如何与彼此和其他基因相互作用。对于TF，首先研究了第二个TF如何共同调控MYC的靶基因。在图2e中示出了对所有可能的三方共同监管关系的解释。我们发现最常见的模式是特征明确的前馈循环（FFL）。在这种情况下，MYC调节另一个TF，并且共同调节targets。关于RBP，发现SUB1的最高共同监管伙伴实际上是MYC。SUB1是许多细胞类型中MYC的直接的target。

4.TF-Gene network rewiring

  图三为TF-基因的重布线网络。如图3a所示，研究者还将细胞类型的网络组织到了层次结构中 。根据图3a可以看出边缘优势最大的基因也在监管体系网络的重新布线上处于血液癌症这些层次结构的最高层。另外，我们发现与驱动癌症基因表达变化更相关的TF也往往是最重要的。MYC是高度重新布线的TF和一种驾驶表达方式的最突出示例。相反，受突变影响更深的TF位于层次结构的底部。图3b中，JUND是CML的最高优势者，其靶标显示基因表达增加。但是，其结合位点几乎不受SV或SNV的影响。该结果表明结合位点突变不会驱动调节变化。我们发现大多数重新接线事件与显着的基因表达和染色质状态变化有关，但不一定与直接的变体诱导的基序丧失或获得事件有关。
  除了通过上面的表达分析进行的TF调节活性改变外，研究者还直接测量了肿瘤正常对的调节边缘改变的分数，以研究TF靶点在肿瘤发生中的变化。研究者称其为重布线索引并根据分数对TF进行排序（图3c）。如图 3c，在白血病中，最著名的癌基因（如MYC和NRF1）是边缘优势最大的基因，而著名的肿瘤抑制因子IKZF1是最显着的优势基因。使用远端，近端和组合网络观察到了类似的重新布线趋势（图 3c）。

5.Oncogenic transformation and cell state

  研究者将来自299种ENCODE细胞类型的RNA-seq数据投影到一个低维空间中，发现各种类型的干细胞形成紧密的簇，如图 4。此外，存在一种趋势，其中从正常细胞到肿瘤细胞的轨迹涉及沿着单个干样成分向干细胞移动。这对于多种不同的癌症都是如此。该观察结果与以前使用表达和甲基化分析的努力一致。值得注意的是，我们从近端和远端染色质数据观察到一致（甚至更强）的模式，可以将其视为观察到的基因表达变化的根本原因。
  我们将ENCODE细胞类型的表达谱（左，poly-A长RNA-seq），近端网络（右数，CTCF ChIP-seq）和远端网络（右边，候选顺式调控元件）投影到较低的位置维度空间。干细胞样细胞类型形成簇，表明干细胞样细胞类型与正常细胞和癌性细胞类型截然不同。此外，我们发现癌细胞类型倾向于更靠近茎样簇。K562中的致癌基因敲除导致与正常细胞类型的转录组相似性更高，而肿瘤抑制基因（TSG）敲除导致与肿瘤细胞类型的相似性更高（与GM12878相比，从左上角第二个）。通常，我们发现致癌基因敲低会导致沿茎状成分向正常状态的轻微回复（左侧，底部第二个）。

6.Variant prioritization and validation

  如图5所示 ，ENCODEC资源组织到一个统一的工作流程中，启用了逐步优先级划分，使我们能够精确定位与肿瘤发生有关的关键调控因子，非编码元件和变异体。具体来说，我们首先重点介绍regulator，这些调节器要么重布线，要么位于集线器中，要么位于层次结构的顶部，要么大大推动了癌症表达的改变。然后，我们优先考虑与这些调节子相关的功能元件，这些功能元件要么承受突变的重担，要么发生大量的染色质变化，要么发生扩展的基因连锁变化。最后，在核苷酸水平上，我们通过估计SNV破坏或引入特定结合位点的能力并评估它们位于优先元件中的程度来确定其优先级。

总结

  最后来总结一下这篇文章。文章描述了一个自定义的ENCODE注释：一种提供包括扩展基因在内的集成网络注释的伴随资源。癌症基因组学是突出资源价值的理想应用，我们将展示它如何根据细胞空间轨迹和网络重新布线来帮助描述致癌性转化。我们还使用专用注释对关键的regulator，元素和变体进行优先级排序。