利用转录组测序数据分析可变剪接的方法

利用转录组测序数据分析可变剪接的方法作者：***来源：《科学与信息化》2020年第08期摘要可变剪接是调节基因表达和产生蛋白组多样性的重要因素，同时参与调控细胞分裂、分化及凋亡等重要生物学过程，异常的可变剪接多与人类疾病有关。

随着新一代测序技术和生物信息学的快速发展，以及先进计算方法的提出，使得我们对可变剪接有了深入的认识。

并且基于剪接机制对于病的靶向药物设计，已得到了有效的临床治疗效果。

本文主要阐述了近年来基于二代测序技术开发的几种识别可变剪接的计算方法，并对未来的发展方向进行展望。

关键词可变剪接;二代测序技术;生物信息学;分析工具可变剪接，又称选择性剪接（Alternative Splicing，AS），是真核生物基因表达的普遍调节机制，是指一个前体mRNA经过不同的剪接形式产生多种不同剪接异构体的过程。

在1978年，Walter Gilbert提出了内含子和外显子命名[1]，不同外显子组合产生特异的异构体。

二代测序技术的迅速发展极大地推动了人类对可变剪接的认识。

现有数据表明，人类大约有92%-94%的基因都会经历某种程度的可变剪接行为，并且在20000多种人类蛋白编码基因中，约37%的基因会编码产生不同的蛋白亚型，这表明可变剪接增加了蛋白质组的多样性和复杂性[2]。

AS对基因的功能起着重要调控作用，同一基因的不同亚型可能参与不同的生物学过程。

例如p53抑癌基因（TP53）在DNA受损细胞的调控中起着核心作用，然而其Δ133β亚型则可以抑制全长p53β亚型5和6从而诱导肿瘤细胞的凋亡[3]。

另外AS几乎参与了所有生物学过程，包括调节细胞的分裂和凋亡、神经系统的发育以及细胞对抗多种环境因素做出的免疫应激反应等[4]。

另一方面，AS的异常调节还与多种遗传性疾病和恶性肿瘤相关，包括神经退行性疾病、心血管疾病和代谢状况等。

据报道，与SNP相关的遗传性疾病多达一半是由于剪接受损引起的[5]。

AS的异常调节对癌症的发生发展有重要的作用，为疾病的发展提供了可能的新颖治疗靶标和生物标志物的来源，而AS位点的预测可以为药物设计提供很好的分子基础。

本文主要阐述了最新的基于二代测序数据计算可变剪接方法的基本原理及应用，并对未来的研究方向做了初步的展望。

1 可变剪接的主要类型目前大量数据研究结果表明，可变剪接主要包括五种形式（图1），分别为外显子跳跃（Skipped Exon，SE），可变5′剪接位点（Alternative 5′ Splice Site，A5SS），可变3′剪接位点（Alternative 3′ Splice Site，A3SS），互斥外显子（Mutually Exclusive Exons，MXE）和内含子保留（Retained Intron，RI）。

此外还有两种不常见的形式：可變的第一个外显子（Alternative First Exon，AFE）和可变的最后一个外显子（Alternative Last Exon，ALE）。

其中SE是人类基因组中最主要的可变剪接模式（约35%），其次是A3SS（16%）与AA5SS （15%），RI发生最少，仅占到全部事件的1%[6]。

2 分析可变剪接的软件可变剪接现象从70年代发现后，其基本的科学问题聚焦于剪接位点的发现、差异分析、调控元件和网络的发现和构建等。

RNA-seq技术的发明，使系统、定量的可变差异剪接分析成为可能。

大量测序数据的可变剪接差异分析需要优秀的统计模型和计算工具，因此一直是生物信息学领域的重要研究课题之一。

2.1 转录组读段回帖序列比对软件STAR（Spliced Transcripts Alignments to a Reference，STAR）[7]是专门针对RNA-seq数据开发的，该算法使用未压缩后缀阵列中的连续最大可比种子进行搜索，接着进行种子聚类（Seed clustering）和缝合过程（Stitching procedure）。

除了典型剪接的无偏从头检测外，STAR还能发现非典型拼接和嵌合（融合）转录本，并且能够比对全长RNA序列。

