背景
真核生物において、一つの遺伝子は選択的スプライシング（Alternative Splicing, AS）によって複数の mRNA アイソフォーム（isoform）を産生しうる。すなわち、pre-mRNA のスプライシング過程で 異なるエクソンの組み合わせをとることで、異なる成熟 mRNA が生成されるのである。この仕組みは、 遺伝子数を増やすことなくプロテオーム（proteome）の多様性を飛躍的に拡大する。
N 個のエクソンを持つ遺伝子は、理論上最大 2^(N-1) 種のスプライスバリアントを産生しうる。 実際には多くの遺伝子が 2〜10 種のアイソフォームを産生するが、ショウジョウバエの Dscam のように 数万種のスプライスバリアント（splice variant）を産生する例外的な遺伝子も存在する。
計算解析が必要な理由
ハイスループットトランスクリプトームシーケンシング（RNA-seq、全長 cDNA シーケンシング）は、 参照ゲノムにマッピングされた数千本の転写産物を生み出す。これらを人手で一つひとつ検査することは 現実的ではない。計算パイプラインに求められるのは以下の処理である。
- 同一遺伝子座に属する転写産物のクラスタリング
- 真のスプライスバリアントとアライメントアーティファクト（alignment artifact）の識別
- 各スプライスバリアントの分類
- 全ゲノム規模での 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 未満）はフィルタリングされ、隣接エクソン同士が 一つのエクソンに統合される。これによりアライメントアーティファクト（alignment artifact）に 対処できる。
統合前: [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)
ドナー部位: genome[intron.start-1 : intron.start-1+6] = "GT...."
アクセプター部位: genome[intron.end-6 : intron.end] = "....AG"
マイナス鎖では座標位置は変わらないが、配列は逆相補鎖をとり、ドナー/アクセプターの ラベルが入れ替わる。
マイナス鎖ゲノム: ...ATCGAT|GT .....intron..... AG|GCTAGC...
(ゲノム 5'->3')
生物学的な読み: 3'<-AG....|.....intron.....|CT....AC<-5'
acceptor donor
標準的スプライス部位（canonical splice sites）は 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, '+']
^最前端開始 ^最後端終了
コアワークフロー
概要
flowchart TD
fasta["ゲノム (FASTA)"]
parse["FASTA パース<br/>→ {chr: seq}"]
gtf["遺伝子 GTF<br/>(参照アノテーション)"]
gtf_parse["GTF パース<br/>エクソン抽出<br/>イントロン導出<br/>アノテーション抽出<br/>スプライス部位抽出"]
gff3["cDNA GFF3<br/>(実験的証拠)"]
gff3_parse["GFF3 パース<br/>アライメント領域抽出<br/>イントロン導出<br/>アノテーション抽出<br/>スプライス部位抽出"]
step1["ステップ1：転写産物クラスタリング<br/>遺伝子-cDNA 重複検出"]
step2["ステップ2：同一転写<br/>産物検出 (グラフクラスタリング)"]
step3["ステップ3：AS Code<br/>生成"]
step4["ステップ4：AS タイプ<br/>分類"]
output["結果出力<br/>+ SVG 可視化"]
fasta --> parse
gtf --> gtf_parse
gff3 --> gff3_parse
parse --> step1
gtf_parse --> step1
gff3_parse --> step1
step1 --> step2
step2 --> step3
step3 --> step4
step4 --> output
ステップ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 エビデンス由来の多数の 転写産物が含まれうる。そのうち一部は同一のスプライスアイソフォームを表し （例：2 本の cDNA がたまたま同じエクソン-イントロン構造にアライメントされた場合）、 別の一部は真に異なるアイソフォームを表す。
遺伝子 G1 クラスタに含まれるもの: T1(遺伝子), T2(遺伝子), T3(cDNA), T4(cDNA), T5(cDNA)
想定される状況:
T1 と T3 は同一アイソフォーム（由来が異なるだけ）
T2 と T5 は同一アイソフォーム
T4 は独立したアイソフォーム
この関係を自動検出する必要がある。
ペアワイズ比較アルゴリズム
二つの転写産物が「同一アイソフォームに属する」とみなされるのは、両者のイントロン構造が 十分に類似している場合、かつその場合に限る。比較手順は以下のとおりである。
ステップ A: 鎖方向の一貫性チェック 二つの転写産物が逆鎖にある場合、同一アイソフォームではありえない。
ステップ B: スプライス部位配列の利用可否チェック スクリプトは両転写産物ともスプライス部位配列レコードを持つか（両方とも spliceSeq を持つ、 または両方とも持たない）をチェックする。canonical と non-canonical の状態は比較しない— レコードの存在のみを判定する。
ステップ C: イントロン重複チェック
スクリプトはネストループを用いて一致するイントロンペアを探す：A の各イントロンについて、 B の全イントロンを走査し、十分なカバレッジで重複するものを一つ見つける。一致が見つかると、 残りのイントロンは延長チェックにより順次比較される。
転写産物 A のイントロン: I1a I2a I3a
転写産物 B のイントロン: I1b I2b
A の各 Iia について、B の全 Ijb を走査:
- Iia と Ijb が重複し、双方向のカバレッジがいずれも閾値（デフォルト 0.9）以上
--> 一致、残りのイントロンについて延長チェックに進む
- 片側カバレッジが (1 - 閾値) を超える
--> 差異が大きすぎる、同一アイソフォームではない、FALSE を返す
- Iia が Ijb より先に終了し、かつ Iia が B の遺伝子範囲に包含される
--> A には B にない余分なイントロンがある、FALSE を返す
ステップ D: 延長チェック
一方の転写産物のイントロンが他方より先に尽きた場合：
ケース1: A のイントロンが B より多い（B が先に尽きた）
A の残りのイントロンが B のエクソン範囲に完全包含されてはならない
（さもなくば A に B が欠く余分なスプライシングがある → 異なるアイソフォーム）
ケース2: B のイントロンが A より多い（A が先に尽きた）
B の残りのイントロンが A のエクソン範囲に完全包含されてはならない
（さもなくば B に A が欠く余分なスプライシングがある → 異なるアイソフォーム）
ケース3: イントロン数が等しくすべて一致
--> 同一アイソフォーム
グラフ構築とクラスタリング
無向グラフを構築する。
- 各ノード = 1 本の転写産物
- 二つの転写産物が「同一アイソフォームに属する」（上記基準による）場合、辺を張る
グラフの連結成分が同一アイソフォームの転写産物ファミリーを表す。
遺伝子 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）のみを 保持する：転写産物範囲がクラスタを包含する（タイプ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'端であるため]
イントロン終了位置 --> ドナー (^)
実装上の注意：マイナス鎖では、スクリプトはコード生成ループ内で鎖を明示的に判定するのではなく、
コード生成前にソート済みイントロン座標リスト（cooris/coorjs）を逆転させることで
これを実現する。ループは常に「偶数インデックス = ドナー、奇数インデックス = アクセプター」
という規則を用い、リストの逆転によりマイナス鎖の生物学的 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 は 6 種の生物学的カテゴリのいずれかにマッピングされる。
| タイプ | コード | 生物学的意味 |
|---|---|---|
| ExonS | 0,X^ / X^,0 | エクソンスキッピング (Exon Skipping)：一方の転写産物が、他方が含むエクソンをスキップする |
| IntronR | 0,X- / X-,0 | イントロン保持 (Intron Retention)：一方の転写産物が、他方がスプライス除去したイントロンを保持する |
| AltD | 1^,2^ / 2^,1^ | 選択的ドナー (Alternative Donor)：二つの転写産物が異なる 5' スプライス部位（ドナー部位）を使用する |
| AltA | 1-,2- / 2-,1- | 選択的アクセプター (Alternative Acceptor)：二つの転写産物が異なる 3' スプライス部位（アクセプター部位）を使用する |
| AltP | X^Y-, W^Z- | 選択的位置 (Alternative Position)：イントロンが相互に重複するが両端とも異なる |
| Other | (その他) | 上記のいずれにも該当しない複雑なイベント |
分類規則の詳細：
一方のコードが "0" の場合:
他方が "^" で終わる: --> ExonS
（^ はドナーを標識し、エクソン-イントロン境界を表す。
これは余分なエクソンの存在を意味する）
他方が "-" で終わる: --> IntronR
（- はアクセプターを標識し、イントロン-エクソン境界を表す。
これはイントロンが保持されたことを意味する）
コードが "1^,2^" または "2^,1^" の場合: --> AltD
コードが "1-,2-" または "2-,1-" の場合: --> AltA
両コードがともにパターン [1-4][\^-][1-4][\^-] にマッチする場合: --> AltP
その他: --> Other
AS イベントタイプ -- 構造の可視化
記号の説明
| 記号 | 意味 |
|---|---|
| ====== | エクソン（exon、実線ブロック） |
| ------ | イントロン（intron、連結線） |
| ^ | ドナースプライス部位（donor、イントロン 5' 端、エクソン→イントロン境界） |
| | | アクセプタースプライス部位（acceptor、イントロン 3' 端、イントロン→エクソン境界） |
| --> | 鎖方向（strand） |
エクソンスキッピング (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===]---[intr]---[===E3===] [===E1===]-------[===E3===]
| on |
| イントロン | イントロンが
| 保持 | スプライス除去
| |
詳細図:
転写産物 A:
[===E1===]---[==intron==]---[===E3===]
イントロンが mRNA の一部として保持
転写産物 B:
[===E1===]---------------------------[===E3===]
イントロンがスプライス除去
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] / [E1]---[E2]--[E3] | 0, 1-2^ | 一方が 0、非ゼロ側に -^ を含む |
| IntronR | [E1]----------[E3] / [E1]--[intr]--[E3] | 0, 1^2- | 一方が 0、非ゼロ側に ^- を含む |
| AltD | [E1]-D1--[E2]---[E3] / [E1]-D2----[E2]-[E3] | 1^, 2^ | 二つのドナー |
| AltA | [E1]---[E2]-A1-[E3] / [E1]---[E2]---A2-[E3] | 1-, 2- | 二つのアクセプター |
| AltP | [E1]-D1..A1-[E2]-[E3] / [E1]-D2....A2--[E3] | 1^2-, 3^4- | 両端とも異なる |
| Other | 複雑な多部位変異 | various | >4 部位または混合パターン |
座標比較エンジン
座標比較はワークフロー全体を貫く基礎的な操作である。二つの区間 [s1, e1] と [s2, e2] が与えられたとき、その重複関係を判定する。
6 種の重複タイプ
タイプ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
// スプライス部位チェック：両方ともレコードを持つか両方とも持たない（存在のみ判定、
// canonical/non-canonical の比較はしない）
IF splice_seq_presence_mismatch(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)
// ネストループ：A の各イントロンについて、B の全イントロンを走査して一致を探す
FOR i IN range(len(introns_a)):
FOR j IN range(len(introns_b)):
overlap = compute_overlap(introns_a[i], introns_b[j])
IF overlap / length(introns_a[i]) > coverage AND
overlap / length(introns_b[j]) > coverage:
// 一致、残りのイントロンについて延長チェック
RETURN extend_overlap(introns_a[i+1:], introns_b[j+1:],
range_a, range_b)
IF overlap / length(introns_a[i]) > 1-coverage OR
overlap / length(introns_b[j]) > 1-coverage:
RETURN FALSE // 差異が大きすぎる
IF introns_a[i] が introns_b[j] より先に終了 AND
introns_a[i] が range_b に包含される: RETURN FALSE
IF introns_b[j] が introns_a[i] より先に終了 AND
introns_b[j] が range_a に包含される: RETURN FALSE
RETURN FALSE
FUNCTION extend_overlap(remaining_a, remaining_b, range_a, range_b):
i = 0
WHILE i < min(len(remaining_a), len(remaining_b)):
overlap = compute_overlap(remaining_a[i], remaining_b[i])
IF overlap / length(remaining_a[i]) < coverage: RETURN FALSE
IF overlap / length(remaining_b[i]) < coverage: RETURN FALSE
i += 1
// 残りのイントロンを検査（相手の範囲に包含 => 異なる）
IF len(remaining_a) > len(remaining_b):
IF intron_contained_in_range(remaining_a[i], range_b): RETURN FALSE
ELSE IF len(remaining_b) > len(remaining_a):
IF intron_contained_in_range(remaining_b[i], 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: クラスタリング
3 本の転写産物はいずれも遺伝子 G1 の範囲 [100, 1100] と重複するため、 cluster = {T1, T2, T3} となる。
ステップ3: 同一転写産物検出
T1 と T2 の比較:
T1 は2つのイントロン, T2 は1つのイントロンを持つ
イントロン1は一致（ともに [201, 499]）
T1 の 2 番目のイントロン [701, 899] と T2 の範囲 [100, 700] の比較:
701 > 700、当該イントロンは T2 の範囲内にない（重複タイプ1）
--> T1 の余分なイントロンは共有領域の外にある --> 同一アイソフォーム
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つのイントロンを持つ
最初のイントロン: T2=[201,499], T3=[201,599]
coverage_T3 = 299/399 = 0.75、(1-0.9)=0.1 より大きいが 0.9 未満
--> 差異が大きすぎる、同一アイソフォームではない
グラフ: T1 ---- T2（同一）, T3（孤立）
ファミリー: {T1, T2}, {T3}
{T1, T2} の代表: T1（より長く、エクソン3つ）
選択的スプライス転写産物: T1, 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