企业官网建设流程全解析

搞转录组数据分析的同行应该都有这种感觉：RNA-seq的常规差异表达分析早就是流水线作业了，但一到可变剪切（Alternative Splicing）层面，很多流程就开始变得不那么“顺手”。可变剪切在真核生物的基因表达调控里太重要了，同一基因通过不同的外显子组合方式产生多种转录本异构体，这些异构体在发育、疾病、药物响应里扮演的角色，往往是普通差异表达分析根本看不出来的。

rMATS是我在实际项目中反复对比之后留下来的工具。它在可变剪切分析领域的地位有点像差异表达里的DESeq2，模型严谨、输出规范、社区用户多，遇到问题基本都能搜到解决方案。这篇文章我就把最近一轮完整跑通的rMATS分析流程做个系统记录，从Linux环境准备、工具安装、BAM和GTF输入文件的准备，到rMATS运行参数调优、五类剪切事件的结果解读、显著事件的筛选、可视化作图，再加上这一路踩过的那些官方文档里不会写的坑，全部按实际执行的顺序拆开来讲。无论你是刚接触生信的新手，还是已经跑过一些流程但卡在结果解读上的人，这篇都能当一份能直接复用的参考。

1. rMATS到底在做什么：五类剪切事件与统计逻辑

1.1 从实际问题理解可变剪切分析的目标

先明确一个容易被新手混淆的点：可变剪切分析回答的，和差异表达分析回答的，完全是两个不同层面的问题。

假设你有两组样本，对照组3个生物学重复，药物处理组3个生物学重复。差异表达分析告诉你的是：某个基因的整体表达量在两组之间有没有显著变化。但一个基因可能产生A、B、C三种转录本，总表达量没变，可B型转录本的比例从20%飙到了80%，这种变化在基因水平上会被完全掩盖。可变剪切分析就是专门盯这种“转录本构成比例的变化”。

rMATS把这种比例叫做inclusion level，也就是某个可变剪切事件中，包含特定外显子或剪接位点的转录本占该基因所有转录本的比例，通常用Ψ（psi）表示。两组之间的Ψ差异，就是可变剪切分析最核心的输出。理解了这一点，后面看结果文件里的IncLevel1、IncLevel2这些列就顺理成章了。

1.2 rMATS的统计模型是怎么工作的

rMATS的算法核心可以概括为：对每个可变剪切事件，分别统计两组样本中支持“包含型”（inclusion）和“跳过型”（skipping）的reads数量，然后对两组样本的inclusion level差异做统计检验。

在read计数上，rMATS把比对到RNA-seq数据的reads分成两类：

• Junction reads：跨越剪接位点的reads，能提供剪接连接的直接证据。
• Exon reads：落在外显子内部的reads，提供外显子表达水平的证据。

rMATS默认只使用junction reads进行计算，这就是JC（Junction Count）结果文件；如果加上exon reads一起计算，就叫JCEC（Junction Count and Exon Count）结果文件。JC结果更保守、假阳性更低，JCEC由于利用了更多信息，在某些情况下能检测到更多事件，但也可能引入更多由reads分布不均造成的噪音。

统计检验方面，rMATS对每个事件构建了一个似然比检验（likelihood ratio test），比较“两组inclusion level相同”和“两组inclusion level不同”这两个模型哪个更符合观测到的reads分布。它把生物学重复之间的变异也纳入模型，这也是rMATS相比早期一些简单比率法工具的优势所在。

1.3 五种可变剪切事件的识别与对应文件

rMATS一共分析五类可变剪切事件，这也是它的核心输出维度。

每种事件在GTF注释里都有特定的“结构特征”。比如SE事件需要找到一段被上下游外显子夹在中间、可能被跳过的外显子；MXE事件则是两个相邻的、通常不会同时出现在同一个转录本里的外显子；RI事件的关键是内含子区域有没有足够的reads覆盖证据。rMATS会根据GTF注释把所有可能的候选事件枚举出来，然后逐一进行reads计数和统计检验。

这个事件类型的设计直接决定了你后续分析的口径。比如你在研究某个神经发育相关基因，文献里反复提到它的exon 6跳跃（SE事件），那你就应该重点看SE.MATS.JC.txt里这个基因对应的那一行，而不是去RI文件里找。理解五类事件的区别，是看懂结果的第一步。

2. Linux环境准备与rMATS安装避坑

2.1 建议的Linux发行版与Python环境隔离方案

rMATS对Linux发行版本身没有特别挑剔的要求，CentOS 7、Ubuntu 18.04/20.04、Debian都跑过。真正让你头疼的往往不是系统本身，而是Python环境混乱。