STAR是目前转录组测序分析最常用的软件之一，由C++编写，运行速度快，同时也是免费开源软件。

TopHat[8]是基于Bowtie[9]短序列比对软件开发的，它可以实现对短序列的快速比对，根据比对的结果可以检测到外显子间的剪接位点。

该算法基本思想是对Bowtie比对结果中大量重叠的片段寻找标志性的剪接位点，并列出它们所有可能的组合情况。

此外对于Bowtie没有比对到基因组上的片段，TopHat会重新建立索引进行比对，并与上个过程得到的所有可能组合进行对比，最终将跨越剪接位点的片段比对到基因组。

2.2 比对后的可变剪接分析rMATS（ replicate multivariate analysis of transcript splicing）是目前最常用的分析RNA-seq 数据可变剪接事件的软件[10]。

该算法依赖于有重复的数据集，利用分层模型同时考虑重复样本的不确定性和变异性。

以SE为例，rMATS首先采用Exon Inclusion level 来定义样本中可变剪接事件的表达量，正常的isoform称为Exon Inclusion Isoform，发生了外显子跳跃的转录本称之为Exon Skipping Isoform，比对到Inclusion; isoform上的reads用L表示，比对到Skipping Isoform上的reads用S表示（图2），则该外显子跳跃的定量如下：Exon Inclusion level为Inclusion isoform所占的比例，计算时用长度（LI和Ls）校正了原始的reads数。

给定阈值c，; 判断两个样本中相对应的inclusion level是否发生了变化：c，c用户可以自定义，取值范围为0～1，表示两个样本中Inclusion level的差值。

如c=0.1，则表示两个样本中该可变剪接事件的Inclusion level相差10%。

其他类型的可变剪接事件也可以划分成上述两种isoforms，进行差异分析比较的是两组样本中inclusion level的差异。

Exon Skipping Isoform為Upstream exon 和Downstream exon直接连接形成;Exon Inclusion Isoform为Upstream exon， Alternative exon和Downstream exon连接形成。

该模型以Likelihood-ratio test计算p值，大大提升了计算速度。

rMATS支持多线程运行且支持两种输入格式：Fastq或者Bam。

根据计算时用到的reads差别，最后会得到两组结果，一种是只用到跨Junction的reads;另一种是比对到剪接位点上的所有reads。

rMATS是目前在RNA-seq数据领域应用最多的分析可变剪接的工具。

DARTS[11]是2019年3月由邢毅教授团队在开发rMATS工具后，新提出的一个基于深度学习检测可变剪接的计算框架，这也是首次将深度学习与贝叶斯假设检验相结合，用于RNA 的可变剪接分析。

DARTS由两部分构成（图3）：深度神经网络（DNN）和贝叶斯假设检验（BHT）。

DNN模型是基于顺式序列特征和RNA结合蛋白（RBP）的mRNA表达水平预测差异剪接的结果;BHT统计模型通过将经验证据整合到特定的RNA-seq数据集中，并结合先验的差异可变剪接概率来进行预测。

在训练过程中，DARTS BHT使用无先验概率信息（DARTS BHT（flat），仅将RNA-seq数据用于推断）分析大规模RNA-seq数据，以生成高置信度差异及未发生剪接的训练标签之间的事件，以此作为DARTS DNN的训练输入。

在应用过程中，已训练完成的DARTS DNN用于预测用户特定数据集中的差异可变剪接，然后与DARTS BHT （DARTS BHT（info））观察到的RNA-seq片段的计数结合到一起，以达到深度学习增强的剪接分析。

DARTS BHT（flat）进行常规分析大规模RNA-seq数据中的可变剪接事件，创建带标签的训练数据，用于训练 DNN模型;新的特定RNA-seq经DNN模型预测作为贝叶斯模型的先验（DARTS BHT（info））;用户的RNA-seq数据则是用于更新先验概率形成后验概率。

顺式序列特征（Cis-sequence）和反式RBP的mRNA水平（Trans-RBP）：DARTS DNN预测差异可变剪接的两个因素。

先验信息（Prior）：DARTS DNN预测的结果。

可能性（Likelihood）：特定RNA-seq数据集观察到的读数。

差异剪接的后验概率（Posterior）由Cis-sequence和Trans-RBP生成的先验概率和Likelihood组成。

