可变剪接分析¶

一句话概述¶

利用rMATS和SUPPA2从RNA-seq数据中检测和量化可变剪接事件（SE/A5SS/A3SS/MXE/RI），通过差异剪接分析发现条件特异性剪接变化，是转录后调控研究的核心分析手段。

核心知识点表格¶

知识点	关键内容	常用工具
可变剪接类型	SE/A5SS/A3SS/MXE/RI五种基本类型	-
PSI值	Percent Spliced In，剪接事件包含率	-
事件检测	从注释或比对中识别剪接事件	rMATS, SUPPA2, Leafcutter
差异剪接	条件间PSI差异的统计检验	rMATS, SUPPA2, DEXSeq
可视化	Sashimi plot展示junction reads	rmats2sashimiplot, ggsashimi
剪接调控	剪接因子motif、RBP binding	rMAPS, SpliceAid
长读长验证	全长转录本直接鉴定isoform	FLAIR, IsoQuant

各步骤详解¶

第一步：理解可变剪接类型¶

白话解释： 一个基因的pre-mRNA可以通过不同方式剪接产生多种成熟mRNA（转录本isoform）。就像一部电影可以有导演剪辑版和院线版。五种基本类型： - SE（Skipped Exon）：某个外显子被跳过 - A5SS（Alternative 5' Splice Site）：同一外显子使用不同的5'端 - A3SS（Alternative 3' Splice Site）：同一外显子使用不同的3'端 - MXE（Mutually Exclusive Exons）：两个外显子互斥，只能选一个 - RI（Retained Intron）：内含子未被剪掉，保留在mRNA中

技术细节： - PSI (Percent Spliced In) = 包含型reads / (包含型reads + 跳过型reads) - ΔPSI = PSI_condition1 - PSI_condition2，表示剪接变化程度 - 显著性标准：通常|ΔPSI| > 0.1 且 FDR < 0.05 - SE约占所有可变剪接事件的40%（人类）

代码示例：

# 基本概念图示
# SE (Skipped Exon):
#   Inclusion: =====[exon1]=====[exon2]=====[exon3]=====
#   Skipping:  =====[exon1]==================[exon3]=====

# A5SS (Alternative 5' Splice Site):
#   Long:  =====[exon_long]===========
#   Short: =====[exon_short]=========

# RI (Retained Intron):
#   Spliced:  =====[exon1]=====[exon2]=====
#   Retained: =====[exon1---intron---exon2]=====

第二步：RNA-seq数据比对¶

白话解释： 可变剪接分析依赖于junction reads（跨越剪接位点的reads）。需要使用支持剪接比对的工具（如STAR、HISAT2），它们能识别reads中的N-CIGAR操作（跨越内含子的gap）。

技术细节： - STAR比对速度快，剪接敏感性高，推荐用于rMATS - 需要两步比对（2-pass mode）提高新剪接位点发现率 - 比对质量直接影响PSI估计的准确性 - 建议使用最新的基因组注释（Ensembl/GENCODE）

代码示例：

# STAR两步比对
# 第一步：构建索引
STAR --runMode genomeGenerate \
    --genomeDir star_index \
    --genomeFastaFiles genome.fa \
    --sjdbGTFfile annotation.gtf \
    --sjdbOverhang 149 \
    --runThreadN 16

# 第二步：第一轮比对（发现新剪接位点）
for sample in ctrl1 ctrl2 ctrl3 treat1 treat2 treat3; do
    STAR --genomeDir star_index \
        --readFilesIn ${sample}_R1.fastq.gz ${sample}_R2.fastq.gz \
        --readFilesCommand zcat \
        --outSAMtype BAM SortedByCoordinate \
        --outSAMstrandField intronMotif \
        --outFilterMultimapNmax 1 \
        --alignSJoverhangMin 8 \
        --alignSJDBoverhangMin 3 \
        --outFileNamePrefix aligned/${sample}_ \
        --runThreadN 8
done

