我用rMATS做可变剪切分析已经有几年了从刚开始对着报错信息手忙脚乱到现在能一套流程稳定跑完十几个样本中间踩过的坑确实不少。市面上关于rMATS的教程很多但大多是告诉你命令怎么敲很少解释为什么要这么敲更不用说那些藏在细节里的坑了。这篇文章我不打算重复官方文档而是把自己从环境配置、文件准备、参数选择到结果解读的完整经验整理出来尽可能说清楚每一步的取舍和原理让第一次接触rMATS的菜鸟也能少走弯路。Linux环境下做可变剪切分析rMATS是绕不开的经典工具。它专门用于RNA-seq数据中差异可变剪切事件的检测支持SE、A5SS、A3SS、MXE、RI五类事件结果是带统计检验的表格方便后续筛选和可视化。适合的场景很明确你已经拿到比对好的BAM文件或者还在上游比对的阶段想比较两组样本比如疾病vs对照、处理vs未处理之间剪切模式的变化。这篇文章从零开始讲透完整流程应该可以帮你在实际项目中直接套用。1. 方案选型为什么是rMATS以及STAR比对为什么更稳1.1 rMATS相比其他可变剪切工具的核心优势可变剪切分析领域并不缺工具SUPPA2、LeafCutter、MISO、MAJIQ都有自己的用户群。我最初也曾纠结选哪个后来在实际项目里把rMATS当成主力主要是它在统计模型和结果可解释性上的平衡做得最好。rMATS用的是一个基于计数矩阵的模型对每个剪切事件把reads映射到包含或排除该外显子的类别上再通过似然比检验判断两组之间有没有显著差异。它输出的IncLevelinclusion level即外显子包含水平非常直观类似PSI值数值范围0到1反映的是某个外显子或剪切位点被使用的相对频率。下游展示时可以直接画sashimi图审稿人看着也舒服。还有一个很实际的优点rMATS对注释文件的要求不算苛刻常规的GTF就能跑。LeafCutter走的是无注释聚类路线能发现新事件但结果解释起来复杂适合做发现型研究SUPPA2计算速度快但依赖转录本表达量估算对定量环节敏感MISO模型经典但速度和扩展性都差一些。rMATS正好卡在“有注释、有统计、结果直观”这个位置日常分析够用发文章也站得住。1.2 上游比对的选择STAR为佳BAM必须按坐标排序rMATS的输入是BAM文件理论上任何比对软件的输出只要符合格式都能用但我还是建议用STAR。原因不是rMATS对STAR有特殊接口而是STAR的剪接比对能力更强识别junction reads的灵敏度和准确性更好这在剪切事件计数环节直接关系结果质量。我以前为了图省事用过HISAT2比对同样一份数据跑下来标准差异事件的数目确实有差距STAR检出的事件更多、在基因组浏览器上人工核对时假阳性也少。用STAR比对时要注意两个点。第一输出BAM必须按坐标排序也就是Sort后生成的Aligned.sortedByCoord.out.bamrMATS要靠坐标索引快速抓取特定区域的reads未排序的BAM会直接报错或者崩溃。第二建议保留包含junction的reads不要用--outSAMstrandField intronMotif这种参数强行加链信息除非你用的是链特异性文库且有特殊需求。rMATS本身能处理链信息它会参考GTF里的strand字段来分配reads不需要你在比对时额外折腾。1.3 单个样本比较与组间比较的取舍rMATS最舒服的使用场景是有生物学重复的分组比较。只有单个样本对单个样本rMATS也能跑但由于没有组内方差估计FDR和PValue很多情况下会是NaN结果基本没法用。如果你实验设计里只有两三个重复也尽量别用单样本模式去硬撑宁可多跑一个样本也不要让后续统计分析陷入尴尬。如果实在只有两个样本比如探索性分析可以用--s1和--s2参数直接指定两个BAM文件路径省去写文本列表的步骤。但要记住这种比较只能做趋势判断写文章时不能宣称有显著差异。2. 环境搭建别让依赖问题消磨掉你的热情2.1 用conda创建独立环境避免版本冲突安装rMATS最省心的方法是conda我会为每个生信项目单独建环境避免软件之间的依赖互相打架。创建命令如下conda create -n rmats python3.7 -y conda activate rmats conda install -c bioconda rmats4.1.2 -y这里特意指定了python3.7因为rMATS的pysam和numpy版本兼容性比较保守在太新的Python版本下有时会出一些莫名其妙的报错。4.1.2是我测试过比较稳定的版本官方后来发布的4.1.2之后更新不多但大版本内部小版本之间的差异也需要注意。装好后先验证一下rmats.py --help如果能正常打印参数说明说明基础环境没问题。我遇到过一个奇怪的情况是conda自动装了最新版rMATS调用时提示No module named numpy这时候别急着重装先看下环境里是否真的缺少numpy如果缺就conda install numpy1.21补上。2.2 常用配套工具清单除了rMATS本身还需要准备samtools、STAR、gffread这些常用工具。samtools用于bam索引和格式检查STAR用于上游比对gffread则是在注释文件格式有问题时的救急工具。建议一并装好conda install -c bioconda samtools star gffread -y2.3 Docker或Singularity作为备选方案有些集群环境不允许随便用conda创建软件环境或者管理员不让装新包这种情况下可以用Docker镜像。rMATS官方提供了Docker镜像用法也很简单docker pull comor/rmats:4.1.2 docker run -v /your/data:/data comor/rmats:4.1.2 rmats.py --b1 /data/b1.txt ...但这套流程有个麻烦需要把宿主机的目录映射到容器里路径容易混乱。我在实际项目中更推荐用conda只有在没有root权限的共享集群上才考虑容器方案。如果集群支持Singularity把它转换一下也能用singularity build rmats.sif docker://comor/rmats:4.1.23. 输入文件制备GTF和BAM里的坑几乎每个人都踩过3.1 GTF注释文件格式检查比想象中重要rMATS分析过程中注释文件决定了事件注释的准确度。不要直接拿一个下载完的GTF文件就开跑先检查文件是否包含必要的gene_id和transcript_id字段。Ensembl和GENCODE的标准GTF通常没问题但如果你从UCSC Table Browser下载或者经过其他工具转换可能丢失transcript_id这时rMATS会报错或者静默地忽略该转录本导致事件检测数量偏少。可以用下面的命令快速检查less GCF_000001405.39_genomic.gtf | head -n 5 # 看看第9列是否同时有gene_id和transcript_id如果发现GTF缺少gene_id和transcript_id可以用gffread从GFF3格式转换gffread genome.gff3 -T -o output.gtf但转换后的GTF也可能存在转录本ID混乱的问题转换完再用R或awk随机抽几行检查一下字段。染色体编号一致性也要注意。如果BAM文件里的染色体是1、2、3这种不带chr前缀的格式GTF里也必须是同样风格二者不一致的话rMATS找read时完全对不上结果里事件数量会少得离谱。检查方式samtools view -H sample.bam | grep ^SQ | head -n 5 grep -v ^# annotation.gtf | cut -f1 | sort -u | head -n 5两边对比一下发现差异就用sed或者samtools reheader统一格式别硬跑。3.2 BAM文件排序、索引和去重这三个问题BAM文件除了按坐标排序外还必须要建索引samtools index sample.bamrMATS读取BAM时需要依托索引来快速定位区域没有索引会直接报错。检查是否已有索引只需看bam同目录下是否有同名.bai文件。去重问题是我踩过最深的坑之一。RNA-seq比对后的BAM一般不建议做PCR去重因为同一转录本的多个reads本来就可能比对到相同位置去重之后junction reads的数量会被严重压缩可变剪切计数就失真了。我第一次用rMATS时上游流程里顺手加了MarkDuplicates这一步最后检出的差异事件少得可怜检查IncLevel数据才发现几乎所有事件的计数都大幅缩水。RNA-seq分析中除非你确信文库存在严重PCR偏好否则保留所有reads。被rMATS忽略的多比对reads同样值得注意。STAR比对时如果允许一个reads比对到多个位置rMATS默认会跳过比对质量低或者多匹配的reads这个策略相对保守能有效减少假阳性但也会损失一部分灵敏度。想保留多比对reads带来的信息也可以但需要你清楚自己数据的具体情况一般标准流程不建议调整太多。3.3 样本分组文件的写法rMATS支持用文本文件指定每个组里的BAM文件路径一行一个。比如b1.txt内容/path/to/control_1.bam /path/to/control_2.bam /path/to/control_3.bamb2.txt写处理组的BAM路径/path/to/treatment_1.bam /path/to/treatment_2.bam /path/to/treatment_3.bam路径用绝对路径最保险相对路径在rMATS内部工作目录变化时容易出问题。文本文件的换行符也要注意Windows下编辑过的文件携带着\r字符会导致rMATS找不到文件用sed -i s/\r$// b1.txt清理一下就好。4. 运行rMATS参数说明和一次完整的实操4.1 推荐的最简运行命令环境准备好、输入文件检查无误后下面的命令可以一次性完成从计数、统计到输出差异事件的全过程rmats.py \ --b1 b1.txt \ --b2 b2.txt \ --gtf annotation.gtf \ -t paired \ --readLength 150 \ --nthread 8 \ --od output \ --tmp tmp解释一下关键参数-t paired测序数据是双端还是单端必须和实际文库一致。--readLength 150测序读长一般填你数据中大多数reads的长度。如果一个文库里有150bp和151bp混着的情况可以加上--variable-read-length参数让rMATS兼容不同长度。--nthread 8线程数不设默认也能跑但时间会慢很多。--tmp中间文件目录跑完以后通常可以删掉但有些调试场景需要保留。输出目录会自动创建里面会生成A3SS、A5SS、MXE、RI、SE五个目录每个目录下都有.MATS.JC.txt和.MATS.JCEC.txt两个结果文件。4.2 关于--novelSS和--cstat的补充说明rMATS有一个--novelSS参数开启后能检测新剪接位点也就是未被注释的外显子边界。但基于注释文件的经典分析模式下开启它可能会引入大量难验证的事件增加后续人工筛选的负担。我个人的经验是标准分析先不开等看过基本结果后再考虑是否探索novel splice site。--cstat参数是显著性阈值默认0.0001这个值控制的是pvalue cutoff的候选事件过滤不是最终FDR阈值。实际筛选时用FDR小于0.05、|IncLevelDifference|大于0.1做主要标准就够了cstat保持默认即可。4.3 运行时间与失败后的检查顺序rMATS跑起来很慢尤其是全基因组数据。我测过一个有6个处理组样本和6个对照组样本、每个样本30M reads的人类全转录组数据16线程跑完大约花了五六个小时。跑之前建议用nohup放到后台避免终端断开导致任务中断nohup bash run_rmats.sh rmats.log 21 运行过程中如果失败首选看log文件末尾的报错信息。按我遇到的情况高频问题依次是GTF字段缺失、BAM坐标系和GTF不一致、样本文件路径出错、内存不够。前两个问题前面已经说过怎么处理路径出错一般检查b1/b2文本文件有没有额外空格或者换行符内存不够就减少线程数或加大机器内存配额。5. 结果文件解读别只盯着PValueIncLevel才是核心5.1 JC与JCEC两种结果文件的区别每个事件类型下都会输出两个字文件命名里带.JC.的是仅用跨越剪切位点的junction reads来计算结果更保守可靠带.JCEC.的会在JC基础上再加入剪切位点附近的exonic reads灵敏度更高但也会受到未成熟信使RNA的影响假阳性概率更高。日常分析我会先看JC结果用它做主筛选再用JCEC的结果交叉验证。如果某个事件在JC里不显著但在JCEC里很显著而且生物学上听起来很有意思那就值得用IGV或sashimi图人工看一眼不要轻易丢。5.2 结果表格里每一列在说什么拿SE事件的结果举例打开SE.MATS.JC.txt后常见的列有IDGeneIDgeneSymbolchrstrandexonStart_0baseexonEnd upstreamESupstreamEEdownstreamESdownstreamEE IncLevel1IncLevel2IncLevelDifferencePValueFDR这里的exonStart_0base和exonEnd是被检测的盒式外显子cassette exon的坐标注意start是0-based的传到IGV或UCSC浏览器里通常能直接用但如果你是拿去做自定义可视化记得做0-based和1-based之间的换算。IncLevel1对应的是b1.txt这一组样本的平均inclusion levelIncLevel2对应b2.txt这一组。IncLevelDifference等于IncLevel1减去IncLevel2正值说明第一组的这个外显子更容易被包含负值说明第二组更容易包含。方向不要搞反我见过不少人把正负号解释反了最后结论全部反了。5.3 差异事件筛选阈值和后续验证常用的筛选条件是FDR 0.05 |IncLevelDifference| 0.1如果组间IncLevel差异绝对值小于0.1即使p值很小生物学意义也不大。剪切调控导致的外显子使用改变如果只有5%的水平验证实验很难重复出来审稿人也容易质疑。筛选时可以用R或者Excel但建议直接在Linux里用awk快速出一份过滤后的列表awk -F\t $20 0.05 $22 0.1 SE.MATS.JC.txt filtered_SE.tsv列数需要根据实际文件头调整第20列一般是FDR第22列是IncLevelDifference跑之前先head -n 1确认一下列号再动手。筛选完并不是终点。我会从过滤后列表里挑top 5到10个事件去IGV里加载BAM文件和GTF人工核对reads覆盖图形是不是真的符合事件类型定义。这个步骤花不了半小时但能帮你大幅提升结论的可信度。6. 可视化实操用rmats2sashimiplot快速出图6.1 软件安装和输入准备拿到结果之后最常用的可视化工具是rmats2sashimiplot。它读取rMATS输出的事件坐标结合BAM文件和GTF画出经典sashimi图。安装方式conda activate rmats pip install rmats2sashimiplot如果安装时提示缺少matplotlib、pysam之类的依赖直接用conda安装即可conda install -c conda-forge matplotlib pysam6.2 从SE结果中提取事件信息并绘图rmats2sashimiplot的输入可以参照官方给的示例需要一个简单的文本文件每行是事件坐标列分别为事件类型、chr、strand、外显子起止和内含子起止等信息。一个比较取巧的办法是直接复用rMATS结果里的坐标按它的示例格式整理。比如对于SE事件可以从SE.MATS.JC.txt中提取chr、strand、exonStart_0base、exonEnd、upstreamEE、downstreamES、PValue、FDR这些列。三组重要的位置分别是上游外显子的3末端upstreamEE、目标外显子的起止、下游外显子的5起始downstreamES。实际绘图命令示例rmats2sashimiplot \ --b1 control.bam \ --b2 treatment.bam \ -t SE \ -e event.txt \ --l1 Control \ --l2 Treatment \ --exon_s 1 \ --intron_s 5 \ -o sashimi_output--exon_s和--intron_s控制图形中外显子和内含子的缩放系数具体值按实际数据调整画出来不满意就调大调小不用太纠结。出图之后要做一个关键的人工确认看看目标外显子的reads分布是不是真的存在包含和排除两种模式。如果sashimi图里本该有junction支持的reads没有那这个事件可能是软件误判建议在结果里剔除。6.3 下游富集分析的简要思路拿到显著差异剪接事件后下一步往往是看看这些事件涉及的基因在什么通路里富集。由于每个基因可能对应多个事件我一般先把事件去重到基因层面再用clusterProfiler做GO和KEGG富集。不同事件类型SE、A5SS等如果一个基因同时出现去重时优先保留FDR最小的那个事件。富集分析我不展开细讲但有个提醒很多可变剪切相关基因并不会在普通差异表达分析里显著变化把它们单独拎出来做富集往往能发现一些不一样的通路。我做过一个肿瘤样本的数据差异表达基因富集不到肿瘤相关通路反倒是差异剪接基因在细胞黏附和免疫应答通路上显著富集后续验证实验也证实了相关剪接因子确实异常表达。这就是可变剪切分析独特的价值。7. 常见问题与排查技巧实录7.1 事件检出数量过少时先检查什么如果跑完结果里A3SS、A5SS这些文件里总共只有几百个事件而同类数据在文献里能检到几千个最常见的原因是GTF和BAM的染色体命名不一致。这个我前面提过但真的太常见了必须再强调一次。其次是STAR比对时用了过于严格的过滤条件导致junction reads数量太少比如用了--outFilterMismatchNmax 0这种参数比对到的reads极少rMATS自然什么也检不出。调出Log.final.out看STAR比对率如果unique mapping rate低于70%就要检查数据质量或者调整比对参数。比对率正常但事件少那就是注释或者组学数据的问题偏多。7.2 结果中PValue和FDR全是NaN这个现象多见于没有重复或者重复数只有两个且组内波动极端的情况。rMATS对每组的重复样本会估计组内方差如果某事件在组内完全无变化或变化为0方差估计就会变成0出现NaN。提高重复数量是终极解法但在数据量固定的前提下可以先只保留那些非NaN的事件再用IncLevelDifference绝对值做趋势筛选配合可视化验证。7.3 同一事件在JC和JCEC里的方向不一致理论上不该发生但如果你遇到优先相信JC结果。JCEC里额外计入的那部分reads有较大可能来自未剪切的pre-mRNA方向不一致说明某个组里pre-mRNA污染比较明显这种情况JC的方向更可信。7.4 跑了一半内存不足被killedrMATS比较吃内存尤其在同时加载大量bam索引信息的时候。降低--nthread并限制机器上的并发任务数量通常能解决。还可以尝试增加--tmp所在分区的磁盘空间rMATS中间会写不少临时文件磁盘满了也会被杀。用watch df -h看下磁盘占用提前清理别的东西。7.5 使用小规模测试数据快速验证流程正式跑全基因组数据之前我强烈建议先拿小规模数据把流程走通。可以只选一条染色体比如chr21或者chr1的一部分先比对、再跑rMATS整个过程几分钟出结果确认输出目录里能生成正常的事件文件再去跑全量数据。这样可以帮你把环境问题、格式问题和路径问题在最早期暴露出来省下后面数小时的等待。8. 我在实际项目中的几点体会做可变剪切真不是简简单单把命令跑通就行。整个流程下来最花时间的往往不是rMATS本身的运行而是前期的数据质控、格式检查以及后期的事件验证。尤其是BAM文件的质量直接决定了后面所有结果的可靠性。如果你发现自己的结果和预期相差很远我会建议先回头检查比对环节而不是急着调整rMATS参数。另外就是样本设计问题。rMATS虽然能处理两两比较但生物学重复的数量和质量是最硬的条件。条件允许的情况下每组至少三个重复尽量让组内一致性高一点否则统计这个环节会非常尴尬。最后再分享一个我习惯用的工作流跑完rMATS后除了过滤显著事件我还会把IncLevelDifference排名靠前但不显著的事件也留一个备份有时候这些边缘事件反而是某些样本特异的、有意思的变化。后续实验验证时这类事件偶尔会给你意外的惊喜。但写文章的时候还是老老实实以显著事件为准边缘信号只能当线索做探索不能当真结论写。这些经验都是拿真金白银的时间和失败的运行换来的。希望这份梳理能让你少踩一些我当年踩过的坑让你把更多精力放在生物学问题的解读上。