rMATS turbo版（目前的主流版本）基于Python 3开发，但它的很多核心模块依赖旧版本库的编译行为。我强烈建议你用miniconda创建一个独立的环境，不要直接装在系统自带的Python里，更不要装在root的全局site-packages里。生信项目多、工具杂，每个工具对Python版本的诉求都不一样，环境隔离必须从第一天就养成习惯。

# 下载并安装miniconda wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh

安装完conda之后，创建rMATS专用环境。我自己用的组合是Python 3.8，这个版本在rMATS上表现最稳。

2.2 rMATS安装的三种方式与踩坑对比

rMATS turbo的安装方式主要有三种：conda安装、pip安装、源码编译。我实践下来，pip安装最省心，conda次之，源码编译适合需要改源码或者对版本有特殊要求的情况。

方式一：pip安装（推荐）

conda create -n rmats python=3.8 -y conda activate rmats pip install rmats

安装完成后验证一下：

rmats.py --version

如果输出版本号，说明安装成功。rMATS turbo的版本号一般类似v4.1.2或者v4.2.0。

方式二：conda安装

conda create -n rmats -c bioconda -c conda-forge rmats

conda安装的好处是依赖关系处理得比较干净，缺点是bioconda源有时更新不及时，拿到的版本可能不是最新的。

方式三：源码编译

git clone https://github.com/Xinglab/rmats-turbo.git cd rmats-turbo python setup.py install

源码编译的坑在于对编译器和python-dev头文件的要求比较严格，如果系统缺少gcc或者Python的头文件，编译会直接报错。

这里有一个我踩过的实际坑：最初我在Python 3.10环境里用pip装了rmats，命令行能正常启动，但跑到统计阶段直接报错，报错信息指向一个隐式类型转换的问题。后来把环境降级到Python 3.8重新安装，问题就消失了。所以如果你遇到莫名其妙的运行错误，第一时间检查Python版本。

2.3 配套工具：samtools与bam处理

rMATS本身不负责序列比对，但它的输入是bam文件，而且必须是按坐标排序并建立索引的bam。所以samtools是实打实的必备工具。

conda install -c bioconda samtools -y

比对工具方面，STAR是目前与rMATS搭配最主流的RNA-seq比对工具。STAR比对速度快、junction支持好，而且可以直接输出按坐标排序的bam。我个人的习惯是，拿到fastq后的第一步处理就是STAR比对加samtools排序索引，一套组合拳打完，直接就能喂给rMATS。

3. 输入数据准备：BAM文件与GTF注释的硬性要求

3.1 BAM文件的质量检查和格式规范

rMATS跑出来的结果质量，很大程度上在输入bam阶段就已经决定了。这个环节有三个容易被忽略的硬性要求。

第一，bam必须按坐标排序（coordinate-sorted），而不是按read名称排序。很多比工具默认输出按坐标排序的bam，比如STAR加--outSAMtype BAM SortedByCoordinate参数后得到的bam，但如果你用的是其他工具或自定义管道，一定要用samtools确认一下：

samtools view -H sample.bam | grep SO # 期望看到 SO:coordinate

第二，bam必须有索引文件。没有bai文件的bam，rMATS在读取时要么报错，要么长时间卡在某个CTGene的读取环节，非常影响排查体验。用samtools index一行搞定：

samtools index sample.bam

第三，生物学重复不可省。rMATS的统计模型需要利用样本间变异来估算离散度参数，如果你只有单个样本和单个样本做比较，rMATS也能跑，但统计效力非常弱，结果基本不可信。所以我一般建议无论实验设计怎么紧张，至少保证每组有两个以上的生物学重复。

3.2 GTF注释文件选对版本是关键中的关键

GTF文件的选择直接决定rMATS能识别出哪些可变剪切事件。我的建议是优先从GENCODE或Ensembl官网下载，并且必须与你比对时用的参考基因组版本严格匹配。

这个“版本一致”的重要性怎么说都不为过。比如你用GRCh38的参考基因组做STAR比对，GTF却用了hg19的版本，染色体命名可能都是chr1和1这种对不上的格式，rMATS跑出来的结果会有一大堆NA，甚至直接报错。另外GTF的时间版本也值得注意，GENCODE会定期更新注释，相同基因组的不同release版本之间差异不小，新旧版本混用可能导致某些基因的事件数明显异常。记录下你用的GTF具体版本号，写论文或报告时也需要用到。

3.3 比对参数设置影响junction reads的识别