DARTS的深度神经网络部分整合了剪接位置附近的序列信息、进化信息、可变剪接产生的RNA二级结构信息和RBP的表达信息等，提高了该模型的准确性。

相较于之前的例如MISO[12]、SpliceTrap[13]、SUPPA[14]等软件，DARTS对测序深度的要求降低。

总的来说，DARTS综合了深度学习和贝叶斯假设检验统计模型的优点，为低测序深度的数据提供了更好的分析可变剪接的方法，拓展了传统RNA-seq可变剪接分析的敏感度和准确度。

3 总结与展望可变剪接（AS）是基因在转录后水平的一项复杂调控机制。

在最近的几十年中，我们对AS的机理和动力学的了解呈指数级增长，人们也越来越认识到AS在人类健康和疾病中发挥的关键作用。

随着RNA-seq和各种计算工具的不断涌现，使得我们可以在全基因组水平检测AS现象。

同时随着近年来单细胞测序技术的发展，可能会在单个细胞中出现新的剪接现象，这需要我们去开发专门针对单细胞测序的更优计算方法。

另外机器学习方法将会加快我们对可变剪接的认识[15]。

最后，可变剪接变体在功能上有助于疾病的耐药性，因此靶向可变剪接可作为某些疾病的新型疗法。

目前在人类某些疾病（如SMA）中已取得了很好的临床治疗效果[16]，但在癌症方面还未有重大突破，相信随着研究的不断深入，在不久的将来也能为癌症提供有效的可选药物。

参考文献[1] Gilbert，Walter. Why genes in pieces？[J]. Nature，1978， 271（5645）：501.[2] Min-Sik， K. A draft map of the human proteome.[J].Nature， 2014，509（7502）：575-581.[3] Aoubala M，Murray-Zmijewski F，Khoury M P，et al. P53 directly transactivatesΔ133p53α， regulating cell fate outcome in response to DNA damage[J]. Cell death and differentiation，2011，18（2）：248-258.[4] Ramanouskaya T V，Grinev V V . The determinants of alternative RNA splicing in human cells[J]. Molecular Genetics & Genomics，2017，（3）：263.[5] Garcia-Blanco M A，Baraniak A P，Lasda E L . Alternative splicing in disease and therapy[J]. Nature Biotechnology，2004，22（5）：535-546.[6] Wang E T，Sandberg R，Luo S，et al. Alternative isoform regulation in human tissue transcriptomes[J]. Nature，2008，456（7221）：470-476.[7] Alexander D ， Davis C A ， Felix S ， et al. STAR： ultrafast universal RNA-seqaligner[J]. Bioinformatics，，2013（1）：1.[8] Trapnell C，Pachter L，Salzberg S L . TopHat： discovering splice junctions with RNA-Seq[J]. Bioinformatics，2009，25（9）：1105-1111.[9] Langmead B，Trapnell C，Pop M，et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[J]. Genome Biology，2009，10（3）：1-10.[10] Shihao Shen，Juw Won Park，Zhixiang Lu，et al. rMATS： Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data[J]. proc natl acad sci u s a，2014，111（51）：5593-601.[11] 佚名. Deep-learning augmented RNA-seq analysis of transcript splicing[J]. Nature Methods，2019，（1）：157.[12] Yang-Ming S I，Yong-Qiang X，Lu C . Differential splicing event analysis of liver tumor-educated blood platelets RNA-seq data with Hisat2 and MISO[J]. journal of inner mongolia university of science and technology，2016，（3）：73.[13] Zhang M Q . SpliceTrap： a method to quantify alternative splicing under single cellular conditions[J]. Bioinformatics，2011，27（21）：3010-3016.[14] more samples. SUPPA： a super-fast pipeline for alternative splicing analysis from RNA-Seq[J]. Biorxiv，2014，（1）：257.[15] Jaganathan K，Kyriazopoulou Panagiotopoulou S，Mcrae J F，et al. Predicting Splicing from Primary Sequence with Deep Learning[J]. Cell，2019，（3）：176.[16] Cleveland，D. W. Antisense oligonucleotide therapy for neurodegenerative disease[J].The Journal of clinical investigation，2006， 116（8）：2290.。