Pipeline Detail

RNA-Seq转录组结构与网络

可变剪接分析 rMATS/SUPPA2

面向 bulk RNA-seq 的可变剪接增强流程，覆盖实验设计、剪接事件类型、STAR 比对、rMATS 事件检测、SUPPA2 PSI/dPSI 分析、可视化和候选事件解释。

创建时间

2026/6/3

分析难度

中级

推荐场景

转录组分析

预计耗时

1-3 天

Metadata

流程元数据

先看应用场景、输入输出和工具依赖，再进入正文命令细节。

Difficulty

中级

Scenario

转录组分析

Estimated Time

1-3 天

Tools

STARSalmonrMATSSUPPA2

Inputs

FASTQBAMGTFTPM

Outputs

Sashimi plotreport

Workflow DAG

流程图

用步骤节点快速理解这个分析从原始数据到结果报告的流转关系。

STEP 1

实验设计与分组文件

→

STEP 2

FASTQ 质控与清洗

→

STEP 3

STAR 剪接感知比对

→

STEP 4

BAM 与 junction QC

→

STEP 5

rMATS 事件级差异剪接

→

STEP 6

Salmon 转录本定量

→

STEP 7

SUPPA2 PSI/dPSI 分析

→

STEP 8

事件过滤与交叉验证

→

STEP 9

Sashimi/PSI/火山图

→

STEP 10

候选事件报告

Protocol

流程文档

正文保留 Markdown 排版、代码语言标识和表格样式，适合边学边复现。

可变剪接分析 rMATS/SUPPA2

一、适用场景

可变剪接分析关注同一个基因在不同条件下是否使用了不同的外显子、剪接位点或转录本结构。它回答的问题不是“基因整体表达是否上调/下调”，而是“基因内部 isoform 使用比例是否改变”。

适用于：

处理组 vs 对照组的外显子跳跃、内含子滞留、可变 5'/3' 剪接位点分析。
疾病、突变体、胁迫处理、发育阶段中的 splicing regulation 研究。
RNA-binding protein 敲除/过表达后的剪接调控事件发现。
与差异表达、转录因子、RBP motif、蛋白结构域变化联合解释。

前提建议：

每组至少 3 个生物学重复。
优先使用 paired-end RNA-seq。
reads 长度建议 100 bp 或以上，越长越利于跨 junction 事件识别。
rMATS 路线需要高质量 BAM；SUPPA2 路线需要可靠转录本注释和表达定量。

二、可变剪接事件类型

事件类型	英文缩写	生物学含义	rMATS 输出
外显子跳跃	SE	一个 cassette exon 在两组间 inclusion 水平改变	`SE.MATS.JC.txt`
可变 5' 剪接位点	A5SS	使用不同 5' splice site	`A5SS.MATS.JC.txt`
可变 3' 剪接位点	A3SS	使用不同 3' splice site	`A3SS.MATS.JC.txt`
互斥外显子	MXE	两个外显子互斥使用	`MXE.MATS.JC.txt`
内含子滞留	RI	intron 在成熟 RNA 中保留	`RI.MATS.JC.txt`
转录本 isoform 切换	ISO	主导转录本比例改变	SUPPA2 transcript-level PSI

三、最佳实践路线

建议同时准备两条互补路线：

STAR + rMATS：以 junction reads 和 exon body reads 为核心，适合标准事件级差异剪接。
Salmon + SUPPA2：以 transcript abundance 和注释事件为核心，适合快速计算 PSI 和 dPSI，也适合 isoform-level 解释。

两者互补：

工具	输入	优势	注意点
rMATS	sorted BAM 或 FASTQ + GTF	事件类型清晰，输出 inclusion level 和 FDR	依赖 junction 支持，BAM 和 read length 要准确
SUPPA2	GTF + transcript TPM	快速、适合多样本 PSI/dPSI 和 transcript-level 分析	依赖转录本注释和定量准确性

四、整体流程图

