干货指南｜一文教会你轻松掌握可变剪切结果的解读与运用

可变剪切（Alternative Splicing, AS）是真核生物基因调控的重要机制：DNA转录形成前体mRNA（pre-mRNA）再加工产生成熟mRNA的过程中，会有选择的对exon/intron及其他元件进行剪切和重组，因此一个基因可以产生多种类型的转录本，进而使得一个基因可以在不同时间或不同环境中翻译出不同的蛋白质，增强生命系统的复杂性和适应性。

可变剪切（AS）事件主要分为5种类型：

1. 外显子跳跃 Skipping exons (SE)

2. 5’端可变剪切 Alternative 5’ (donor) splice sites (A5SS)

3. 3’端可变剪切 Alternative 3’ (acceptor) splice sites (A3SS)

4. 内含子保留 Retained introns (RI)

5. 互斥外显子 Mutually exclusive exon usage (MXE)

还有两种不太常见的类型为：第一个外显子剪切 Alternative first exon (AF)和最后一个外显子剪切Alternative last exon (AL)。

目前已经开发出了多种计算方法来识别和量化RNA-seq数据中的可变剪切，需要知道的是，可变剪切一般通过计算支持剪切的reads数量，实现对可变剪切程度或比例的量化，进而用于单个样本或者组间样本的比较。简单理解便是以下公式：

PSI = splice_in / (splice_in + splice_out)

在转录组分析中，二代转录组会使用rMATS进行分析[1]，反映剪切水平的指标为IncLevel，其计算逻辑与上述公式基本一致，根据剪接事件所产生同工型的有效长度以及支持可变剪切事件的reads数量进行计算，结果反映样本的最终mRNA转录本中平均包含各自外显子的频率的信息，也就是剪切事件exon inclusion isoform在两个isoform总count数的比值。三代转录组中会使用SUPPA2进行分析[2]，反映剪切水平的指标为PSI（percent spliced in），计算逻辑大致与rMATS相似，根据发生可变剪切与正常剪切产生的isoform的表达量TPM值计算。下面小编将简要解读可变剪切的分析结果并提供一些可变剪切挖掘及验证思路以供参考。

rMATS是一款适用于二代转录组数据的差异可变剪切分析软件，可以对可变剪切事件进行分类鉴定以及差异分析。rMATS可以识别的可变剪切类型为SE、A5SS、A3SS、RI和MXE。

rMATS可以识别的可变剪切类型

rMATS定量可变剪切基于两种定量方式：Junction count only和Reads on target and junction counts，在结果中分别为包含*.JC和*.JCEC的表格。JC指Junction Count，表示跨越剪切位点的reads数目；JCEC是Junction Count和Exon Count的合并，不仅考虑跨越剪接位点的reads（Junction Count），还同时考虑未跨越剪接位点的reads（Exon Count）。如果只是比较两组样品间可变剪切的差异一般使用JC的结果即可。另外，结果中JC和JCEC的表格均可直接作为rmats2sashimiplot可变剪切可视化的输入文件。

在单样本的可变剪切事件鉴定结果中，第一列是AS事件的编号，从0开始；第二列（GeneID）是基因的编号；第三列（geneSymbol）是基因名称，如果GFF/GTF文件不包含geneSymbol信息，此列为NA；第四列（chr）是基因所在染色体的编号；第五列（strand）是基因所在链的方向。随后几列是剪接位点的位置信息。具体每种事件的剪接位点位置信息类型可以参考：

https://github.com/Xinglab/rmats-turbo/blob/v4.3.0/README.md#output