rMATS事件检测依赖bam中的junction reads，因此比对这一步对最终结果的影响非常大。

我常用STAR比对，关键的参数设置如下：

STAR --runThreadN 16 \ --genomeDir /path/to/star_index \ --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outSAMunmapped Within \ --outSAMattributes NH HI AS NM MD

有几个细节需要留意：

• --outSAMtype BAM SortedByCoordinate一步到位，得到排序bam。
• --outSAMattributes里的NH（number of hits）属性对rMATS去多比对reads很重要，最好不要省。
• 如果后续还要跑其他工具比如RSEM、StringTie，可以参考STAR手册加上--quantMode，但单纯为rMATS准备输入时不需要。
• 双端测序时，--outSAMprimaryFlag AllBestScore这类参数会影响read的primary标记，统一使用STAR默认设置即可，不要随意调。

建索引时还有一点，STAR的基因组索引要和GTF配套。STAR建索引时使用--sjdbGTFfile参数加入GTF，这能显著提升junction reads的比对准确性。

4. rMATS实战运行：命令、参数与输出

4.1 样本列表文件的准备格式

rMATS通过两个文本文件来指定对照组和处理组的bam路径。每个文件每行一个bam文件的绝对路径，没有任何额外格式要求。

# b1.txt（对照组） /home/project/data/ctrl_rep1.bam /home/project/data/ctrl_rep2.bam /home/project/data/ctrl_rep3.bam # b2.txt（处理组） /home/project/data/treat_rep1.bam /home/project/data/treat_rep2.bam /home/project/data/treat_rep3.bam

这里有个非常容易踩的坑：一定要用绝对路径。我自己就曾经在脚本里用了相对路径，rMATS启动时的当前工作目录一变，程序直接报找不到文件，排查了半天才发现是路径问题。另外，bam文件名里不要带空格或特殊符号，这会干扰rMATS内部的shell调用来处理。

4.2 运行命令逐参数详解

rMATS的完整运行命令如下：

rmats.py --b1 b1.txt --b2 b2.txt \ --gtf gencode.v38.annotation.gtf \ -t paired --readLength 150 \ --nthread 8 --od output --tmp tmp

这个命令里的每个参数都值得展开讲清楚：

• --b1和--b2：指定两组样本的bam列表文件。
• --gtf：GTF注释文件的路径。
• -t：测序类型，paired表示双端，single表示单端。这个参数千万别设错，错了直接影响比对事件的检测。
• --readLength：测序读长。双端150bp就写150。rMATS会用它来校正一些read长度相关的计数偏差。
• --nthread：线程数。我一般按物理核心数的一半设置，8到16之间比较常见。线程设太高时CPU上下文切换开销明显，加速效果反而下降。
• --od：输出目录，注意要预先创建好，rMATS不会自动创建。
• --tmp：临时文件目录，rMATS运行时会产生大量中间文件，建议专门指定一个目录，方便运行结束后清理。

还有几个进阶参数，按需使用：

• --statoff：跳过统计检验阶段，只做事件检测和reads计数。当数据量特别大、想做预筛选时很有用。
• --cstat：指定统计量类型，0表示输出P值，1表示输出FDR（默认是1）。
• --readLength和--libType配合使用可以更精细控制文库类型，不过绝大多数场景用默认设置就够了。

4.3 资源消耗与运行监控

rMATS的运行时间受数据量、参考基因组大小、事件候选数量和线程数共同影响。以3对3的人类转录组样本、每个样本约30M对reads为例，8线程环境下通常需要跑3到6小时。整个运行过程比较“安静”，没有实时进度条，但你可以通过查看tmp目录里的临时文件增长来确认程序确实在干活。

# 查看内存和CPU占用 top -u $USER # 查看临时目录大小变化 du -sh tmp/

一个容易被忽略的资源参数是磁盘空间。rMATS的tmp目录在运行期间会写入大量中间文件，有些项目能占到几十GB。我建议在运行前用df -h确认磁盘剩余空间在50GB以上，否则跑到一半磁盘满了，整个任务直接崩溃，前功尽弃。

4.4 输出目录结构解析

rMATS运行完成后，输出目录里会生成如下文件：

output/ ├── 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 ├── fromGTF.SE.txt ├── fromGTF.A5SS.txt ├── ... ├── JC.raw.input.txt └── readCount/

fromGTF.*.txt系列文件是从GTF注释中鉴定出的所有候选事件，不管有没有reads支持都列在里面；readCount目录里则是每个事件的具体reads计数矩阵。JC和JCEC文件就是我们做下游分析的主要对象。