flowchart TD
    A[样本设计与分组] --> B[FASTQ 质控与清洗]
    B --> C[STAR two-pass 或剪接感知比对]
    C --> D[BAM sort/index 与 junction QC]
    D --> E[rMATS 检测 SE/A5SS/A3SS/MXE/RI]
    B --> F[Salmon 转录本定量]
    F --> G[SUPPA2 generateEvents]
    G --> H[SUPPA2 psiPerEvent]
    H --> I[SUPPA2 diffSplice]
    E --> J[显著事件过滤]
    I --> J
    J --> K[Sashimi plot / PSI 箱线图 / dPSI 火山图]
    K --> L[候选基因与事件报告]

五、项目目录建议

alternative_splicing_project/
├── 00_metadata/
│   ├── sample_info.csv
│   ├── group1_bams.txt
│   └── group2_bams.txt
├── 01_raw_data/
├── 02_clean_data/
├── 03_star/
│   ├── bam/
│   └── logs/
├── 04_rmats/
│   ├── output/
│   └── tmp/
├── 05_salmon/
├── 06_suppa2/
│   ├── events/
│   ├── psi/
│   └── diff/
├── 07_visualization/
└── report/

sample_info.csv 示例：

sample_id	group	batch	read1	read2
Ctrl_1	Ctrl	B1	Ctrl_1_R1.fq.gz	Ctrl_1_R2.fq.gz
Ctrl_2	Ctrl	B1	Ctrl_2_R1.fq.gz	Ctrl_2_R2.fq.gz
Treat_1	Treat	B1	Treat_1_R1.fq.gz	Treat_1_R2.fq.gz
Treat_2	Treat	B1	Treat_2_R1.fq.gz	Treat_2_R2.fq.gz

六、FASTQ 质控与清洗

mkdir -p 02_qc_raw/fastqc 02_qc_raw/multiqc