在组间比较的差异可变剪切结果中，前十一列和单样本的可变剪切事件结果相似，从第十二列起有所差别，第十二列是ID，同第一列；第十三列（IJC_SAMPLE_1）是分组1的inclusion junction（IJC）的count数；第十四列（SJC_SAMPLE_1）是分组1的skipping junction（SJC）的count数；第十五列（IJC_SAMPLE_2）是分组2的IJC的count数；第十六列（SJC_SAMPLE_2）是分组2的SJC的count数，同组样本的结果以逗号分隔。第十七列（lncFormLen）是AS事件Exon Inclusion Isoform的有效长度；第十八列（SkipFormLen）是AS事件Exon Skipping Isoform的有效长度；第十九列（PValue）是两组样品可变剪切的显著差异指标；第二十列（FDR）是对PValue的校正值；第二十一列（lncLevel1）是分组1中AS事件Exon Inclusion Isoform在两个Isoform总count数的比值；第二十二列（IncLevel2）是组2中可变剪切事件Exon Inclusion Isoform在两个Isoform总count数的比值；第二十三列（IncLevelDifference）是lncLevel1均值与IncLevel2均值的差值，其中IncLevel1对应处理组，IncLevel2对应对照组。更为具体的可以参考：

https://www.jianshu.com/p/99a626391b04

SUPPA2是一款适用于三代转录组数据的可变剪切分类鉴定和差异分析的软件，同样是基于AS事件定量。SUPPA2可以识别的可变剪切类型为SE、A5、A3、RI、MX、AF和AL。

SUPPA2可以识别的可变剪切类型

软件输出结果中包含单样品可变剪切鉴定和差异可变剪切鉴定结果。单样本可变剪切鉴定结果为*.psi的表格，第一列（gene_id）是发生可变剪切事件的基因ID，第二列（event_id）是发生的可变剪切的事件名称，组成为：基因ID;可变剪切类型;基因所在染色体名称;剪切位点;正负链，第三列（alternative_transcripts）是该基因通过可变剪切形成的转录本的ID；第四列（total_transcripts）是该基因在剪切位点区域通过正常剪切和可变剪切形成的全部转录本的ID，第五列（PSI）是发生可变剪切的程度，这里的计算方法是PSI=发生该可变剪切转录本TPM值/(发生该可变剪切转录本TPM值+未发生该可变剪切转录本TPM值)，如果PSI为NAN或者0，则认为该基因在该样本中未发生可变剪切。另外，包含*psi.add_info的表格为对发生可变剪切事件的基因添加功能注释的结果，包含*.pie的图为不同类型的可变剪切事件在该样本中的数量分布比例情况。

在组间比较的差异可变剪切结果中，文件名为*gene.dpsi.add_info的表格为各AS事件在组间差异程度的分析，前五列与单样本结果相似，第一列（gene_id）是发生可变剪切事件的基因ID，第二列（event_id）是发生的可变剪切事件的名称，组成为：基因ID;可变剪切类型;基因所在染色体名称;剪切位点;正负链，第三列（alternative_transcripts）是该基因通过可变剪切形成的转录本的ID；第四列（total_transcripts）是该基因在剪切位点区域通过正常剪切和可变剪切形成的全部转录本的ID，*PSI列为对应样品发生可变剪切的程度，dpsi列是 psi difference(ΔPSI)，即两组PSI值取平均值后相减，dpsi不为NAN或者0表明两组间可变剪切存在差异；diffAS-pvalue列是两组样品可变剪切的显著差异指标，一般pvalue＜0.05，认为是显著差异。gene_name 列为geneSymbol，仅在参考基因组GFF/GTF注释文件中包含相关信息时添加，随后几列是基因的相关功能注释。

真核生物基因转录加工过程往往受到不同环境或者条件的影响，通过可变剪切形成不同的转录本，进而编码不同的蛋白质，因此从可变剪切产生的转录本本身出发可以深入挖掘剪切转录本的功能或者结构，例如Corre等[3]通过对李斯特菌感染肠上皮细胞进行ONT全长转录组测序，发现细菌感染产生毒素可以诱导关键调控基因CIRBP和CLK1受到可变剪切的影响，导致其产生NMD靶标等包含不同外显子结构的转录本，编码不同的毒素结合蛋白以响应细菌感染刺激；又如Cao等[4]利用ONT全长转录组测序技术鉴定幼年和成年大鼠心脏细胞不同发育时期全长Tpm1转录本，发现RBFOX2蛋白可以特异性调控Tpm1末端外显子的可变剪切，影响Tpm1的多聚腺苷酸化，从而产生具有不同3’末端的转录本，影响细胞生长和发育。

