我最早碰rMATS的时候其实被可变剪切这四个字吓住了。那时候我手里只有一堆RNA-seq的BAM文件导师丢过来一句跑一下差异可变剪切看看然后甩给我一篇rMATS的参考文献。当时的第一反应是这玩意儿是能跑的后来下载了工具、装完依赖、转完GTF、三次跑崩之后我才真正意识到——rMATS本身并不复杂真正会让新手卡住的全是那些藏在水面下的环境问题和输入文件规范。这篇博文就基于我完整的实操记录来写从零开始把rMATS可变剪切分析跑通再把输出结果逐字段拆开讲清楚最后附上我踩过的坑和排查思路。无论你是刚入门生信的研究生还是已经跑过比对、正准备接触可变剪切分析的同学这篇文章应该能帮你少走一大段弯路。1. rMATS到底在算什么东西可变剪切事件与统计模型的底层逻辑1.1 五种可变剪切事件类型先认识你的分析对象可变剪切翻译成人话就是同一个基因转录出前体mRNA之后外显子可以被不同方式连接最后产生多种成熟mRNA翻译出功能可能完全不同的蛋白质。rMATS要做的就是在两组样本之间找到那些剪切模式发生显著变化的基因事件。rMATS把事件划分为五个经典类型事件类型缩写通俗解释Skipped ExonSE一个外显子在某些样本中被保留在另一些样本中被跳过Alternative 5 Splice SiteA5SS5端剪接位点发生改变外显子边界移动Alternative 3 Splice SiteA3SS3端剪接位点发生改变外显子边界移动Mutually Exclusive ExonsMXE两个外显子互斥只能有一个保留在成熟mRNA中Retained IntronRI内含子在部分样本中未被切除被保留在成熟mRNA中SE是最经典也最常见的可变剪切事件也是绝大多数文章里第一个放出来的图。A5SS和A3SS看的是剪接位点的滑动MXE比较典型但数量通常较少RI虽然在有些场景下很有生物学意义但因为内含子保留往往和RNA未成熟降解有关分析时需要额外谨慎。rMATS做的事情本质上就是对这类事件做量化 统计检验。量化指标叫做inclusion level也叫PSIPercent Spliced In指的是一条包含某个外显子/剪接位点的转录本占该基因所有转录本的比例。两组样本的PSI差异就是rMATS输出的核心结果。1.2 rMATS的统计模型为什么值得信任rMATS用的统计模型是一个分层贝叶斯框架。简单说它考虑到每个样本内部reads计数的抽样变异也考虑到同组样本之间的生物学变异。它不会像早期某些工具那样直接把所有样本的reads堆在一起算一个比例而是把样本间的离散度也纳入了估计。这一步太关键了。比如你对照组有3个样本处理组有3个样本一个基因在组内有很大的样本间差异如果忽略这层变异很容易跑出一堆假阳性。rMATS输出的PValue就是基于这个模型得到的FDR则是对多重检验校正后的值。另外要记住一个点rMATS的JC和JCEC两种结果分别代表两种不同严格程度的计数方式。JCJunction Counts只统计剪接位点处的junction reads要求比对到跨越外显子边界的reads上。JCECJunction Counts and Exon Counts在JC基础上额外增加了落在外显子内部的reads。你的reads是150bp还是100bp对JC与JCEC的选择没有绝对的对错之分但通常建议先看JC因为junction reads是剪切事件最直接的证据。2. 装环境、转GTF、备BAM跑通rMATS前最容易翻车的三个环节2.1 环境安装直接上conda别自己折磨自己rMATS本身是用Python 2写的而且依赖Boost C库手动编译rMATS的时候经常会遇到boost版本对不上、gcc版本太新导致编译失败之类的问题。我一开始就是自己下载源码包编译折腾了将近一天最后老老实实换成了conda。推荐直接用bioconda安装conda create -n rmats python2.7 -y conda activate rmats conda install -c bioconda rmats -y装完验证一下rmats.py --help如果能正常打印出帮助信息说明安装就没问题了。conda会自动把Boost、samtools等依赖打包处理好省去了一大堆手工编译的麻烦。需要rMATS的turbo版本时也可以用conda安装rmats2samtools或rMATS-turbo。不过我早期先用标准版跑通了全部流程这里就以标准版为例来写。2.2 转GTFEnsembl注释文件的格式陷阱rMATS对GTF文件的要求比较严格。它需要每个转录本或外因子的属性里包含gene_id和transcript_id而且Ensembl的GTF里gene_id的格式可能是gene_id ENSG00000141510;同时还有gene_name、gene_biotype之类的附加字段。实际运行中我遇到过一个报错提示找不到transcript_id最后发现是rMATS解析GTF时对属性字段的顺序和格式有依赖。一个稳妥的解决办法是先用rmats2samtools工具集里的脚本对GTF做一次规范化转换rmats2samtools --gtf input.gtf --out output.gtf如果不想额外装这个工具也可以直接用GFFread把GTF转成更干净的形式gffread input.gtf -T -o output.gtf转换后的GTF建议做一个简单检查确认里面同时存在exon、CDS、transcript行并且gene_id和transcript_id字段完整再进入rMATS才不会出幺蛾子。2.3 BAM文件的准备比对参数与排序规范rMATS的输入本质上是一组BAM文件加GTF注释文件。它不做从头比对所以你的上游比对质量会直接影响rMATS的准确性。我用的上游比对工具是STAR两个关键参数必须注意。第一个是--outSAMstrandField intronMotifSTAR输出时会根据内含子基序补上XS标签这对于rMATS这类依赖剪接位点信息的工具很有帮助。第二个是--outSAMtype BAM SortedByCoordinate直接输出按坐标排序的BAM省得后面再手动排序。如果用的是HISAT2也要确保最终BAM经过samtools sort排序并建立索引samtools sort - 8 -o sample.sorted.bam sample.bam samtools index sample.sorted.bamrMATS对BAM文件的染色体命名和GTF的命名要求一致。我踩过一个大坑比对参考基因组用的是UCSC的染色体命名chr1、chr2GTF却用的是Ensembl的注释1、2结果rMATS跑出来的事件数直接少了80%因为大部分reads根本匹配不上注释。提示跑之前全局检查一次BAM头部和GTF的染色体命名格式务必要统一。这一步排查只要两分钟却能省下后面半天时间。3. 第一个完整运行配置文件、readLength与命令参数详解3.1 group文件怎么写b1.txt与b2.txt的格式规范rMATS通过两个文本文件来指定对照组和处理组的BAM路径。每行一个样本路径写绝对路径一个典型的内容长这样/path/to/control_rep1.sorted.bam /path/to/control_rep2.sorted.bam /path/to/control_rep3.sorted.bam/path/to/treat_rep1.sorted.bam /path/to/treat_rep2.sorted.bam /path/to/treat_rep3.sorted.bam文件里千万不要有多余的空行路径里也不要带空格。曾经因为样本目录名里带了一个空格rMATS报了各种诡异的错误最后查了半天才发现是路径解析的问题。另外建议b1.txt和b2.txt的样本顺序保持一致比如对照组三个重复在前处理组三个重复在后这样后续看IncLevel1和IncLevel2列时会更容易对应。3.2 核心参数逐个拆解readLength、gtf、od与tmp一个标准的运行命令长这样rmats.py --b1 b1.txt --b2 b2.txt \ --gtf gencode.v38.annotation.gtf \ --od output_dir --tmp tmp_dir \ --readLength 150 \ --nthread 10 \ --tstat 4每个参数的坑我都踩过一个个说--readLength指的是测序读长。这一步直接决定rMATS对reads的过滤策略。如果你用的是150bp双端测序就填150如果是100bp就填100。填错了最典型的现象是事件检出的数量异常多或者异常少因为rMATS会用readLength来判断一个reads是否可能跨越某个剪接位点。实在不确定的时候可以用samtools view看一下BAM里的一条记录samtools view sample.sorted.bam | head -5第10列就是read长度取多数reads的长度填进去即可。--gtf指定注释文件路径。--od是输出目录--tmp是临时文件目录两个都要预先创建好rMATS不会自动创建不存在的目录。--nthread是并行线程数建议不超过你机器CPU核数的70%。--tstat是最后一步统计检验时的线程数通常给4到8就够。如果机器内存不够大跑RI事件时会非常吃内存这时候可以把tstat调小一点。3.3 跑起来之后怎么看进度rMATS运行过程分为两个阶段第一阶段是对每个事件统计reads计数输出各个事件的原始计数文件第二阶段是计算inclusion level并做统计检验。如果运行过程中日志卡在某个地方很久很可能不是死机而是在处理比较大的染色体或者密集的剪切事件区域。可以用如下命令观察CPU和内存占用top -u $USER标准版rMATS默认会先把事件分类信息写到临时目录最后再汇总。如果一切顺利输出目录里会出现SE.MATS.JC.txt SE.MATS.JCEC.txt A5SS.MATS.JC.txt A5SS.MATS.JCEC.txt A3SS.MATS.JC.txt A3SS.MATS.JCEC.txt MXE.MATS.JC.txt MXE.MATS.JCEC.txt RI.MATS.JC.txt RI.MATS.JCEC.txt看到这十个文件说明rMATS已经完整跑通了。4. 输出文件到底怎么读JC与JCEC差异、每个字段逐列拆解4.1 JC和JCEC差的那一个E到底差在哪同一个事件类型会出现两个文件比如SE.MATS.JC.txt和SE.MATS.JCEC.txt。这两个文件的行数一样但计数方式不同。JCJunction Counts只统计跨越剪接位点的junction reads这类reads直接证明了外显子的连接关系特异性最高。JCECJunction Counts and Exon Counts则额外统计了完全落在组成型外显子和可变外显子内部的reads然后把这些reads也按一定规则分摊到inclusion或skipping的计数中。JCEC相当于把外显子本体表达量也拉进了估算会让PSI估计更稳定但也可能引入一些来自未剪切pre-mRNA的背景噪音。大多数论文用的是JC结果作为主结果JCEC作为稳健性验证。我个人的习惯是先看JC的显著事件列表再用JCEC的结果做一次交叉验证两边都显著的才进入后续分析。这比单纯用某一套结果要可靠得多。4.2 关键字段逐个解析从ID到IncLevelDifference打开SE.MATS.JC.txt表头大概有这些列列名含义IDrMATS给每个事件分配的唯一编号GeneID基因的Ensembl IDgeneSymbol基因符号chr染色体strand正负链exonStart_0base可变外显子起始位置0-base坐标exonEnd可变外显子结束位置upstreamEE上游组成型外显子的结束位置downstreamES下游组成型外显子的起始位置IJC_SAMPLE_1对照组样本中支持inclusion的junction countsSJC_SAMPLE_1对照组样本中支持skipping的junction countsIncFormLeninclusion isoform的有效长度SkipFormLenskipping isoform的有效长度PValue未校正的p值FDR多重检验校正后的p值IncLevel1对照组每个样本的inclusion levelPSIIncLevel2处理组每个样本的inclusion levelPSIIncLevelDifferenceIncLevel2的平均值减去IncLevel1的平均值IncLevel1和IncLevel2这两列值得仔细说。它们并不是一个值而是以逗号分隔的多个值每个值对应该组内的一个样本。样本顺序和你写在b1.txt、b2.txt里的顺序严格一致。如果某个样本在当前事件上没有足够reads支持那个位置的值为NaN。筛选时只要发现NaN比例过高就说明这个事件在部分样本中没有表达量即便是显著事件也要打个问号。IncLevelDifference是核心筛选字段。它代表处理组平均PSI减去对照组平均PSI正数表示处理组inclusion level升高也就是外显子更多地被保留负数则表示处理组倾向于跳过这个外显子。4.3 我先看一眼数据再谈筛选拿到文件之后我建议先用简单的shell命令快速扫一遍分布wc -l SE.MATS.JC.txt head -1 SE.MATS.JC.txt再用awk看一下显著的条目数awk -F\t NR1 $13 0.05 {print $0} SE.MATS.JC.txt | wc -l注意这里的$13是PValue列如果你的表头顺序不同列号会有变化。更稳妥的方式是用data.table或者pandas来读按列名筛选。5. 筛选显著事件与sashimiplot可视化从数据表到论文图5.1 我的筛选策略FDR 差异幅度 计数支持度很多教程只会告诉你FDR小于0.05就是显著。如果你真的直接把这个标准套上去大概率会得到几百个显著事件但其中有不少是极低表达量带来的噪音。我实际用的筛选标准是三个条件同时满足FDR 0.05。这是统计学显著性的基础门槛。abs(IncLevelDifference) 0.1。这一条的意义是FDR小只能说明差异不是随机波动并不代表差异足够大。0.1的PSI变化意味着至少10%的转录本发生了剪切模式的改变生物意义上通常比较可信。如果你的样本量很大FDR很容易小但IncLevelDifference只有0.02这种事件写进文章容易被审稿人质疑。IJC和SJC的计数支持度。在两组样本中至少有一组的平均计数不低于5。过滤掉那些几乎不表达的事件。用R或Python筛完以后我一般还会用IncLevel1和IncLevel2画出每个事件的PSI分布图肉眼确认一下组间差异不是单个样本拉动的结果。在只有两三个生物学重复的实验里这一步尤其重要。5.2 rmats2sashimiplot把显著事件变成sashimi图拿到了显著事件列表下一步是可视化。sashimi图可以展示剪切事件的isoform结构以及两组样本在剪切位点处的reads覆盖度也是论文里最常见的可变剪切展示方式。推荐用rmats2sashimiplot。安装方式同样推荐condaconda install -c bioconda rmats2sashimiplot -y画图前需要准备两样东西一是rMATS输出的SE.MATS.JC.txt文件二是一份BAM文件路径的配置文件。这个配置文件的格式长这样control_rep1:/path/to/control_rep1.sorted.bam control_rep2:/path/to/control_rep2.sorted.bam treat_rep1:/path/to/treat_rep1.sorted.bam treat_rep2:/path/to/treat_rep2.sorted.bam然后对指定事件画图。比如想要画SE事件中的一个显著事件命令如下rmats2sashimiplot --b1 bam1.txt --b2 bam2.txt \ -t SE -e 123 \ --l1 Control --l2 Treat \ --exonStart 1000 --exonEnd 2000 \ --annotation gencode.v38.annotation.gtf \ --output-dir sashimi_output \ --read-length 150这其中的-e 123指的是SE.MATS.JC.txt里你要画的第123行事件。--exonStart和--exonEnd是可变外显子的坐标范围用来限定画图的窗口宽度。--read-length同样要填150和rMATS运行时保持一致。画完以后sashimi图里上方是isoform模型下方是reads覆盖度曲线junction reads会用弧线连接起来弧线上标的数字就是该处的junction read计数。这两组之间的差异一眼就能看明白。有一个小习惯我在跑sashimiplot时会把--exonStart和--exonEnd适当往外扩一点让图里带上上下游组成型外显子的部分这样审稿人看图时能直接理解可变外显子的位置关系。5.3 画图之前顺手做的坐标检查rmats2sashimiplot对坐标要求很严格如果在画图时遇到Event not found或Error parsing coordinate之类的提示先回到SE.MATS.JC.txt里看那几列坐标。Sashimiplot要求的坐标是基于事件表格里exonStart_0base和exonEnd的如果是从其他文件里copy的坐标要确认到底是不是0-base否则容易产生偏移。另外画图用的BAM必须和rMATS运行时用的BAM完全一致最好在画图前重新samtools index一次避免因为索引文件被改动而导致画图中断。6. 跑批时碰到的问题排查记录几条真实报错的定位思路6.1 AN ERROR OCCURRED while reading GTFGTF解析失败这个报错我遇到过两次第一次是因为GTF里缺少transcript_id属性第二次是因为GTF文件用文本编辑器打开并重新保存后编码被改成了带BOM的UTF-8格式导致rMATS解析首行时把BOM字符当成字段内容。排查思路很简单head -3 your_annotation.gtf如果第一行开头有一个肉眼看不见的\ufeff用sed去掉即可sed -i 1s/^\xef\xbb\xbf// your_annotation.gtf同时确认GTF的第2列source、第3列feature结构完整。rMATS主要依赖feature为exon的行如果GTF里exon行本身就缺失后面的事件识别会完全跑空。6.2 事件数异常稀少染色体命名不一致我最初跑的GTF来自Ensembl当时的版本没有chr前缀而STAR比对用的索引是UCSC来源的BAM里全是chr1、chr2。rMATS正常启动了没有任何报错但输出的事件数量少得离谱——SE事件只有几十个正常应该有几万个。这种情况最让人头疼因为所有环节都是成功结束的只有结果数量和预期严重不符。排查方法是在比对后的BAM里随机提取一条记录的染色体名再在GTF里查一下同一个基因是否存在确认命名体系是否一致。samtools view sample.sorted.bam | head -1 | cut -f3 grep -m1 gene_id \ENSG your_annotation.gtf | cut -f1两个输出的前缀一对比问题立刻暴露。6.3 readLength参数填错导致的PSI值大面积异常有一批数据实际是双端150bp测序我不小心填成了100。跑完之后发现输出的事件数量比预期的多了很多而且很多IncLevelDifference的值算出来特别大接近0.9甚至1.0。排查时先怀疑是GTF的问题后来检查BAM里read length才反应过来把readLength改回150重新跑结果明显恢复正常。这里的原理是rMATS会把reads比对到外显子连接处或内部readLength影响的是某个reads是否能被判定为有效跨越事件连接处。填短的readLength会让rMATS认为更多reads满足跨越条件从而增加有效计数填长则相反。务必用samtools确认真实读长后再填别偷懒。6.4 内存不足导致RI事件计算中断标准版rMATS对RI事件的处理会消耗较多内存尤其是在注释文件较大、染色体覆盖范围广的情况下。表现为日志在处理某个染色体时报MemoryError或进程直接被OOM Killer杀掉。排查思路不是去调整rMATS内部参数而是从运行策略上绕开。你可以把RI事件单独跑因为RI在植物和动物研究里的事件数量通常相对有限单独跑可以避免峰值内存叠加。也可以考虑使用rMATS-turbo版本它的内存占用优化会好很多适合大规模样本。6.5 结果文件行数不一致同一事件在两套输出里的排序问题有时候你会看到SE.MATS.JC.txt和SE.MATS.JCEC.txt的行数不一样少一行或多一行。这通常和运行时某些事件在某一种计数方式下得不到有效reads有关。rMATS的JC和JCEC虽然针对相同事件但两套结果文件的ID是一一对应的行数不一致时别急着去合并两个文件先按ID做匹配再merge。我的做法是用Python的pandas以ID列为key做left join。如果一定要用Excel操作记得把ID列设置成文本格式否则Excel可能自动转成科学计数法导致ID失真。7. 把rMATS结果推进下游分析GO/KEGG富集与课题故事线7.1 从基因列表到富集分析的基本套路筛选出显著可变剪切事件后下一步逃不开的问题是这些剪切变化影响了哪些生物学通路。此时你需要从显著事件表里提取GeneID或geneSymbol去重后作为基因列表做GO和KEGG富集分析。推荐用clusterProfilerR包或者在线工具DAVID。clusterProfiler处理Ensembl ID时会自动帮你做ID转换但注意它需要联网获取注释数据。如果服务器没有外网条件建议提前下载好OrgDb包或KEGG数据文件在本地完成富集。一个常见的坑是可变剪切显著基因和差异表达基因DEG列表往往不重合这就对了。可变剪切变化代表的是转录本层面的调控与表达量层面的调控是两个维度。如果你发现两组列表高度重合反而要检查一下是不是rMATS的计数受表达量影响太大产生了表达量驱动的假阳性。7.2 挑选候选事件时的生物学判断富集结果出来以后我一般会从显著事件里挑出富集通路中那些既在通路里有明确功能、又在剪切模式上有较大差异的事件作为后续湿实验验证的候选。挑选时我还会重新回到IGV里人工看一遍。IGV加载BAM文件和GTF注释直接在基因组浏览器里看reads覆盖和junction情况能确认这个事件在真实数据中是否清晰可见。sashimi图是成品展示IGV则是自查手段。很多时候sashimi图给出的junction read计数是准的但IGV里一拉发现该区域比对质量很差或者有重复序列干扰这样的候选事件我会直接放弃。7.3 从单个事件到全局模式的验证如果处理后得到的显著剪切事件数量比较多可以考虑做一个全局剪切模式的汇总比如统计所有显著事件中inclusion上调与下调的比例、在不同染色体上的分布、事件类型偏好性。这些不仅仅是论文里的描述性统计也能帮你判断rMATS结果是否真的具有生物学逻辑。比如在某一组处理条件下如果SE事件整体偏向于inclusion水平下降意味着处理条件下更多外显子被跳过可能暗示该处理抑制了剪切复合体的某些核心组分的功能。这种模式如果能和已知的生物学背景对应上故事线就顺了。8. 关于rMATS分析流程最后再分享几个实际体感整个rMATS流程跑下来我最深的感受是这个工具对新手最大的门槛不在算法本身而在Linux环境、上游比对、注释格式这些看起来不起眼的细节上。你花两个小时把环境装好花半小时确认GTF和BAM格式统一后面跑起来就会顺很多。另外rMATS跑完只是分析的开始真正花时间的往往是从显著事件表到最终生物学结论的过程。我建议你从一开始就给每个样本做好命名规范把b1.txt、b2.txt和输出目录按项目维度管理否则等你要复现结果、补充分析的时候光是对照文件路径就能花掉一下午。最后再多说一句如果是刚接触Linux的同学别在rMATS卡住的时候盲目加参数或者反复重跑。先学会看日志、学会用samtools view和head去检查文件内容很多问题都能自己定位。rMATS的报错信息虽然不友好但只要你能定位到是哪一类输入文件出了问题问题就解决了一大半。