fastqc -t 8   -o 02_qc_raw/fastqc   01_raw_data/*.fq.gz

multiqc 02_qc_raw/fastqc   -o 02_qc_raw/multiqc

mkdir -p 02_clean_data 02_qc_clean/fastp

for sample in Ctrl_1 Ctrl_2 Ctrl_3 Treat_1 Treat_2 Treat_3
do
  fastp     -i 01_raw_data/${sample}_R1.fq.gz     -I 01_raw_data/${sample}_R2.fq.gz     -o 02_clean_data/${sample}_R1.clean.fq.gz     -O 02_clean_data/${sample}_R2.clean.fq.gz     --detect_adapter_for_pe     --length_required 30     --thread 8     --html 02_qc_clean/fastp/${sample}.html     --json 02_qc_clean/fastp/${sample}.json
done

质控重点：

每个样本总 reads 数是否接近。
接头污染是否明显。
GC 分布是否异常。
paired-end 两端质量是否一致。
可变剪接项目要特别关注 read length，因为 junction-spanning reads 对事件识别很关键。

七、STAR 剪接感知比对

1. 构建索引

mkdir -p 03_star/index

STAR   --runThreadN 16   --runMode genomeGenerate   --genomeDir 03_star/index   --genomeFastaFiles ref/genome.fa   --sjdbGTFfile ref/genes.gtf   --sjdbOverhang 149

sjdbOverhang 通常设置为 read length - 1，例如 PE150 设置为 149。

2. 比对并输出 sorted BAM

mkdir -p 03_star/bam 03_star/logs

for sample in Ctrl_1 Ctrl_2 Ctrl_3 Treat_1 Treat_2 Treat_3
do
  STAR     --runThreadN 16     --genomeDir 03_star/index     --readFilesIn       02_clean_data/${sample}_R1.clean.fq.gz       02_clean_data/${sample}_R2.clean.fq.gz     --readFilesCommand zcat     --outFileNamePrefix 03_star/bam/${sample}_     --outSAMtype BAM SortedByCoordinate     --outSAMstrandField intronMotif     --twopassMode Basic

  samtools index 03_star/bam/${sample}_Aligned.sortedByCoord.out.bam
done

建议：

rMATS 可从 BAM 输入，BAM 应统一使用同一参考基因组和同一 GTF。
如果后续使用 IGV/Sashimi plot，保留 indexed BAM 很重要。
--twopassMode Basic 有助于发现 novel junction，但项目中要记录是否开启，保证批次一致。

八、rMATS 差异剪接分析

1. 准备 BAM 分组文件

rMATS 的 --b1 和 --b2 接收逗号分隔的 BAM 路径。

ls 03_star/bam/Ctrl_*_Aligned.sortedByCoord.out.bam | paste -sd, - > 00_metadata/group1_bams.txt
ls 03_star/bam/Treat_*_Aligned.sortedByCoord.out.bam | paste -sd, - > 00_metadata/group2_bams.txt

文件内容示例：

03_star/bam/Ctrl_1_Aligned.sortedByCoord.out.bam,03_star/bam/Ctrl_2_Aligned.sortedByCoord.out.bam,03_star/bam/Ctrl_3_Aligned.sortedByCoord.out.bam

2. 运行 rMATS

mkdir -p 04_rmats/output 04_rmats/tmp

python rmats.py   --b1 00_metadata/group1_bams.txt   --b2 00_metadata/group2_bams.txt   --gtf ref/genes.gtf   --od 04_rmats/output   --tmp 04_rmats/tmp   -t paired   --readLength 150   --nthread 12   --libType fr-unstranded

关键参数：

参数	含义
`--b1` / `--b2`	两个比较组的 BAM 列表
`--gtf`	注释文件，必须与比对参考版本一致
`-t paired`	paired-end 数据；单端数据用 `single`
`--readLength`	read 长度，必须按实际数据填写
`--libType`	文库链特异性，如 `fr-unstranded`、`fr-firststrand`、`fr-secondstrand`
`--od`	结果输出目录
`--tmp`	临时目录

3. rMATS 输出解读

常见结果文件：

文件	说明
`SE.MATS.JC.txt`	exon skipping 事件，仅 junction counts
`SE.MATS.JCEC.txt`	exon skipping 事件，junction + exon counts
`A5SS.MATS.JC.txt`	alternative 5' splice site
`A3SS.MATS.JC.txt`	alternative 3' splice site
`MXE.MATS.JC.txt`	mutually exclusive exons
`RI.MATS.JC.txt`	retained intron

重要字段：

字段	解释
`IncLevel1`	group1 每个样本的 inclusion level
`IncLevel2`	group2 每个样本的 inclusion level
`IncLevelDifference`	group1 - group2 的 PSI 差异
`PValue`	原始 P 值
`FDR`	多重检验校正后的显著性
`IJC_SAMPLE_1/2`	inclusion junction counts
`SJC_SAMPLE_1/2`	skipping junction counts

九、rMATS 事件过滤与汇总

import pandas as pd
from pathlib import Path

event_types = ["SE", "A5SS", "A3SS", "MXE", "RI"]
out_dir = Path("04_rmats/output")
summary = []

for event in event_types:
    path = out_dir / f"{event}.MATS.JC.txt"
    df = pd.read_csv(path, sep="	")
    df["event_type"] = event

    sig = df[
        (df["FDR"] < 0.05)
        & (df["IncLevelDifference"].abs() >= 0.1)
    ].copy()

    sig.to_csv(out_dir / f"{event}.significant.tsv", sep="	", index=False)

    summary.append({
        "event_type": event,
        "total_events": len(df),
        "significant_events": len(sig),
        "up_in_group1": (sig["IncLevelDifference"] > 0).sum(),
        "up_in_group2": (sig["IncLevelDifference"] < 0).sum(),
    })

pd.DataFrame(summary).to_csv(out_dir / "rmats_event_summary.tsv", sep="	", index=False)

推荐阈值：

FDR < 0.05
abs(IncLevelDifference) >= 0.1
mean junction count >= 10

IncLevelDifference 的方向要小心：

rMATS 中通常是 group1 - group2。
如果 --b1 是 Ctrl，--b2 是 Treat，则正值表示 Ctrl inclusion 更高。
报告中必须写清楚方向。

十、SUPPA2 分析路线

SUPPA2 适合从 GTF 生成事件，再基于 transcript TPM 计算 PSI 和 differential splicing。

1. Salmon 转录本定量

mkdir -p 05_salmon ref/salmon_index

salmon index   -t ref/transcripts.fa   -i ref/salmon_index   -p 12

for sample in Ctrl_1 Ctrl_2 Ctrl_3 Treat_1 Treat_2 Treat_3
do
  salmon quant     -i ref/salmon_index     -l A     -1 02_clean_data/${sample}_R1.clean.fq.gz     -2 02_clean_data/${sample}_R2.clean.fq.gz     --validateMappings     --gcBias     --seqBias     -p 12     -o 05_salmon/${sample}
done

2. 生成 transcript TPM 矩阵

import pandas as pd
from pathlib import Path

samples = ["Ctrl_1", "Ctrl_2", "Ctrl_3", "Treat_1", "Treat_2", "Treat_3"]
tpm_tables = []

for sample in samples:
    quant = pd.read_csv(Path("05_salmon") / sample / "quant.sf", sep="	")
    tpm = quant[["Name", "TPM"]].rename(columns={"TPM": sample})
    tpm_tables.append(tpm)

expr = tpm_tables[0]
for tpm in tpm_tables[1:]:
    expr = expr.merge(tpm, on="Name", how="outer")

expr = expr.fillna(0)
expr.to_csv("06_suppa2/transcript_tpm.tsv", sep="	", index=False)

SUPPA2 表达矩阵第一列通常是 transcript ID，后面是样本 TPM。

3. 生成剪接事件

mkdir -p 06_suppa2/events

python suppa.py generateEvents   -i ref/genes.gtf   -o 06_suppa2/events/annotation   -e SE SS MX RI FL   -f ioe

常见事件：

SUPPA2 类型	含义
`SE`	skipped exon
`SS`	alternative splice site
`MX`	mutually exclusive exon
`RI`	retained intron
`FL`	alternative first/last exon

4. 计算 PSI

mkdir -p 06_suppa2/psi

for event_file in 06_suppa2/events/*.ioe
do
  base=$(basename ${event_file} .ioe)

  python suppa.py psiPerEvent     -i ${event_file}     -e 06_suppa2/transcript_tpm.tsv     -o 06_suppa2/psi/${base}
done

5. 差异剪接分析

准备分组表达矩阵：

import pandas as pd

expr = pd.read_csv("06_suppa2/transcript_tpm.tsv", sep="	")

ctrl_cols = ["Name", "Ctrl_1", "Ctrl_2", "Ctrl_3"]
treat_cols = ["Name", "Treat_1", "Treat_2", "Treat_3"]

expr[ctrl_cols].to_csv("06_suppa2/ctrl_tpm.tsv", sep="	", index=False)
expr[treat_cols].to_csv("06_suppa2/treat_tpm.tsv", sep="	", index=False)

运行 diffSplice：

mkdir -p 06_suppa2/diff

python suppa.py diffSplice   -m empirical   -gc   -i 06_suppa2/events/annotation_SE_strict.ioe   -p 06_suppa2/psi/annotation_SE_strict.psi   -e 06_suppa2/ctrl_tpm.tsv 06_suppa2/treat_tpm.tsv   -o 06_suppa2/diff/SE_Ctrl_vs_Treat

十一、rMATS 与 SUPPA2 结果整合

整合思路：

层级	rMATS	SUPPA2
事件坐标	明确 genomic coordinate	基于 transcript event
显著性	FDR	p-value / empirical model
效应量	IncLevelDifference	dPSI
可视化	Sashimi plot	PSI distribution

优先级建议：

rMATS 中 FDR < 0.05 且 abs(IncLevelDifference) >= 0.1 的事件。
SUPPA2 中 abs(dPSI) >= 0.1 且显著的事件。
两个工具方向一致的事件作为高置信候选。
与差异表达、RBP motif、功能结构域相关的事件优先解释。

十二、可视化

1. dPSI 火山图

library(tidyverse)
library(ggrepel)

se <- read.delim("04_rmats/output/SE.MATS.JC.txt")

plot_df <- se |>
  mutate(
    neg_log10_fdr = -log10(ifelse(FDR == 0, 1e-300, FDR)),
    regulation = case_when(
      FDR < 0.05 & IncLevelDifference >= 0.1 ~ "Higher inclusion in group1",
      FDR < 0.05 & IncLevelDifference <= -0.1 ~ "Higher inclusion in group2",
      TRUE ~ "Not significant"
    )
  )

ggplot(plot_df, aes(IncLevelDifference, neg_log10_fdr, color = regulation)) +
  geom_point(alpha = 0.7, size = 1.4) +
  geom_vline(xintercept = c(-0.1, 0.1), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  theme_bw() +
  labs(x = "ΔPSI / IncLevelDifference", y = "-log10(FDR)")

2. PSI 箱线图

library(tidyverse)

event_id <- 1
inc1 <- strsplit(se$IncLevel1[event_id], ",")[[1]] |> as.numeric()
inc2 <- strsplit(se$IncLevel2[event_id], ",")[[1]] |> as.numeric()

psi_df <- tibble(
  group = c(rep("Ctrl", length(inc1)), rep("Treat", length(inc2))),
  psi = c(inc1, inc2)
)

ggplot(psi_df, aes(group, psi, fill = group)) +
  geom_boxplot(width = 0.5, outlier.shape = NA) +
  geom_jitter(width = 0.08, size = 2) +
  theme_bw() +
  labs(y = "PSI", x = NULL)

3. Sashimi plot

推荐使用 IGV、rmats2sashimiplot 或 ggsashimi 对候选事件画 junction 支持图。

rmats2sashimiplot   --b1 00_metadata/group1_bams.txt   --b2 00_metadata/group2_bams.txt   -t SE   -e 04_rmats/output/SE.MATS.JC.txt   --l1 Ctrl   --l2 Treat   --exon_s 1   --intron_s 5   -o 07_visualization/sashimi_SE

十三、结果解释原则

指标	推荐解释
`FDR`	差异剪接显著性，常用 `< 0.05`
`IncLevelDifference`	PSI 差异，常用 `abs >= 0.1`
junction counts	支撑事件的 reads 数，太低不可靠
事件类型	SE、RI、A5SS 等类型的生物学意义不同
基因表达量	低表达基因的剪接事件要谨慎

注意：

可变剪接事件不一定伴随基因总表达变化。
需要区分“gene-level DEG”和“isoform/splicing-level change”。
报告中必须写清楚比较方向，例如 Ctrl - Treat 还是 Treat - Ctrl。
候选事件最好用 Sashimi plot 或 RT-PCR 进一步验证。

十四、常见问题与排查

问题	可能原因	处理建议
显著事件很少	reads 深度不足、重复少、剪接变化弱	检查 junction reads、降低解释预期
RI 事件异常多	intronic reads、pre-mRNA、rRNA 或比对问题	检查文库类型和 intronic mapping
两组方向解释反了	`--b1` / `--b2` 顺序混淆	在报告中固定方向并记录
SUPPA2 与 rMATS 不一致	事件定义和输入层级不同	优先看高 read support 事件
Sashimi plot 不支持结果	BAM 未 index 或坐标版本不一致	检查 BAM bai、GTF 和参考版本
FDR 显著但 PSI 很小	大样本下微弱变化显著	加入 `abs(dPSI)` 阈值

十五、主要交付物

FASTQ/MultiQC 质控报告
STAR sorted BAM 和 BAM index
STAR junction/mapping summary
rMATS 五类事件结果表
rMATS significant event tables
Salmon transcript TPM matrix
SUPPA2 event annotation、PSI、dPSI 结果
dPSI 火山图
PSI 箱线图
候选事件 Sashimi plots
高置信差异剪接事件清单
与 DEG/功能注释联动的解释报告

十六、参考资料

rMATS-turbo 官方仓库：提供 rMATS 构建、运行参数、输入输出和事件类型说明。
SUPPA/SUPPA2 官方仓库与教程：提供 generateEvents、psiPerEvent、diffSplice 等命令说明。
STAR 官方手册：用于 RNA-seq 剪接感知比对和 junction 支持。
Salmon 官方文档：用于 transcript-level quantification。