# 收集新发现的剪接位点
cat aligned/*_SJ.out.tab | awk '$7 >= 3' | cut -f1-6 | sort -u > novel_junctions.tab

# 第三步：第二轮比对（使用合并的剪接位点）
for sample in ctrl1 ctrl2 ctrl3 treat1 treat2 treat3; do
    STAR --genomeDir star_index \
        --readFilesIn ${sample}_R1.fastq.gz ${sample}_R2.fastq.gz \
        --readFilesCommand zcat \
        --outSAMtype BAM SortedByCoordinate \
        --sjdbFileChrStartEnd novel_junctions.tab \
        --outFilterMultimapNmax 1 \
        --outFileNamePrefix aligned_pass2/${sample}_ \
        --runThreadN 8

    # 索引BAM
    samtools index aligned_pass2/${sample}_Aligned.sortedByCoord.out.bam
done

第三步：rMATS差异剪接分析¶

白话解释： rMATS是最广泛使用的差异可变剪接检测工具。它统计每个剪接事件在两组样本中的PSI值差异，用似然比检验判断差异是否显著。支持两种计数模式：JC（只用junction reads）和JCEC（junction + exon body reads）。

技术细节： - 使用多元均匀分布(multivariate uniform)建模生物学变异 - 似然比检验(LRT)比较固定和自由模型 - JC模式更保守但更准确，JCEC模式sensitivity更高 - 支持成对(paired)设计 - 需要比对到同一参考注释

代码示例：

# 安装rMATS
# 注意：rMATS-turbo 是当前官方推荐版本（速度提升 100x），安装：conda install -c bioconda rmats-turbo
conda install -c bioconda rmats

# 准备BAM文件列表
echo "aligned_pass2/ctrl1_Aligned.sortedByCoord.out.bam,aligned_pass2/ctrl2_Aligned.sortedByCoord.out.bam,aligned_pass2/ctrl3_Aligned.sortedByCoord.out.bam" > b1.txt
echo "aligned_pass2/treat1_Aligned.sortedByCoord.out.bam,aligned_pass2/treat2_Aligned.sortedByCoord.out.bam,aligned_pass2/treat3_Aligned.sortedByCoord.out.bam" > b2.txt

# 运行rMATS
rmats.py \
    --b1 b1.txt \
    --b2 b2.txt \
    --gtf annotation.gtf \
    -t paired \           # -t paired 表示双端测序数据（paired-end reads），非配对实验设计
    --readLength 150 \
    --nthread 16 \
    --od rmats_output \
    --tmp rmats_tmp \
    --variable-read-length  # 允许可变read长度

# 输出文件（每种事件类型各一个）：
# SE.MATS.JC.txt         - Skipped Exon (仅junction counts)
# SE.MATS.JCEC.txt       - Skipped Exon (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

# 筛选显著差异剪接事件
awk -F'\t' 'NR==1 || ($20 < 0.05 && ($23 > 0.1 || $23 < -0.1))' \
    rmats_output/SE.MATS.JCEC.txt > significant_SE.txt

# 列说明：
# $20 = FDR
# $23 = IncLevelDifference (ΔPSI)
# $24 = IncLevel1 (组1的PSI值，逗号分隔)
# $25 = IncLevel2 (组2的PSI值，逗号分隔)

第四步：SUPPA2差异剪接分析¶

白话解释： SUPPA2是另一种主流差异剪接工具，特点是使用转录本定量结果（如Salmon/kallisto的TPM值）而非BAM文件，速度极快。它用经验分布比较PSI值差异的显著性。

技术细节： - 基于转录本定量（TPM/counts），不需要BAM - 支持两种模式：基于注释的事件检测(splicing events)和基于转录本的isoform分析 - 使用经验累积分布函数估计p-value - 支持更复杂的实验设计（时间序列、多因素） - 比rMATS快很多，适合大规模分析

代码示例：

# 安装SUPPA2
pip install SUPPA2  # PyPI 包名为 SUPPA2（大写），旧版 SUPPA 已废弃

# 步骤1: 从GTF生成事件
suppa.py generateEvents \
    -i annotation.gtf \
    -o events \
    -f ioe \
    -e SE SS MX RI FL

# 输出: events_SE_strict.ioe, events_A3_strict.ioe等

# 步骤2: 转录本定量（使用Salmon）
for sample in ctrl1 ctrl2 ctrl3 treat1 treat2 treat3; do
    salmon quant -i salmon_index \
        -l A \
        -1 ${sample}_R1.fastq.gz \
        -2 ${sample}_R2.fastq.gz \
        -o salmon_output/${sample} \
        -p 8
done

# 步骤3: 合并TPM表达矩阵
# 使用salmon的quant.sf文件合并
python3 << 'EOF'
import pandas as pd
import os

samples = ['ctrl1','ctrl2','ctrl3','treat1','treat2','treat3']
tpm_data = {}
for s in samples:
    df = pd.read_csv(f"salmon_output/{s}/quant.sf", sep="\t")
    tpm_data[s] = df.set_index("Name")["TPM"]

tpm_matrix = pd.DataFrame(tpm_data)
tpm_matrix.to_csv("transcript_tpm.txt", sep="\t")
EOF

# 步骤4: 计算PSI值
suppa.py psiPerEvent \
    -i events_SE_strict.ioe \
    -e transcript_tpm.txt \
    -o SE_psi

# 步骤5: 差异剪接检测
# 分割PSI和TPM矩阵为两组
head -1 SE_psi.psi > header.txt
# 准备ctrl和treat的PSI文件
paste <(cut -f1 SE_psi.psi) <(cut -f2,3,4 SE_psi.psi) > ctrl_psi.txt
paste <(cut -f1 SE_psi.psi) <(cut -f5,6,7 SE_psi.psi) > treat_psi.txt

# 运行差异分析
suppa.py diffSplice \
    -m empirical \
    -i events_SE_strict.ioe \
    -p treat_psi.txt ctrl_psi.txt \
    -e treat_tpm.txt ctrl_tpm.txt \
    -gc \                   # gene correction
    -o SE_diff

# 输出: SE_diff.dpsi (ΔPSI值和p-value)
# 筛选显著事件
awk '$2 != "nan" && ($2 > 0.1 || $2 < -0.1) && $3 < 0.05' \
    SE_diff.dpsi > significant_SE_suppa.txt

第五步：可视化（Sashimi Plot）¶

白话解释： Sashimi plot是展示可变剪接的标准图形。上面显示reads覆盖深度曲线，弧线表示junction reads（跨越外显子-外显子连接的reads），弧线上的数字表示支持该连接的reads数量。通过比较不同条件的sashimi plot可以直观看到剪接变化。

代码示例：

# 方法1: rmats2sashimiplot
conda install -c bioconda rmats2sashimiplot

rmats2sashimiplot \
    --b1 aligned_pass2/ctrl1.bam,aligned_pass2/ctrl2.bam,aligned_pass2/ctrl3.bam \
    --b2 aligned_pass2/treat1.bam,aligned_pass2/treat2.bam,aligned_pass2/treat3.bam \
    -t SE \
    -e rmats_output/SE.MATS.JCEC.txt \
    --l1 Control --l2 Treatment \
    -o sashimi_plots/ \
    --event-type SE \
    --exon_s 1 --intron_s 5

# 方法2: ggsashimi (基于ggplot2)
conda install -c bioconda ggsashimi

# 准备配置文件
echo -e "ctrl1\taligned_pass2/ctrl1.bam\tControl" > sashimi_config.tsv
echo -e "ctrl2\taligned_pass2/ctrl2.bam\tControl" >> sashimi_config.tsv
echo -e "treat1\taligned_pass2/treat1.bam\tTreatment" >> sashimi_config.tsv
echo -e "treat2\taligned_pass2/treat2.bam\tTreatment" >> sashimi_config.tsv

ggsashimi.py \
    -b sashimi_config.tsv \
    -c chr1:1000000-1050000 \
    -g annotation.gtf \
    -o sashimi_output.pdf \
    --aggr mean_j \
    --alpha 0.25 \
    -M 10

第六步：功能解读与剪接调控分析¶

白话解释： 发现差异剪接事件后，需要理解它的功能意义：该事件是否影响蛋白质结构域？是否受某个剪接因子调控？可以通过motif分析在剪接位点附近寻找RNA结合蛋白的结合位点。

代码示例：

# R: 差异剪接事件功能分析
library(tidyverse)
library(clusterProfiler)
library(org.Hs.eg.db)

# 读取rMATS结果
se_results <- read.table("rmats_output/SE.MATS.JCEC.txt", 
    header=TRUE, sep="\t")

# 筛选显著事件
sig_events <- se_results %>%
    filter(FDR < 0.05, abs(IncLevelDifference) > 0.1)

# 提取涉及的基因
sig_genes <- unique(sig_events$geneSymbol)

# GO富集分析
ego <- enrichGO(gene = sig_genes,
    OrgDb = org.Hs.eg.db,
    keyType = "SYMBOL",
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05)

# 蛋白结构域影响预测
# 检查跳过的外显子是否编码已知domain
library(biomaRt)
mart <- useMart("ensembl", dataset="hsapiens_gene_ensembl")

# 获取外显子对应的蛋白域信息
domain_info <- getBM(
    attributes=c("ensembl_exon_id","pfam","pfam_start","pfam_end"),
    filters="ensembl_exon_id",
    values=sig_events$exonID,
    mart=mart)

# 剪接调控motif分析
# 提取差异剪接外显子上下游序列，搜索RBP motif
# 使用rMAPS工具或MEME

实战命令¶

#!/bin/bash
# === 可变剪接分析完整流程 ===
set -e

GTF="annotation.gtf"
GENOME="genome.fa"
THREADS=16

# 1. STAR索引
STAR --runMode genomeGenerate --genomeDir star_idx \
    --genomeFastaFiles $GENOME --sjdbGTFfile $GTF --runThreadN $THREADS

# 2. 两步比对
for s in ctrl1 ctrl2 ctrl3 treat1 treat2 treat3; do
    STAR --genomeDir star_idx --readFilesIn ${s}_R1.fq.gz ${s}_R2.fq.gz \
        --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate \
        --outFileNamePrefix pass1/${s}_ --runThreadN $THREADS
done

cat pass1/*_SJ.out.tab | awk '$7>=3' | cut -f1-6 | sort -u > merged_SJ.tab

for s in ctrl1 ctrl2 ctrl3 treat1 treat2 treat3; do
    STAR --genomeDir star_idx --readFilesIn ${s}_R1.fq.gz ${s}_R2.fq.gz \
        --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate \
        --sjdbFileChrStartEnd merged_SJ.tab \
        --outFileNamePrefix pass2/${s}_ --runThreadN $THREADS
    samtools index pass2/${s}_Aligned.sortedByCoord.out.bam
done

# 3. rMATS分析
echo "pass2/ctrl1_Aligned.sortedByCoord.out.bam,pass2/ctrl2_Aligned.sortedByCoord.out.bam,pass2/ctrl3_Aligned.sortedByCoord.out.bam" > b1.txt
echo "pass2/treat1_Aligned.sortedByCoord.out.bam,pass2/treat2_Aligned.sortedByCoord.out.bam,pass2/treat3_Aligned.sortedByCoord.out.bam" > b2.txt

rmats.py --b1 b1.txt --b2 b2.txt --gtf $GTF -t paired \
    --readLength 150 --nthread $THREADS --od rmats_out --tmp rmats_tmp

# 4. 统计结果
for type in SE A5SS A3SS MXE RI; do
    sig=$(awk -F'\t' 'NR>1 && $20<0.05 && ($23>0.1 || $23<-0.1)' rmats_out/${type}.MATS.JCEC.txt | wc -l)
    echo "${type}: ${sig} significant events"
done

面试常问点¶

Q1: rMATS和SUPPA2的选择依据？ A: rMATS基于BAM文件计数junction reads，统计模型严格（似然比检验），适合小样本。SUPPA2基于转录本定量TPM，速度快，适合大样本和复杂设计。二者互补，重要发现建议用两种方法验证。

Q2: PSI值的计算原理？ A: PSI = 支持包含(inclusion)的reads数 / (支持包含 + 支持跳过(exclusion)的reads数)。对SE事件：inclusion reads = 跨上游外显子→目标外显子 + 目标外显子→下游外显子的junction reads；exclusion reads = 跨上游→下游外显子的junction reads。

Q3: 差异表达和差异剪接是一回事吗？ A: 不是。差异表达(DE)是基因总体表达量变化；差异剪接(DS)是同一基因不同isoform比例变化（总量可能不变）。一个基因可以无DE但有DS，或有DE无DS。应该分开分析。

Q4: 如何验证差异剪接事件？ A: 1）RT-PCR设计跨越可变外显子的引物验证；2）长读长测序（PacBio/Nanopore）直接观察全长isoform；3）Western blot检测蛋白isoform变化；4）Sashimi plot可视化。

Q5: 可变剪接分析需要多少测序深度？ A: 建议每个样本至少80-100M PE reads。junction reads只占总reads的一小部分（~5-10%），低深度时很多事件因reads不足无法可靠量化。

易错点¶

比对工具选错：不能用BWA等不支持剪接比对的工具，必须用STAR/HISAT2等splice-aware aligner。
Read length参数错误：rMATS的--readLength必须与实际read长度一致，否则junction计数会出错。
忽略生物学重复：至少需要3个生物学重复才能区分生物学变异和真正差异。无重复只能做描述性分析。
ΔPSI阈值过低：|ΔPSI| < 0.05的变化通常无生物学意义。建议阈值0.1-0.2。
混淆JC和JCEC模式：JC只用junction reads（更保守），JCEC加入exon body reads（更sensitive但可能引入bias）。报告时需说明使用哪种。
忽略低覆盖事件过滤：reads太少的事件PSI估计不可靠，应过滤（如IJC+SJC >= 10）。

补充知识¶