背景
在真核生物中，一个基因可以通过可变剪接（Alternative Splicing, AS）产生多种 mRNA 亚型 – 即在 pre-mRNA 剪接过程中，通过不同的外显子组合方式，生成不同的 成熟 mRNA。这一机制在不增加基因数目的前提下，极大扩展了蛋白质组的多样性。
一个含有 N 个外显子的基因，理论上最多可产生 2^(N-1) 种剪接变体。实际中， 大多数基因产生 2-10 种亚型，但个别基因（如果蝇 Dscam）可产生数万种。
为什么需要计算分析
高通量转录组测序（RNA-seq、全长 cDNA 测序）会产生数千条转录本比对到参考基因组 上的结果。人工逐个检查每个基因的剪接变异是不现实的。计算流程需要：
- 将属于同一基因位点的转录本聚类
- 区分真正的剪接变体与比对假象
- 对每种剪接变异进行分类
- 对全基因组的 AS 事件进行定量和汇总
AS Code 的概念
AS Code 体系（Sammeth et al., 2008）提供了一种紧凑、无歧义的记号来描述任意 可变剪接事件。其核心思想是：
对于任意一对重叠的转录本，它们之间不同的剪接位点定义了 AS 事件。将这些差异 位点按位置编号，并标记为供体（^）或受体（-），即可得到唯一描述该事件结构的 编码。
这一编码使得 AS 事件可以被系统地分类为具有生物学意义的类别。
输入数据与预处理
所需输入文件
+-------------------+------------------+-----------------------------------+
| 文件 | 格式 | 内容 |
+-------------------+------------------+-----------------------------------+
| 基因组 | FASTA | 参考基因组序列 |
| 基因注释 | GTF | 基因模型及外显子坐标 |
| cDNA 比对 | GFF3 (cDNA_match)| 全长 cDNA 比对结果 |
+-------------------+------------------+-----------------------------------+
基因注释提供参考转录本结构；cDNA 比对提供实验观测到的转录本结构，可能揭示注释中 未收录的新剪接变体。
解析基因/cDNA 注释
GTF（基因）和 GFF3（cDNA）文件均被解析以提取每个转录本的外显子坐标：
GTF 行: chr1 ensembl exon 1000 1200 . + . gene_id "G1"; transcript_id "T1";
GFF3 行: chr1 est cDNA_match 1000 1200 . + . ID=T2;
提取结果:
转录本 T1 位于 chr1: 外显子 [1000, 1200], 正链
转录本 T2 位于 chr1: 比对区 [1000, 1200], 正链
每个转录本存储为有序的外显子坐标列表：
转录本 T1: [[1000, 1200, '+'], [3000, 3500, '+'], [5000, 5500, '+']]
外显子1 外显子2 外显子3
从外显子推导内含子
内含子是同一转录本上相邻外显子之间的间隔。对于按起始位置排序的外显子：
外显子: [1000, 1200] [3000, 3500] [5000, 5500]
exon1 exon2 exon3
| | |
内含子: | [1201, 2999] | [3501, 4999] |
| intron1 | intron2 |
规则: intron = [exon_i.end + 1, exon_(i+1).start - 1, strand]
短内含子（默认 < 9 bp）被过滤，相邻外显子合并为一个外显子。这可以处理比对假象：
合并前: [1000, 1200] [1205, 1300] [3000, 3500]
exon1 短内含子 exon2
合并后: [1000, 1300] [3000, 3500]
合并后 exon1 exon2
提取剪接位点序列
对每个内含子，从基因组中提取两端各 6 bp 的序列：
基因组: ...ATCGAT|GT .....intron..... AG|GCTAGC...
donor acceptor
(6bp) (6bp)
Donor 位点: genome[intron.start-1 : intron.start-1+6] = "GT...."
Acceptor 位点: genome[intron.end-6 : intron.end] = "....AG"
负链上，坐标位置不变，但序列取反向互补，且 donor/acceptor 标签互换：
负链基因组: ...ATCGAT|GT .....intron..... AG|GCTAGC...
(基因组 5'->3')
生物学读取: 3'<-AG....|.....intron.....|CT....AC<-5'
acceptor donor
规范剪接位点为 GT/AG（主要）和 GC/AG（次要）。
计算基因范围
对每个基因，基因组跨度为其所有转录本中最左端外显子起始到最右端外显子终止：
基因 G1 的转录本:
T1: [1000, 1200] [3000, 3500] [5000, 5500]
T2: [1000, 1200] [3000, 3500]
T3: [1100, 1200] [3000, 3500] [5000, 5800]
基因范围: [1000, 5800, '+']
^最前起始 ^最后终止
核心工作流
总览
+------------------+
| 基因组 (FASTA) |
+--------+---------+
|
v
+------------------+ +------------------+ +------------------+
| 基因 GTF | | cDNA GFF3 | | 解析 FASTA |
| (参考注释) | | (实验证据) | | -> {chr: seq} |
+--------+---------+ +--------+---------+ +--------+---------+
| | |
v v |
解析 GTF 解析 GFF3 |
提取外显子 提取比对区 |
推导内含子 推导内含子 |
计算基因范围 计算基因范围 |
提取剪接位点序列 提取剪接位点序列 |
| | |
+-----------+-----------+ |
| |
v |
+---------------------+ |
| 步骤1: 转录本聚类 |<-----------------------+
| 基因-cDNA 重叠检测 |
+----------+----------+
|
v
+---------------------+
| 步骤2: 同一转录本 |
| 检测 (图聚类) |
+----------+----------+
|
v
+---------------------+
| 步骤3: AS Code |
| 生成 |
+----------+----------+
|
v
+---------------------+
| 步骤4: AS 类型 |
| 分类 |
+----------+----------+
|
v
+---------------------+
| 输出结果 |
| + SVG 可视化 |
+---------------------+
步骤1 – 转录本聚类
目标: 将每条 cDNA 转录本分配到其可能来源的基因。
方法: 检查 cDNA 的基因组范围是否与某个基因的范围重叠。
基因 G1 范围: [1000, 5800]
cDNA T4 范围: [1050, 5500] --> 重叠 --> T4 归入 G1 簇
cDNA T5 范围: [7000, 9000] --> 不重叠 --> T5 为新基因
基因 G2 范围: [6000, 8000]
cDNA T5 范围: [7000, 9000] --> 重叠 --> T5 归入 G2 簇
对每个基因簇，收集：
- 属于该基因的所有参考基因转录本
- 与该基因范围重叠的所有 cDNA 转录本
额外检查：
- 若 cDNA 的链方向与基因链方向冲突，标记为方向错误
- 没有任何 cDNA 覆盖的基因标记为"未验证"
- 不与任何基因重叠的 cDNA 标记为"新基因"
步骤2 – 同一转录本检测
目标: 在每个基因簇内，识别哪些转录本实际上是同一剪接亚型（应合并）， 仅保留真正的可变剪接差异。
这是算法最复杂的步骤，使用基于图的聚类方法。
问题本质
一个基因簇可能包含来自基因注释和 cDNA 证据的许多转录本。其中一些代表同一 剪接亚型（例如两条 cDNA 恰好比对到相同的外显子-内含子结构），而另一些则 代表真正不同的亚型。
基因 G1 簇包含: T1(基因), T2(基因), T3(cDNA), T4(cDNA), T5(cDNA)
可能的情况:
T1 和 T3 是同一亚型（仅来源不同）
T2 和 T5 是同一亚型
T4 是独立亚型
需要自动检测这种关系。
两两比较算法
两条转录本被认为"属于同一亚型"，当且仅当它们的内含子结构足够相似。 比较过程如下：
步骤 A: 检查链方向一致性 若两条转录本位于相反链，不可能是同一亚型。
步骤 B: 检查剪接位点一致性 若两条转录本都有内含子，其剪接位点序列（供体/受体）必须兼容。若一方在 某位置使用规范位点，另一方使用非规范位点，则不是同一亚型。
步骤 C: 逐内含子重叠检查
从 5’ 端开始，逐对比较两条转录本的内含子：
转录本 A 内含子: I1a I2a I3a
转录本 B 内含子: I1b I2b
比较 I1a 与 I1b:
- 若两者重叠，且双向覆盖率均 >= 阈值（默认 0.9）
--> 匹配，继续比较 I2a 与 I2b
- 若单侧覆盖率 > (1 - 阈值)
--> 差异过大，不是同一亚型，返回 FALSE
- 若 I1a 在 I1b 之前结束，且 I1a 被包含在 B 的基因范围内
--> A 有 B 没有的额外内含子，返回 FALSE
步骤 D: 延伸检查
当一条转录本的内含子先于另一条用完时：
情况1: A 的内含子比 B 多
A 剩余的内含子不能被 B 的外显子范围完全包含
（否则 A 有 B 缺失的额外剪接 --> 不同亚型）
情况2: 内含子数量相等且全部匹配
--> 同一亚型
图构建与聚类
构建无向图：
- 每个节点 = 一条转录本
- 若两条转录本"属于同一亚型"（按上述标准），则连一条边
图的连通分量代表同一亚型的转录本家族：
基因 G1 的图:
T1 ---- T3 T2 ---- T5 T4
(连通) (连通) (孤立)
连通分量:
家族1: {T1, T3} -- T1 为基因转录本, T3 为 cDNA（已验证）
家族2: {T2, T5} -- T2 为基因转录本, T5 为 cDNA（已验证）
家族3: {T4} -- T4 仅有 cDNA 证据（未验证）
在每个家族中，选择总外显子长度最长的转录本作为代表。其余转录本视为冗余， 从可变剪接分析中移除。
合并后剩余的转录本构成该基因的真正可变剪接变体集合。
步骤3 – AS Code 生成
目标: 对每对可变剪接转录本，计算一个紧凑编码，精确描述其剪接差异。
算法
给定两条转录本 A 和 B（均有内含子，位于同一链）：
第1步: 收集两条转录本的所有内含子坐标
转录本 A 内含子: [1201, 2999, +] [3501, 4999, +]
转录本 B 内含子: [1201, 2499, +] [3501, 4999, +]
第2步: 按基因组重叠聚类内含子
相互重叠的内含子被归入同一簇。每个簇代表一个可能存在剪接变异的区域。
簇1: {[1201,2999], [1201,2499]} -- 重叠的内含子
簇2: {[3501,4999], [3501,4999]} -- 完全相同的内含子（无变异）
第3步: 过滤簇
仅保留簇范围被两条转录本范围同时包含的簇（重叠类型3或4）。这确保我们 仅分析两条转录本共同覆盖区域内的变异。
第4步: 提取差异剪接位点
在每个簇内，收集所有内含子边界位置（起始和终止）。仅在一条转录本中出现 的位置为"差异位点"。
簇1 的边界:
A: start=1201, end=2999
B: start=1201, end=2499
共享: 1201（两者均从此起始）
差异: 2999（仅 A）, 2499（仅 B）
第5步: 过滤微小差异
若两个差异位点之间的距离小于阈值（默认 3 bp），视为测序/比对噪声并移除。
第6步: 分配 AS Code
将剩余差异位点按基因组坐标排序（负链则逆序）。编号 1, 2, 3, … 并标记 每个位点为供体（^）或受体（-）：
排序后的差异位点: 2499, 2999
位置 2499: 这是 B 内含子的终止端 --> 受体 --> 1-
位置 2999: 这是 A 内含子的终止端 --> 受体 --> 2-
转录本 A: 包含位置 2999 --> 编码 "2-"
转录本 B: 包含位置 2499 --> 编码 "1-"
AS Code: 2-,1-
供体/受体的判定取决于该位置是内含子起始端（供体, ^）还是终止端（受体, -）， 以及链方向：
正链:
内含子起始位置 --> 供体 (^)
内含子终止位置 --> 受体 (-)
负链:
内含子起始位置 --> 受体 (-) [因为基因组5'端是生物学的3'端]
内含子终止位置 --> 供体 (^)
AS Code 格式
一对转录本的 AS Code 写为：
编码_A,编码_B
其中每个编码由以下部分拼接: 位置编号 + 类型符号
位置编号: 1, 2, 3, ...（在事件中的相对顺序）
类型符号: ^（供体）或 -（受体）
若某条转录本在该事件中无差异位点，其编码为 "0"。
示例：
1^,2^ -- A 使用供体位点1, B 使用供体位点2
1-,0 -- A 使用受体位点1, B 无差异位点（外显子跳跃）
1^2-,3^4- -- A 使用位点1(供体)和2(受体), B 使用位点3(供体)和4(受体)
步骤4 – AS 类型分类
AS Code 被映射为六种生物学类别之一：
+----------+-----------+--------------------------------------------------+
| 类型 | 编码 | 生物学含义 |
+----------+-----------+--------------------------------------------------+
| ExonS | 0,X^ | 外显子跳跃: 一条转录本跳跃了另一条包含的外显子 |
| | X^,0 | |
+----------+-----------+--------------------------------------------------+
| IntronR | 0,X- | 内含子保留: 一条转录本保留了另一条剪接掉的内含子 |
| | X-,0 | |
+----------+-----------+--------------------------------------------------+
| AltD | 1^,2^ | 可变供体: 两条转录本使用不同的 5' 剪接位点 |
| | 2^,1^ | (供体位点) |
+----------+-----------+--------------------------------------------------+
| AltA | 1-,2- | 可变受体: 两条转录本使用不同的 3' 剪接位点 |
| | 2-,1- | (受体位点) |
+----------+-----------+--------------------------------------------------+
| AltP | X^Y-, | 可变位置: 内含子相互重叠但两端均不同 |
| | W^Z- | |
+----------+-----------+--------------------------------------------------+
| Other | (其余) | 不符合以上类别的复杂事件 |
+----------+-----------+--------------------------------------------------+
分类规则详解：
若一方编码为 "0":
若另一方以 "^" 结尾: --> ExonS （^ 标记供体，表示外显子-内含子边界，
意味着存在额外外显子）
若另一方以 "-" 结尾: --> IntronR（- 标记受体，表示内含子-外显子边界，
意味着保留了内含子）
若编码为 "1^,2^" 或 "2^,1^": --> AltD
若编码为 "1-,2-" 或 "2-,1-": --> AltA
若双方编码均匹配模式 [1-4][+/-][1-4][+/-]: --> AltP
其余: --> Other
AS 事件类型 – 结构可视化
符号说明
====== 外显子（实心方块）
------ 内含子（连接线）
^ 供体剪接位点（内含子 5' 端, 外显子->内含子边界）
| 受体剪接位点（内含子 3' 端, 内含子->外显子边界）
--> 链方向
外显子跳跃 (ExonS)
一条转录本包含一个盒式外显子（cassette exon），另一条完全跳过该外显子。
基因: [===E1===]-----------[===E2===]-----------[===E3===]
^donor ^acceptor
转录本 A (跳跃): 转录本 B (包含):
[===E1===]---------------------------[===E3===]
| 跳过的 E2 |
[===E1===]-------[===E2===]--------[===E3===]
AS Code: 0, 1^2-
解读:
转录本 A: 无差异位点 --> 编码 "0"
转录本 B: E1-E2 边界的供体(^), E2-E3 边界的受体(-)
--> 位置1和2 --> 编码 "1^2-"
非零编码以 "^" 开头但包含 "-" --> ExonS
内含子保留 (IntronR)
一条转录本保留了另一条剪接掉的内含子。
基因: [===E1===]---[===E2===]---[===E3===]
转录本 A (已剪接): 转录本 B (保留内含子):
[===E1===]-------[===E3===] [===E1===]---[intr]---[===E3===]
| | on |
| 内含子被剪接 | 内含子 |
| | 保留 |
详细视图:
转录本 A:
[===E1===]---------------------------[===E3===]
内含子被剪接
转录本 B:
[===E1===]---[==intron==]---[===E3===]
内含子保留为 mRNA 的一部分
AS Code: 0, 1-2^
解读:
转录本 A: 无差异位点 --> 编码 "0"
转录本 B: 内含子起始端的受体(-), 内含子终止端的供体(^)
--> 位置1和2 --> 编码 "1-2^"
非零编码以 "-" 开头 --> IntronR
与 ExonS 的关键区别:
IntronR: "0" 方转录本缺少内含子边界
非零方有内含子边界（先受体后供体: -...^）
ExonS: "0" 方转录本缺少外显子边界
非零方有外显子边界（先供体后受体: ^...-）
可变供体位点 (AltD)
两条转录本共享同一受体位点，但使用不同的供体位点，导致 5’ 剪接位置不同。
基因: [===E1===]---[===E2===]---[===E3===]
转录本 A (近端供体): 转录本 B (远端供体):
[===E1===]---[==E2==]---[===E3===] [===E1===]------[E2]---[===E3===]
^D1 ^D2 ^D2
详细视图（放大变异区域）:
转录本 A: 转录本 B:
[===E1===] [===E1===]
| D1 | D2
| v | v
+--.....intron.....---+ +--..........intron..........---+
| ^donor_A | | ^donor_B |
| | | |
[==E2==] [===E3===] [E2] [===E3===]
^同一受体 ^同一受体
AS Code: 1^, 2^
解读:
两个不同的供体位点，按位置编号1和2
转录本 A 使用供体1（近端）--> 编码 "1^"
转录本 B 使用供体2（远端）--> 编码 "2^"
双方编码均为单个供体 --> AltD
可变受体位点 (AltA)
两条转录本共享同一供体位点，但使用不同的受体位点，导致 3’ 剪接位置不同。
基因: [===E1===]---[===E2===]---[===E3===]
转录本 A (近端受体): 转录本 B (远端受体):
[===E1===]---[==E2==]---[===E3===] [===E1===]---[E2]------[===E3===]
^同一供体 ^同一供体 ^
A1 A2
详细视图（放大变异区域）:
转录本 A: 转录本 B:
[===E3===] [===E3===]
^A1 受体_A ^A2 受体_B
| |
+--.....intron.....---+ +--..........intron..........---+
| ^同一供体 | | ^同一供体 |
| | | |
[===E1===] [==E2==] [===E1===] [E2]
AS Code: 1-, 2-
解读:
两个不同的受体位点，按位置编号1和2
转录本 A 使用受体1（近端）--> 编码 "1-"
转录本 B 使用受体2（远端）--> 编码 "2-"
双方编码均为单个受体 --> AltA
可变位置 (AltP)
两个内含子相互重叠，但供体和受体位点均不同。
基因: [===E1===]---[===E2===]---[===E3===]---[===E4===]
转录本 A: 转录本 B:
[===E1===]----[==E2==]---[===E3===] [===E1===]---[E2]------[E3]---[===E4===]
^D_A ^A_A ^D_B ^A_B
详细视图（放大变异区域）:
转录本 A: 转录本 B:
^D_A ^D_B
| |
+--.....intron_A.....---+ +--........intron_B..........---+
| | | |
^A_A ^A_B
(同一区域, 两端边界均不同)
重叠示意图:
转录本 A 内含子: |<------D_A..............A_A------>|
转录本 B 内含子: |<---D_B...............A_B------------>|
^ ^
供体不同 受体不同
AS Code: 1^2-, 3^4-
解读:
共4个差异位点，按位置排序:
1: A 的供体 (^) 2: A 的受体 (-)
3: B 的供体 (^) 4: B 的受体 (-)
转录本 A: 位点1和2 --> 编码 "1^2-"
转录本 B: 位点3和4 --> 编码 "3^4-"
双方编码各有2个位点，且 ^ 和 - 混合 --> AltP
全部 AS 类型汇总
类型 结构 AS Code 关键特征
------ -------------------------------- -------- -------------------
ExonS [E1]----------[E3] 0, 1^2- 一方为0,
[E1]---[E2]--[E3] 非零方含 ^-
IntronR [E1]----------[E3] 0, 1-2^ 一方为0,
[E1]--[intr]--[E3] 非零方含 -^
AltD [E1]-D1--[E2]---[E3] 1^, 2^ 两个供体
[E1]-D2----[E2]-[E3]
AltA [E1]---[E2]-A1-[E3] 1-, 2- 两个受体
[E1]---[E2]---A2-[E3]
AltP [E1]-D1..A1-[E2]-[E3] 1^2-, 两端均不同
[E1]-D2....A2--[E3] 3^4-
Other 复杂多位点变异 various >4个位点或
混合模式
坐标比较引擎
坐标比较是贯穿整个流程的基础操作。给定两个区间 [s1, e1] 和 [s2, e2]， 判定其重叠关系。
六种重叠类型
类型1: 无重叠 (B 在 A 之前) 类型6: 无重叠 (A 在 B 之前)
s2----e2 s1----e1 s1----e1 s2----e2
| | | | | | | |
类型2: B 与 A 左侧重叠 类型5: B 与 A 右侧重叠
s2----s1----e2----e1 s1----s2----e1----e2
| | | | | | | |
| 重叠区域 | | 重叠区域 |
类型3: B 包含 A 类型4: A 包含 B
s2----s1----------e1----e2 s1----s2----------e2----e1
| | | | | | | |
| | A 在 B 内 | | | B 在 A 内 |
重叠长度:
类型1,6: 0
类型2: e2 - s1 + 1
类型3: e1 - s1 + 1 (= A 的长度)
类型4: e2 - s2 + 1 (= B 的长度)
类型5: e1 - s2 + 1
在流程中的应用
- 聚类: 检查 cDNA 范围是否与基因范围重叠（类型 != 1,6 即重叠）
- 同一转录本检测: 逐内含子检查重叠，计算覆盖率（类型2-5，计算 overlap / intron_length）
- AS Code 生成: 检查簇范围是否被两条转录本范围同时包含（类型3或4表示簇在转录本内部）
- 基因范围计算: 找到所有外显子中最外侧的坐标
输出规格
文件列表
输出目录结构:
alternative.splice.result/
|
|-- transcript.cluster.list.txt 基因 -> cDNA 簇映射
|-- alternative.splice.list.txt 基因 -> AS 转录本列表
|-- proved.transcript.list.txt 基因转录本 -> 匹配的 cDNA
|-- unproved.transcript.list.txt 无基因转录本匹配的 cDNA
|-- unproved.gene.list.txt 无 cDNA 证据的基因
|-- novel.gene.list.txt 不匹配任何已知基因的 cDNA
|-- error.orient.list.txt 链方向冲突的转录本
|
|-- acceptor.list.txt 受体位点序列（可选）
|-- donor.list.txt 供体位点序列（可选）
|
|-- splice.ascode.list.txt AS Code 详情（类 GTF 格式）
|-- splice.ascode.stat.txt AS 类型统计（可选）
|
|-- gene.cluster.picture/ SVG 可视化（可选）
|-- GENEID.svg
AS Code 输出格式
类 GTF 格式，附带自定义属性：
chr1 Undefined as_event 1201 2999 . + . transcript_id T1,T2; gene_id G1; \
structure 2-,1-; splice_chain 2999-,2499-
字段说明:
- 第1列: 染色体/scaffold
- 第2列: “Undefined”（占位符）
- 第3列: “as_event”
- 第4-5列: AS 事件的基因组范围（起始, 终止）
- 第6列: “.”
- 第7列: 链方向
- 第8列: “.”
- 第9列: 属性
- transcript_id: 被比较的两条转录本
- gene_id: 所属基因
- structure: AS Code（编码_A,编码_B）
- splice_chain: 坐标级编码（位置_A,位置_B）
AS 类型统计格式
AS_Number ExonS 1523
AS_Number IntronR 847
AS_Number AltD 2105
AS_Number AltA 1893
AS_Number AltP 672
AS_Number Other 234
###########################
Gene_Number ExonS 892
Gene_Number IntronR 501
Gene_Number AltD 1203
Gene_Number AltA 1067
Gene_Number AltP 398
Gene_Number Other 156
###########################
Code_Number 0,1^2- 892
Code_Number 0,1-2^ 501
Code_Number 1^,2^ 1203
Code_Number 1-,2- 1067
Code_Number 1^2-,3^4- 398
...
- AS_Number: 各 AS 类型的事件总数
- Gene_Number: 各 AS 类型涉及的独立基因数
- Code_Number: 各 AS Code 模式的出现次数
完整算法伪代码
主流程
FUNCTION main(genome_fasta, gene_gtf, cdna_gff3, params):
genome = parse_fasta(genome_fasta) // {chr: sequence}
gene = parse_gtf(gene_gtf, feature='exon') // 提取外显子坐标
gene.introns = derive_introns(gene.exons, min_length=params.min_intron)
gene.ranges = compute_gene_ranges(gene.exons)
gene.splice_seqs = extract_splice_sites(gene.introns, genome)
cdna = parse_gff3(cdna_gff3, feature='cDNA_match')
cdna.introns = derive_introns(cdna.exons, min_length=params.min_intron)
cdna.ranges = compute_gene_ranges(cdna.exons)
cdna.splice_seqs = extract_splice_sites(cdna.introns, genome)
// 步骤1: 将 cDNA 聚类到基因
clusters = {}
FOR each geneID, gene_range IN gene.ranges:
FOR each cdnaID, cdna_range IN cdna.ranges:
IF overlap(gene_range, cdna_range) > 0:
clusters[geneID].add(cdnaID)
clusters[geneID].add(all_transcripts_of(geneID))
// 步骤2: 检测同一转录本
FOR each geneID IN clusters:
G = empty_graph()
FOR each pair (ti, tj) IN clusters[geneID]:
IF is_same_transcript(ti, tj, params.coverage):
G.add_edge(ti, tj)
families = connected_components(G)
FOR each family IN families:
representative = longest_transcript(family)
alternative_splice[geneID].add(representative)
// 步骤3: 生成 AS Code
FOR each geneID IN alternative_splice:
FOR each pair (ti, tj) IN alternative_splice[geneID]:
as_code = compute_as_code(ti, tj, params)
IF as_code is not empty:
store as_code for (ti, tj)
// 步骤4: 分类 AS 类型
FOR each as_code IN all_as_codes:
as_type = classify_as_code(as_code)
count[as_type] += 1
genes[as_type].add(geneID)
OUTPUT all results
同一转录本检测
FUNCTION is_same_transcript(transcript_a, transcript_b, coverage):
IF strand(a) != strand(b): RETURN FALSE
IF splice_sites_incompatible(a, b): RETURN FALSE
introns_a = sorted(a.introns)
introns_b = sorted(b.introns)
range_a = (a.first_exon.start, a.last_exon.end)
range_b = (b.first_exon.start, b.last_exon.end)
i = 0
WHILE i < min(len(introns_a), len(introns_b)):
overlap = compute_overlap(introns_a[i], introns_b[i])
IF overlap / length(introns_a[i]) < coverage: RETURN FALSE
IF overlap / length(introns_b[i]) < coverage: RETURN FALSE
i += 1
// 检查剩余内含子
IF len(introns_a) > len(introns_b):
FOR j FROM i TO len(introns_a):
IF intron_contained_in_range(introns_a[j], range_b): RETURN FALSE
ELSE IF len(introns_b) > len(introns_a):
FOR j FROM i TO len(introns_b):
IF intron_contained_in_range(introns_b[j], range_a): RETURN FALSE
RETURN TRUE
AS Code 生成
FUNCTION compute_as_code(transcript_a, transcript_b, params):
all_introns = introns_a + introns_b
clusters = cluster_overlapping_introns(all_introns)
code_a = ""
code_b = ""
FOR each cluster IN clusters:
IF cluster NOT internal_to_both(cluster, range_a, range_b): CONTINUE
boundaries = collect_all_boundaries(cluster)
diff_sites = boundaries 中仅出现在一条转录本中的位置
diff_sites = filter_small_differences(diff_sites, params.asvaryedge)
IF diff_sites 为空: CONTINUE
SORT diff_sites by position（负链则逆序）
site_number = 1
FOR each site IN diff_sites:
symbol = "^" IF site 是供体 ELSE "-"
IF site 属于 transcript_a:
code_a += str(site_number) + symbol
ELSE:
code_b += str(site_number) + symbol
site_number += 1
IF code_a 为空: code_a = "0"
IF code_b 为空: code_b = "0"
RETURN (code_a, code_b)
AS 类型分类
FUNCTION classify_as_code(code_a, code_b):
IF code_a == "0" OR code_b == "0":
non_zero = code_a IF code_b == "0" ELSE code_b
IF non_zero 以 "^" 结尾: RETURN "ExonS"
IF non_zero 以 "-" 结尾: RETURN "IntronR"
IF (code_a == "1^" AND code_b == "2^") OR
(code_a == "2^" AND code_b == "1^"):
RETURN "AltD"
IF (code_a == "1-" AND code_b == "2-") OR
(code_a == "2-" AND code_b == "1-"):
RETURN "AltA"
IF 双方编码均匹配模式 /^[1-4][\^-][1-4][\^-]$/:
RETURN "AltP"
RETURN "Other"
关键参数及其影响
+-------------------+---------+-------------------------------------------+
| 参数 | 默认值 | 影响 |
+-------------------+---------+-------------------------------------------+
| coverage | 0.9 | 两个内含子被视为"匹配"的最低重叠比例。 |
| | | 较低的值更宽松（更多转录本被合并为 |
| | | "同一亚型"）。 |
+-------------------+---------+-------------------------------------------+
| minintronlength | 9 | 短于此长度的内含子被移除，相邻外显子 |
| | | 合并。防止比对间隙产生假内含子。 |
+-------------------+---------+-------------------------------------------+
| asvaryedge | 3 | 差异剪接位点间距在此范围内则被过滤。 |
| | | 减少不精确比对引入的噪声。 |
+-------------------+---------+-------------------------------------------+
| canonical | FALSE | 若为 TRUE，仅考虑具有规范剪接位点 |
| | | (GT/AG 或 GC/AG) 的内含子。 |
| | | 过滤非规范剪接事件。 |
+-------------------+---------+-------------------------------------------+
| type | exon | GTF feature 类型。"exon" 分析完整转录本 |
| | | 结构，"CDS" 仅分析编码区。 |
+-------------------+---------+-------------------------------------------+
完整计算示例
输入
基因 G1 位于 chr1 正链，有3条转录本：
T1 (基因): [100, 200] [500, 700] [900, 1100]
外显子1 外显子2 外显子3
T2 (基因): [100, 200] [500, 700]
外显子1 外显子2
T3 (cDNA): [100, 200] [600, 700] [900, 1100]
外显子1 外显子2 外显子3
步骤1: 推导内含子
T1 内含子: [201, 499, +] [701, 899, +]
T2 内含子: [201, 499, +]
T3 内含子: [201, 599, +] [701, 899, +]
步骤2: 聚类
三条转录本均与基因 G1 范围 [100, 1100] 重叠，因此 cluster = {T1, T2, T3}。
步骤3: 同一转录本检测
比较 T1 与 T2:
T1 有2个内含子, T2 有1个内含子
内含子1匹配（均为 [201, 499]）
T1 的第二个内含子 [701, 899] 在 T2 的范围 [100, 700] 内
--> T1 有 T2 缺失的额外剪接 --> 不是同一亚型
比较 T1 与 T3:
内含子1: T1=[201,499], T3=[201,599]
overlap = 499-201+1 = 299
coverage_A = 299/299 = 1.0 >= 0.9
coverage_B = 299/399 = 0.75 < 0.9 --> 不是同一亚型
比较 T2 与 T3:
T2 有1个内含子, T3 有2个内含子
数量不同且第一个内含子不满足覆盖率
--> 不是同一亚型
图: T1 -- T2 -- T3（无边）
三条均为独立的剪接变体。
步骤4: T1 与 T3 的 AS Code
所有内含子: [201,499,+], [701,899,+], [201,599,+], [701,899,+]
按重叠聚类:
簇 [201, 599]: {[201,499], [201,599]} -- T1 内含子1 与 T3 内含子1 重叠
簇 [701, 899]: {[701,899], [701,899]} -- 完全相同，无变异
处理簇 [201, 599]:
边界: 201(T1,T3), 499(T1), 599(T3)
共享: 201
差异: 499（仅 T1）, 599（仅 T3）
过滤: |599-499| = 100 > 3, 保留
按位置排序: 499, 599
位置 499: T1 内含子的终止端 --> 受体 --> 1-
位置 599: T3 内含子的终止端 --> 受体 --> 2-
T1 编码: "1-"
T3 编码: "2-"
AS Code: 1-,2- --> AltA（可变受体）
这在生物学上是合理的: T1 和 T3 共享位置 201 处的同一供体位点，但 T1 的 内含子在 499 处终止，T3 的内含子在 599 处终止，即 T3 使用了更远端的受体位点。
T1: [===E1===]---[===E2===]---[===E3===]
^D ^A1
T3: [===E1===]------[==E2==]---[===E3===]
^D ^A2
同一供体 (^D), 不同受体 (A1 vs A2) --> AltA