5. 结果深度解析：从统计学显著到生物学意义

5.1 SE.MATS.JC.txt字段逐列解读

以最核心的SE.MATS.JC.txt为例，展开每个字段的含义。

其中的IncLevel这一列最能直观体现可变剪切变化的方向和幅度。比如IncLevel1是0.85,0.82,0.88，IncLevel2是0.20,0.25,0.18，说明这个外显子本来在对照组里大部分转录本都包含，到了处理组几乎全被跳过了，差异非常大。IncLevelDifference就是-0.65左右，负号代表处理组的inclusion水平降低。

5.2 筛选显著事件的实操策略

拿到结果文件后，筛选显著差异的可变剪切事件，我一般分三步走。

第一步，过滤掉统计检验不稳定的行。rMATS结果中PValue和FDR列有时会出现NA，这些通常是因为事件的reads数太少，统计检验不可靠，直接过滤掉。

# 以SE.MATS.JC.txt为例，过滤掉PValue或FDR为NA的行 head -1 SE.MATS.JC.txt > SE.filtered.txt awk 'NR>1 && $20!="NA" && $21!="NA"' SE.MATS.JC.txt >> SE.filtered.txt

第二步，按阈值筛选显著事件。常用的组合是FDR < 0.05且|IncLevelDifference| > 0.1。如果是做初筛，可以放宽到0.05；如果要锁定核心候选事件做实验验证，我建议收紧到|IncLevelDifference| > 0.2，这个标准的生物学差异更明显。

# 筛选FDR < 0.05且差异绝对值大于0.2的事件（SE文件中FDR是第21列，IncLevelDifference是第22列） awk 'NR==1 || ($21<0.05 && $22>0.2) || ($21<0.05 && $22<-0.2)' SE.filtered.txt > SE.significant.txt

第三步，检查IncLevel本身的数值范围。如果一个事件在两组中都是极端值，比如IncLevel1在0.98左右，IncLevel2在0.95，虽然统计上可能显著，但生物学意义的解读要谨慎。反过来，如果IncLevel从0.5跳到0.1，这种变化对蛋白功能的影响通常远超统计数字本身。

5.3 可视化：sashimi plot与IGV验证

筛选出显著事件后，可视化验证是必不可少的。最常用的是sashimi plot，它能直观展示两组样本在某个基因位点的reads覆盖和剪接连接情况。rMATS官方配套工具rmats2sashimiplot可以很方便地画。

# 安装rmats2sashimiplot pip install rmats2sashimiplot # 画sashimi图 rmats2sashimiplot --b1 ctrl_rep1.bam,ctrl_rep2.bam,ctrl_rep3.bam \ --b2 treat_rep1.bam,treat_rep2.bam,treat_rep3.bam \ -e SE.MATS.JC.txt \ --event_type SE \ --l1 Control --l2 Treatment \ -o sashimi_output

sashimi图上半部分是基因结构示意图，下半部分是两组样本的reads覆盖曲线，弧线代表junction reads的跨外显子连接。某条弧线的厚度直接反映该剪接连接的相对丰度。通过sashimi图能很直观地看到，处理组中支持外显子跳跃的junction弧线明显变粗，而支持包含型的弧线变细甚至消失。

除了自动化绘图，我强烈建议对最终候选事件逐个用IGV做人工检查。IGV加载bam文件和GTF注释，手动查看目标基因位点的reads分布，可以排除一些软件自动判定中的假阳性。比如某个“显著的SE事件”可能实际上被附近的重复区域或假基因干扰，人工检查一眼就能发现。

5.4 从定量结果到功能解读

统计显著的可变剪切事件，后续怎么和生物学功能挂钩，是这个环节最有意思也最考验功底的部分。

• 如果目标基因是转录因子或信号通路关键分子，检索文献看这个外显子是否编码蛋白的功能结构域，比如DNA结合域、跨膜域、激酶活性位点等。
• 用生物信息学工具预测剪接异构体是否可能触发无义介导的mRNA降解（NMD），这决定了剪接变化是“改变功能”还是“直接降低表达”。
• 把显著事件涉及的基因做GO/KEGG富集分析，看是否集中在某些特定通路里，比如凋亡、免疫应答、细胞周期等。
• 和大规模公共数据库交叉验证，比如TCGA中同一基因在肿瘤样本里是否存在类似的剪接异常。

5.5 另外一个容易忽略的检验思路

