企业官网建设流程全解析

我用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映射到包含或排除该外显子的类别上，再通过似然比检验判断两组之间有没有显著差异。它输出的IncLevel（inclusion 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.bam，rMATS要靠坐标索引快速抓取特定区域的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 python=3.7 -y conda activate rmats conda install -c bioconda rmats=4.1.2 -y

这里特意指定了python=3.7，因为rMATS的pysam和numpy版本兼容性比较保守，在太新的Python版本下有时会出一些莫名其妙的报错。4.1.2是我测试过比较稳定的版本，官方后来发布的4.1.2之后更新不多，但大版本内部小版本之间的差异也需要注意。

装好后先验证一下：

rmats.py --help

如果能正常打印参数说明，说明基础环境没问题。我遇到过一个奇怪的情况是conda自动装了最新版rMATS，调用时提示No module named 'numpy'，这时候别急着重装，先看下环境里是否真的缺少numpy，如果缺就conda install numpy=1.21补上。

2.2 常用配套工具清单

除了rMATS本身，还需要准备samtools、STAR、gffread这些常用工具。samtools用于bam索引和格式检查，STAR用于上游比对，gffread则是在注释文件格式有问题时的救急工具。建议一并装好：

conda install -c bioconda samtools star gffread -y

2.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.2

3. 输入文件制备：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.bam

rMATS读取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.bam

b2.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 2>&1 &

运行过程中如果失败，首选看log文件末尾的报错信息。按我遇到的情况，高频问题依次是：GTF字段缺失、BAM坐标系和GTF不一致、样本文件路径出错、内存不够。前两个问题前面已经说过怎么处理，路径出错一般检查b1/b2文本文件有没有额外空格或者换行符，内存不够就减少线程数或加大机器内存配额。

5. 结果文件解读：别只盯着PValue，IncLevel才是核心

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后常见的列有：

ID，GeneID，geneSymbol，chr，strand，exonStart_0base，exonEnd upstreamES，upstreamEE，downstreamES，downstreamEE IncLevel1，IncLevel2，IncLevelDifference，PValue，FDR

这里的exonStart_0base和exonEnd是被检测的盒式外显子（cassette exon）的坐标，注意start是0-based的，传到IGV或UCSC浏览器里通常能直接用，但如果你是拿去做自定义可视化，记得做0-based和1-based之间的换算。

IncLevel1对应的是b1.txt这一组样本的平均inclusion level，IncLevel2对应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 pysam

6.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 跑了一半内存不足被killed

rMATS比较吃内存，尤其在同时加载大量bam索引信息的时候。降低--nthread并限制机器上的并发任务数量通常能解决。还可以尝试增加--tmp所在分区的磁盘空间，rMATS中间会写不少临时文件，磁盘满了也会被杀。用watch df -h看下磁盘占用，提前清理别的东西。

7.5 使用小规模测试数据快速验证流程

正式跑全基因组数据之前，我强烈建议先拿小规模数据把流程走通。可以只选一条染色体（比如chr21或者chr1的一部分）先比对、再跑rMATS，整个过程几分钟出结果，确认输出目录里能生成正常的事件文件，再去跑全量数据。这样可以帮你把环境问题、格式问题和路径问题在最早期暴露出来，省下后面数小时的等待。

8. 我在实际项目中的几点体会

做可变剪切真不是简简单单把命令跑通就行。整个流程下来，最花时间的往往不是rMATS本身的运行，而是前期的数据质控、格式检查以及后期的事件验证。尤其是BAM文件的质量，直接决定了后面所有结果的可靠性。如果你发现自己的结果和预期相差很远，我会建议先回头检查比对环节，而不是急着调整rMATS参数。

另外就是样本设计问题。rMATS虽然能处理两两比较，但生物学重复的数量和质量是最硬的条件。条件允许的情况下，每组至少三个重复，尽量让组内一致性高一点，否则统计这个环节会非常尴尬。

最后再分享一个我习惯用的工作流：跑完rMATS后，除了过滤显著事件，我还会把IncLevelDifference排名靠前但不显著的事件也留一个备份，有时候这些边缘事件反而是某些样本特异的、有意思的变化。后续实验验证时，这类事件偶尔会给你意外的惊喜。但写文章的时候，还是老老实实以显著事件为准，边缘信号只能当线索做探索，不能当真结论写。

这些经验都是拿真金白银的时间和失败的运行换来的。希望这份梳理能让你少踩一些我当年踩过的坑，让你把更多精力放在生物学问题的解读上。