可变剪切结构分析^[3]

其次，可变剪切差异也是一个切入点，尽管基因普遍发生可变剪切，但是在不同处理或者不同阶段剪切的程度或者比例可能存在差异，由此可以探究剪切存在差异的基因以及产生的转录本富集参与哪些通路或者功能，例如Wu等[5]对不同处理下的人类主动脉平滑肌细胞进行ONT 全长转录组测序，鉴定分析发现差异可变剪切事件大多参与RNA代谢相关通路，并且组间显示出特定的剪切模式；又如Zhou等^[6]利用Direct RNA测序技术，鉴定分析中华鳖不同性别之间的可变剪切差异，发现A3、A5和AF等可变剪切事件数量在由雌性到假雄性的性逆转过程中显著减少，而在由雄性到假雌性的性逆转过程中显著增加。

差异可变剪切分析^[5]

此外，通过将可变剪切与基因表达或者修饰等多层面数据结合[7]，可以揭示更为复杂的生物学调控机制。例如与表达关联，通过比较基因水平无显著差异而转录本水平存在差异的基因，来探究可变剪切是否存在差异，同时还可以分析差异表达基因与差异可变剪切基因之间的交集，探讨这些基因可能参与的信号通路或发挥的功能，从而更好地理解转录本或基因的使用情况；此外，对于那些未发生差异可变剪切但存在差异表达的基因，可以从转录后调控的角度进行研究，如关注poly(A)尾巴长度或RNA修饰的变化，以进一步解释可变剪切事件的发生与变化。

可变剪切与表达及修饰关联^[7]

那如何对可变剪切进行后续的验证呢？可变剪切最常用且最简单的验证方法就是反转录PCR，其原理为在剪切事件发生区域设计特异引物以扩增特定片段，通过电泳分离产物并观察条带大小差异，从而验证是否存在不同的剪切形式。基本步骤为根据剪切位点信息确定发生的事件及借助IGV浏览器确定对应区域的序列，在序列区域附近设计特定引物进行PCR扩增，确定片段长度。当然也可以结合其他实验手段验证可变剪切的功能或者作用机制。

可变剪切验证^[8]

参考文献：

[1] Shen S, Park J W, Lu Z, et al. rMATS: robust and flexible detection of differential alternative splicing from replicate RNA-Seq data[J]. Proceedings of the National Academy of Sciences, 2014, 111(51): E5593-E5601.

[2] Trincado J L, Entizne J C, Hysenaj G, et al. SUPPA2: fast, accurate, and uncertainty-aware differential splicing analysis across multiple conditions[J]. Genome Biology, 2018, 19: 1-11.

[3] Corre M, Boehm V, Besic V, et al. Alternative splicing induced by bacterial pore-forming toxins sharpens CIRBP-mediated cell response to Listeria infection[J]. Nucleic Acids Research, 2023, 51(22): 12459-12475.

[4] Cao J, Routh A L, Kuyumcu‐Martinez M N. Nanopore sequencing reveals full‐length Tropomyosin 1 isoforms and their regulation by RNA‐binding proteins during rat heart development[J]. Journal of Cellular and Molecular Medicine, 2021, 25(17): 8352-8362.

[5] Wu H, Lu Y, Duan Z, et al. Nanopore long-read RNA sequencing reveals functional alternative splicing variants in human vascular smooth muscle cells[J]. Communications Biology, 2023, 6(1): 1104.

[6] Zhou T, Chen G, Chen M, et al. Direct full-length RNA sequencing reveals an important role of epigenetics during sexual reversal in Chinese soft-shelled turtle[J]. Frontiers in Cell and Developmental Biology, 2022, 10: 876045.

[7] Ma Q, Gui Y, Ma X, et al. N6-methyladenosine writer METTL16-mediated alternative splicing and translation control are essential for murine spermatogenesis[J]. Genome Biology, 2024, 25(1): 193.

[8] Zhang H, Shen X, Sun S, et al. Integrated transcriptome and proteome analysis provides new insights into camptothecin biosynthesis and regulation in Camptotheca acuminata[J]. Physiologia Plantarum, 2023: e13916.