rMATS是基于参考基因组的注释来枚举事件的，因此它的检测范围受限于GTF注释的完整度。如果关心的物种注释质量很差，可以考虑先做一次转录本组装（比如用StringTie）来补充注释，再跑rMATS。不过这样做要谨慎，组装出来的转录本可靠性参差不齐，容易引入大量假阳性事件，建议还是优先使用成熟物种的高质量注释。

6. 常见问题与排查技巧实录

6.1 安装与环境问题速查

6.2 运行阶段的高频报错

运行中最常见的报错集中在BAM文件和GTF的兼容性上。

报错一：染色体命名不一致

GTF里的染色体名是chr1，但bam里的参考序列名是1，或者反过来。rMATS会找不到对应染色体的注释，报错信息可能不那么直观，但结果通常是一大堆事件检测失败。解决办法是用samtools给bam重新头：

# 查看bam里的染色体名 samtools view -H sample.bam | grep @SQ # 查看GTF前几行的染色体名 head -5 annotation.gtf

如果两者确实不一致，可以用samtools reheader或者awk批量替换GTF里的染色体名，统一成同一套命名。

报错二：GTF文件有未排序的记录

rMATS要求GTF按染色体和起始位置排序，如果不满足，它会直接抛出格式错误。用下面的命令对GTF排序：

# 先按染色体，再按起始位置排序 sort -k1,1 -k4,4n annotation.gtf > annotation.sorted.gtf

报错三：tmp目录写入失败或磁盘满

No space left on device是最直接的提示，但有时rMATS会卡在某个文件写入阶段不退出。排查方法就是定期用df -h看磁盘容量，给tmp目录留够空间。

6.3 结果异常的处理思路

结果几乎全是NA怎么办？

首先检查bam和GTF的染色体命名是否一致、基因组版本是否匹配；其次检查bam里有没有足够的junction reads，用samtools view看一眼bam大小和reads量，如果比对率极低或数据量太小，rMATS自然统计不出有意义的结果。

两组样本的IncLevel差异不显著但生物学表型很明确？

这时候要考虑是不是测序深度不够。可变剪切检测对测序深度的要求比基因表达定量高得多，尤其是那些表达量中低等的基因，junction reads本来就少，统计检验的效力跟不上。一个可行的补救办法是提高重复数，4个或5个生物学重复能明显提升rMATS检出力。另外也可以尝试JCEC文件，它在某些场景下能利用exon reads的信息来补充junction reads的不足。

筛选出来的显著事件特别多或特别少？

特别多的情况多半是FDR阈值放得太宽，或者两组样本之间的整体剪接模式差异本身就很大。特别少的情况除了数据质量问题，也可能是GTF注释不全、事件候选数量太少。对比一下fromGTF文件里的候选事件总数，如果连几千个都没有，说明注释确实存在瓶颈。

6.4 实际项目中值得保留的几个习惯

跑rMATS跑多了，我慢慢形成了一套固定的操作习惯，在这里分享给大家参考。

• 每次运行前把版本信息、GTF版本、参数命令完整记录到一个文本文件里，方便论文方法部分直接引用。
• 最终筛选出的显著事件列表，一定要导出成Excel，同时保留IncLevel的原始值和差异值，方便后续做热图展示和给合作者查看。
• 对所有下游分析输入文件（筛选后的事件列表、BAM路径等）都用绝对路径，脚本里不要用~或者相对路径。
• 在跑大样本量的rMATS任务之前，先挑一个染色体做小规模测试，确认参数和输入文件都没问题，再全量跑，能省去很多返工时间。

6.5 关于重复样本数量的一点补充

rMATS的统计模型对重复数的要求比普通差异表达分析更高。2个重复能跑，但结果方差估计通常不太稳定；3个重复是比较常见的配置，大多数文章都认可；4个及以上会更稳，但边际收益递减。如果你的数据是混合样本池（pooled sample），每个池里有多个个体的RNA混合，那么rMATS的统计检验可能不够合理，因为混合样本丢失了生物学重复之间的变异信息，P值会偏向显著，所有事件看起来都在差异表达。这种数据最好换用其他工具或者重新设计实验。

我个人在实际项目中还有一个心得：rMATS的结果只是候选列表，它告诉我们“哪些事件可能变了”，但到底“这个剪接变化对细胞功能有没有影响”，必须回到具体的基因和生物学背景中去验证。筛选完事件之后，挑两三个关键基因做一下qRT-PCR验证或者半定量PCR跑胶，成本不高，但对结论的说服力提升是决定性的。剪接水平的验证引物设计技巧也不少，感兴趣的话以后可以再写一篇详细展开。