- 著者: Alexander Dobin, Carrie A Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, Thomas R Gingeras
- Corresponding author: Alexander Dobin (Cold Spring Harbor Laboratory, Cold Spring Harbor, NY, USA)
- 雑誌: Bioinformatics
- 発行年: 2013
- Epub日: 2012-10-25
- Article種別: Original Article
- PMID: 23104886
背景
RNA-Seq解析におけるリードとゲノムのアライメントは、ミスマッチやインデルの許容、およびスプライシングによる非連続的な転写産物構造を持つリードの高精度なマッピングという二重の課題を抱えている。既存のRNA-Seqアライナーの多くは、DNAアライナーにスプライス認識ロジックを後付けする方式を採用しており、例えば Trapnell et al. 2009 や GSNAP (Wu and Nacu 2010)、MapSplice (Wang et al. 2010)、RUM (Grant et al. 2011) などが挙げられる。このアプローチは、特に非正規スプライスやキメラ（融合）転写産物のマッピング精度に限界があり、マッピング速度も不足していた。さらに、これらのツールはリード長が200 bp以下であることを想定して設計されており、第三世代シーケンシング技術（PacBioやOxford Nanopore）によって生成される長リードには不適であった。ENCODE Transcriptomeのような大規模コホート研究では、480 billion readsという膨大なデータ量に対し、計算資源そのものがボトルネックとなり、効率的な解析が困難であった。これらの課題を克服する新しいアライメントツールの開発が不足しており、既存手法では大規模データセットの解析が手薄であるという課題が残されていた。特に、超高速かつ高精度なアライメントを同時に実現する汎用的なツールは未解明な点が多かった。
目的
本研究の目的は、ENCODE Transcriptomeのような膨大なRNA-Seqデータを効率的に処理可能とする新規RNA-SeqアライナーSTAR (Spliced Transcripts Alignment to a Reference) を開発することである。STARは、(1) 圧倒的な高速性、(2) 正規および非正規スプライス、キメラ（融合）転写産物の一括検出、(3) リード長による制約を持たない汎用設計（ショートリードからロングリード、第三世代シーケンスデータを含む）という要件を満たすことを目指した。
結果
マッピング速度とスケーラビリティ: STARは、ENCODE Transcriptome設定（ヒトゲノムに対し5億5000万の2 × 76 bpペアエンドリード）において、12コアサーバー上で1時間あたり5億4990万リードペアをマッピングし、全マッピングを完了した (Table 1)。これは、TopHat2と比較して約50倍、GSNAP (Genomic Short-read Nucleotide Alignment Program) やMapSpliceと比較しても10〜30倍以上の高速化に相当する。STARはスレッド数に対してほぼ線形のスループットスケーリングを示し、スレッド数を6から12に増やした場合でもスレッドあたりのマッピング速度の低下は約10%に留まった。例えば、6スレッドでは3億920万リードペア/時、12スレッドでは5億4990万リードペア/時であった (Table 1)。この速度向上は、大規模なRNA-seqデータセットの解析において極めて重要である。
シミュレーションデータにおけるマッピング精度: シミュレーションリード（既知スプライスサイト付与）に対するマッピングリコールと精度において、STARはTopHat2を上回るスプライスジャンクション検出精度を示した (Fig 2)。正規ジャンクションでは、リコール率96%以上、精度95%以上を達成した。最低検出閾値（ジャンクションあたり1リード）において、STARは高い感度を維持しつつ、最も低い偽陽性率を示した。例えば、偽陽性率1%の時点で、STARは他のアライナーと比較して高い真陽性率を示した (Fig 2)。この結果は、STARが低カバレッジ領域でも高精度なスプライス検出が可能であることを示唆する。
実験データにおけるマッピング精度: K562細胞株のENCODEロングRNA-seqデータセット（n=4000万の2 × 76リード）を用いた評価では、STARとGSNAPが最も高いリードのアライメント率（両者とも94%）を示した。スプライスジャンクション検出の精度指標（GENCODE 7 アノテーションに対する）では、STAR、RUM (RNA-seq Unified Mapper)、TopHat2が同様の性能を示した (Fig 3a-c)。特に、STARは最も低い偽陽性率（排他的に検出されたジャンクションの割合）を示し、同時にクラス内で2番目に高い偽感度を達成した (Fig 3e, f)。例えば、検出閾値5リードの時点で、STARの排他的に検出されたジャンクションの割合は約10%であり、他のアライナーよりも低い値を示した (Fig 3e)。
新規スプライスジャンクションの実験的検証: ENCODEプロジェクトの一環として、H1ES細胞 (n=1920 junctions) およびHUVEC細胞 (n=1920 junctions) から得られた1960個の予測新規遺伝子間スプライスジャンクションをRoche 454 RT-PCRシーケンスで実験的に検証した。その結果、80〜90%のジャンクションが成功裏に検証され（実在するbona fideスプライス）、STARのマッピング戦略の高い精度が実証された (Table 2)。特に、2リードのみでサポートされる候補ジャンクションであっても、H1ESで72.4%、HUVECで74.0%という高い検証率を維持した (Table 2)。これは、npIDRの予測値よりも高い検証率であり、STARの低カバレッジ領域における検出能力の高さを示唆する。
非正規スプライスおよびキメラ（融合）転写産物の検出: STARは、GT/AG以外のGC/AG、AT/ACなどの非正規モチーフを持つジャンクションも適切に検出可能であり、低発現または希少なスプライスイベントの発見に貢献する。また、複数の染色体や逆鎖をまたぐシードクラスターからキメラアライメントを生成し、BCR-ABLなどの癌融合遺伝子検出にネイティブに対応する。この機能は、癌ゲノム解析における融合遺伝子スクリーニングに有用である。
ロングリード対応: 第三世代シーケンシング技術によるロングリード（数kb）も、同一のフレームワーク内でマッピング可能である。シミュレートされたPacBioリード（平均1.5 kb）の95%以上を正確に配置できることを示した。これにより、リード長による設計上の上限を持たない「ユニバーサルアライナー」としての特性が実証された。
メモリ使用量: 高速性の代償として、ヒトゲノムのインデックスに約27 GBのRAMを必要とし、マッピング実行時にもサフィックスアレイ全体をメモリに常駐させる必要がある (Table 1)。これは、HISAT2の階層型FM-index（約8 GB）と比較して大きな要件である。しかし、STARにはスパースSAを使用するオプションがあり、その場合、ヒトゲノムで15.6 GBのRAM消費に抑えつつ、マッピング速度は約25%低下するものの、アライメント精度は維持される。
考察/結論
先行研究との違い: これまでのRNA-SeqアライナーがDNAアライナーの拡張として開発されたのと異なり、STARは非連続的なRNA-Seqリードを直接参照ゲノムにアライメントするためにゼロから設計された。これにより、既存のアライナーと比較してマッピング速度を最大450倍向上させつつ、アライメントの感度と精度も改善した。
新規性: 本研究で初めて、未圧縮サフィックスアレイを用いた逐次MMP検索とシードクラスタリング/スティッチングを組み合わせたアプローチが、超高速かつ高精度なRNA-Seqアライメントを実現することを新規に示した。また、正規および非正規スプライスジャンクション、キメラ（融合）転写産物のde novo検出、およびフルレングスRNAシーケンスのマッピングを単一のツールで可能にした点は、これまで報告されていない汎用性を提供する。
臨床応用: STARの高速性と高精度は、ENCODE Transcriptomeのような大規模コホート研究におけるRNA-Seqデータ解析を現実的な時間枠で完遂可能にし、ENCODE Phase III、GTEx、TCGA RNA-Seqパイプラインの標準アライナーとして採用された。特に、キメラ転写産物の検出能力は、癌における融合遺伝子検出など、臨床応用への大きな臨床的意義を持つ。
残された課題: STARの主なlimitationとして、ヒトゲノムのインデックスに約30 GBという比較的大きなRAM要件が挙げられる。これは、計算資源が制約される環境では選択肢を限定する可能性がある。また、超長リード（100 kb以上）の最適性においてはBLATやminimap2に劣る可能性があり、全長転写産物再構成にはStringTieやCufflinksなどの下流ツールとの組み合わせが必要である。今後の検討課題として、メモリ効率のさらなる改善や、より多様なシーケンシングプラットフォームへの最適化が挙げられる。
方法
STARは、主に2つの主要なステップから構成されるアライメントアルゴリズムを採用している。
Indexing: 参照ゲノムと既知のスプライスジャンクションGTF (Gene Transfer Format) ファイルから、非圧縮サフィックスアレイ (SA) を構築する。圧縮FM-indexとは異なり、SAは定数時間でのランダムアクセスを可能にし、検索の高速化に最適化されている。ヒトゲノムの場合、約30 GBのメモリを消費するが、これはBWAよりも大きい。
Step 1: Maximal Mappable Prefix (MMP) シード検索: 各リードの先頭から、SA上で正確に一致する最大マッパブルプレフィックスを順次検索する。ミスマッチ、ジャンクション、ポリAテール、アダプター、またはスプライス境界に遭遇するとMMPは終端し、その位置からリードの残りの部分に対して再度MMP検索を開始する「逐次シード拡張」を行う。これにより、スプライスされたリードでもエクソン部分ごとに独立してシードが取得され、本質的にスプライス認識型のアライメントが可能となる。このアプローチは、スプライスジャンクションの正確な位置をリード配列内で見つける自然な方法であり、リード配列を任意に分割するスプリットリード法よりも優位性がある。スプライスジャンクションは、スプライスジャンクションの遺伝子座や特性に関する事前の知識なしに、単一のアライメントパスで検出される。
Step 2: シードクラスタリングとスティッチング: ゲノムにアライメントされたシード群をゲノム位置でクラスタリングし、近接するシード同士を「スティッチング」して最終的なアライメントを構築する。スティッチングは、局所的なSmith-Watermanに相当するスコアリング（マッチ、ミスマッチ、ギャップ、ジャンクションペナルティ）を用いて最適経路を決定する。正規スプライスモチーフ（GT/AG、GC/AG、AT/AC）には低いペナルティが課せられ、非正規モチーフやキメラジャンクションには高いペナルティが課せられる。複数の染色体や逆鎖をまたぐシードクラスターは、キメラ（融合）アライメントとして報告される。ペアエンドRNA-seqリードのメイトからのシードは、同時にクラスタリングおよびスティッチングされ、各ペアエンドリードは単一のシーケンスとして扱われる。これにより、アルゴリズムの感度が向上し、一方のメイトからの1つの正しいアンカーだけでリード全体を正確にアライメントできる。
2-pass mode: 1パス目で検出された新規スプライスジャンクションを2パス目のインデックスに追加することで、低リードカバレッジ領域における新規ジャンクション検出感度を向上させる。
Output: SAM/BAM形式のアライメントファイル、スプライスジャンクションのタブ区切りファイル、キメラジャンクションのタブ区切りファイル、およびWiggleトラックを出力する。
実験条件: シミュレーションデータは、1000万の2 × 100nt Illuminaライクなリードをマウス転写産物から生成し、様々なゲノム変異とシーケンスエラーを導入した。実験データには、K562細胞株のENCODEロングRNA-seqデータセット（n=4000万の2 × 76リード）を使用した。新規スプライスジャンクションの検証には、H1ES (human embryonic stem cells) および HUVEC (human umbilical vein endothelial cells) 細胞株から得られたポリA+ロングRNA転写産物を用いた。統計解析は、npIDR (non-parametric irreproducible discovery rate) アプローチを用いて、複製間の再現性を評価した。データ比較にはStudent’s t-testおよびMann-Whitney U testが使用された。