RNA-seq 可变剪接分析与 PSI 解读
可变剪接指同一基因的前体 mRNA 选择不同外显子或剪接位点，产生多种成熟转录本。可变剪接分析比较两组样本的剪接模式，核心指标 PSI（percent spliced in）是某事件中包含型转录本的比例，组间差值记为 ΔPSI。rMATS 基于 BAM 中的剪接读段统计已注释和新事件；SUPPA2 基于转录本定量计算 PSI；leafcutter 不依赖注释，比较内含子切除比例。典型数据是深度较高的双端 bulk RNA-seq，常用于剪接因子敲低和肿瘤与正常组织比较。
什么时候用
- 剪接因子（如 SR 蛋白、hnRNP 家族）敲低或过表达后鉴定受调控事件
- 比较肿瘤与正常组织，寻找肿瘤特异的剪接异构体
- 基因表达量不变但功能改变，怀疑异构体比例发生切换
- 发现未注释的新剪接位点或内含子保留事件
- 做剪接 QTL（sQTL）前计算各样本的内含子使用比例
不适合的情况
- 读长短、测序深度低的数据：剪接读段太少，PSI 估计方差很大
- 3′ 端测序和多数液滴法单细胞数据：读段基本不覆盖剪接位点
- 需要完整异构体结构：短读长只能看局部事件，全长结构需长读长测序
需要准备的数据
- BAM 文件
- rMATS、leafcutter 需要 STAR 或 HISAT2 比对的排序 BAM，各样本使用同一基因组和比对参数
- GTF 注释
- 定义已知事件，须与比对时的版本一致；rMATS 会在注释基础上结合读段发现新剪接位点
- 转录本定量结果
- SUPPA2 需要 Salmon 或 kallisto 的转录本 TPM，转录本 ID 与生成事件所用 GTF 一致
- 分组与文库信息
- 每组样本列表、读长、单端或双端、链方向，rMATS 需要在命令中显式指定
分析步骤
- 1
评估数据是否适合
确认读长、双端与否和测序深度，查看 STAR 的 SJ.out.tab 中剪接读段数量。每组通常建议至少 3 个生物学重复。用 IGV 抽查几个已知剪接基因，确认剪接读段清晰可见。
- 2
统一比对
所有样本用相同参数比对，STAR 可用两遍模式（--twopassMode Basic）提高新剪接位点检出。各样本读长应一致；数据经过修剪、长度不一时，后续 rMATS 需允许可变读长。
- 3
运行事件检测
rMATS 输入两组 BAM、GTF、读长和 libType，输出 SE、A5SS、A3SS、MXE、RI 五类事件。SUPPA2 由 GTF 生成事件后用 diffSplice 比较。leafcutter 从剪接读段聚类内含子。
- 4
选择计数方式
rMATS 同时给出 JC（仅跨剪接位点读段）和 JCEC（再加外显子内部读段）两套结果。JC 直接反映剪接方式，常作为主结果；内含子保留事件的读段多落在内含子内，可参考 JCEC。
- 5
筛选显著事件
按 FDR 和 |ΔPSI| 双重阈值筛选，同时要求包含和跳跃读段之和达到最低支持数，去掉由少量读段决定的事件。统计各类事件数量，外显子跳跃（SE）通常占多数。
- 6
可视化与功能解读
用 rmats2sashimiplot 或 ggsashimi 绘制 sashimi 图展示各组剪接读段。对受影响基因做富集分析，并检查事件是否造成移码、引入提前终止密码子或改变蛋白结构域。
- 7
实验验证
挑选 ΔPSI 大、读段支持充分的事件，在两侧组成型外显子上设计引物做 RT-PCR。凝胶上包含型与跳跃型条带的比例变化方向，应与测序得到的 ΔPSI 一致。
关键参数与决策
ΔPSI 阈值
常用 |ΔPSI| > 0.1 作为经验阈值，结合 FDR < 0.05。剪接因子敲低等强扰动可提高到 0.2 以减少假阳性；ΔPSI 很小的事件即使显著，RT-PCR 也很难验证出来。
最低读段支持
建议要求每个样本中该事件的包含与跳跃读段之和达到 10 条左右（经验值），或至少在大多数样本中满足。低覆盖事件的 PSI 可能只由一两条读段决定，波动极大。
读长与可变读长
rMATS 默认只使用与 --readLength 一致的读段。数据经修剪后长度不一时需加 --variable-read-length，否则大量读段被丢弃，检出事件明显变少却不容易察觉。
libType
--libType 需与文库一致：非链特异用 fr-unstranded，dUTP 法链特异文库用 fr-firststrand。设错会影响反义重叠基因区域的读段归属。
对注释的依赖
rMATS 和 SUPPA2 以注释为主，结论偏向已知事件；关注未注释剪接或涉及多个外显子的复杂事件时，可选 leafcutter 或 MAJIQ 这类不预设事件类型的方法。
常见错误
- 把基因表达变化当剪接变化：基因整体下调不改变 PSI，但低表达会让 PSI 估计不稳，筛选时要同时看基因表达量和读段支持数。
- 两组样本读长或测序批次不同：读长影响剪接读段的检出概率，会造成系统性 PSI 差异，合并不同来源的公共数据时尤其常见。
- 只按 FDR 排序挑事件：高表达基因上极小的 ΔPSI 也能得到很小的 FDR，生物学意义有限，应同时考虑 ΔPSI 大小和读段数。
- 直接比较不同工具的事件列表：rMATS 与 SUPPA2 对同一事件的坐标定义和类型名称不同，未统一坐标就求交集会低估一致性。
- 忽视内含子保留的技术来源：核内未完成剪接的前体 RNA、总 RNA 文库和 RNA 降解都会抬高内含子读段，RI 事件需谨慎解读。
- 弄反 IncLevelDifference 的方向：rMATS 中该值是第一组均值减第二组均值，b1、b2 的输入顺序决定正负，解释前要核对。
结果怎么看
事件类型统计表
按五类事件列出显著事件数和 PSI 升高、降低的数量。若 RI 事件远多于其他类型，先检查文库类型和核 RNA 比例，再考虑是否为真实的生物学调控。
差异事件明细表
rMATS 输出 IJC、SJC 读段数，IncLevel1、IncLevel2（每个样本的 PSI）、IncLevelDifference 和 FDR。看同组样本 PSI 是否接近，组内离散很大说明该事件结论不稳。
sashimi 图
展示外显子覆盖度和弧线标注的剪接读段数。显著事件应能看到组间弧线读段比例的明显变化；如果弧线上的读段只有个位数，结论需要更多证据支持。
PSI 聚类热图
用显著事件的 PSI 对样本聚类，同组样本应聚在一起。若样本按测序批次而非分组聚类，说明剪接结果受技术因素影响较大，需要回查批次设计。
和相近方法怎么选
rMATS 与 SUPPA2
rMATS 直接统计 BAM 中的剪接读段，结果可追溯到具体读段，适合中高深度数据；SUPPA2 基于转录本定量，速度快，可复用 Salmon 结果，但依赖注释完整性，转录本定量误差会传递到 PSI。
leafcutter 与基于事件的方法
leafcutter 不预设事件类型，以内含子簇为单位检验切除比例变化，适合发现未注释剪接和复杂事件，也常用于 sQTL；缺点是结果不直接对应 SE、RI 这类直观类型，解释需要额外注释。
剪接事件与异构体切换分析
DEXSeq 检验外显子使用差异，IsoformSwitchAnalyzeR 关注完整异构体占比切换并预测功能后果。关心局部剪接事件用 rMATS 一类；关心主导异构体及其蛋白功能变化用异构体切换分析。
常见问题
PSI 是什么意思，怎么计算？
PSI 是某剪接事件中包含型的比例。以外显子跳跃为例，PSI 等于包含读段数除以包含与跳跃读段数之和，两类读段都按各自的有效长度归一。PSI 为 0.8 表示约 80% 的转录本含该外显子。
rMATS 用 JC 还是 JCEC 的结果？
多数研究以 JC 为主，因为跨剪接位点的读段直接证明剪接方式；JCEC 额外计入外显子内部读段，统计功效更高，但会受外显子自身表达影响。论文中写明使用哪一套，关键事件可两套都查看。
可变剪接分析需要多少测序深度？
没有统一门槛。剪接分析需要的深度明显高于基因表达定量，经验上常建议双端、读长 100 bp 以上并保证较高的读段数。深度有限时，可只分析高表达基因上的事件，并提高最低读段支持要求。
单细胞数据能做可变剪接分析吗？
10x 等液滴法的 3′ 或 5′ 测序读段集中在转录本一端，基本不能做常规剪接分析。Smart-seq2 等全长方法可以，但单个细胞读段少，通常按细胞类型合并后再计算 PSI。
显著事件太多，怎么挑重点？
先按 |ΔPSI|、读段支持数和组内一致性收紧过滤，再优先看造成移码、引入提前终止密码子或影响已知功能域的事件，以及在独立数据集中也能重复的事件，然后再安排 RT-PCR 验证。
在 Codar 里怎么做
Codar 先与你确认读长、文库链方向、FDR 与 ΔPSI 阈值和最低读段支持数，再读取已有的比对 BAM 运行可变剪接分析，统计各类事件，计算样本 PSI 和组间 ΔPSI 并做显著性筛选。交付事件类型统计、显著事件明细表、PSI 聚类图、可复现脚本和方法说明。
我有 SRSF1 敲低和对照各 3 个样本的 STAR 比对 BAM，双端 150 bp，请用 rMATS 找差异剪接事件，ΔPSI 阈值设 0.1