- 著者: Bo Li, Colin N Dewey
- Corresponding author: Colin N Dewey (Department of Computer Sciences, University of Wisconsin-Madison, Madison, WI, USA)
- 雑誌: BMC Bioinformatics
- 発行年: 2011
- Epub日: 2011-08-04
- Article種別: Original Article
- PMID: 21816040
背景
RNA-Seqは、マイクロアレイに代わる転写産物定量手法として急速に普及しているが、その転写量推定における主要課題は、リードが複数の遺伝子やアイソフォームに多重にマッピングされる点にある。特に近縁アイソフォーム間では配列が共通する領域が多いため、リードをユニークに割り付けることは困難であり、ユニークマッピングリードのみを用いる定量はバイアスを生むことが知られている。例えば、以前の研究では、リードマッピングの不確実性を統計モデルで適切に考慮することが、最も正確な発現量推定を達成するために重要であることが示されている (Li et al. Bioinformatics 2010)。
さらに、参照ゲノムが存在しない非モデル生物では、de novoトランスクリプトームアセンブリの結果を直接定量に利用する必要があるが、アイソフォーム同士が同一遺伝子由来かを判別できないため、既存の定量手法が前提とする「参照ゲノムへのアライメント」が成立しないという課題があった。多くの先行研究では、参照ゲノムへのアライメントを前提としていたため、このような状況での定量は困難であった。例えば、CufflinksやrQuantなどのツールはゲノムアライメントを必要とし、スプライシングやポリアデニル化によるリードのアライメントの複雑さに対処する必要があった。また、既存のEMアルゴリズム (Expectation-Maximization algorithm) に基づいたアイソフォーム定量法も提案されていたが、当時のモデルはシングルエンド・固定長・クオリティスコアなしという制約があり、現代のRNA-Seq実験プロトコルに対応しきれていない点が未解明であった。
このような状況において、多重マッピングリードの統計的処理、参照ゲノム非依存性、および多様なリード仕様への対応が不足しており、より汎用性の高い定量ツールの開発が求められていた。特に、リードマッピングの不確実性を考慮しない手法は、データの一部を破棄し、バイアスのある推定値を生むことが指摘されており、アイソフォームレベルの推定には不正確さをもたらすことが知られている (Wang et al. 2010)。この知識ギャップを埋めることが、RNA-Seqデータからの正確な転写産物定量において重要な課題であった。既存のツールでは、多重マッピングリードの処理が不十分であったり、参照ゲノムの存在を前提としたりするため、非モデル生物や複雑なトランスクリプトームの解析において、信頼性の高い発現量推定が不足していた。このため、より包括的かつ高精度な定量ソフトウェアの開発が急務であった。
目的
本研究の目的は、de novoトランスクリプトームを含む任意の参照トランスクリプトセットを入力として、遺伝子およびアイソフォームの発現量を高精度で推定する、ユーザーフレンドリーなRNA-Seq定量ソフトウェアパッケージ RSEM (RNA-Seq by Expectation Maximization) を提供することである。具体的には、多重マッピングリードをEMアルゴリズムで適切に分配し、ペアエンドリード (PE)、可変長リード、クオリティスコア、およびフラグメント長分布をモデルに統合する。
さらに、95%信頼区間 (CI) と事後平均推定値 (PME: posterior mean estimate) を出力すること、および異なるリード仕様 (リード数、リード長、シングルエンド (SE) vs ペアエンド) が定量精度に与える影響を体系的に評価することを目的とする。これにより、RNA-Seq実験の費用対効果の高い設計に関する実践的な指針を提供し、参照ゲノムが利用できない種や、複雑なトランスクリプトームを持つサンプルにおいても、信頼性の高い発現量定量を実現することを目指す。RSEMは、既存のツールが抱える参照ゲノム依存性や多重マッピングリード処理の限界を克服し、より幅広いRNA-Seq解析に対応可能な汎用ツールとなることを目指した。
結果
シミュレーションデータにおける定量精度: マウス由来のシミュレートデータ (2000万 SE 35bpリード) を用いた評価では、RSEMはCufflinks、IsoEM、rQuant、およびRSEM v0.6と比較して、同等以上の精度を示した (Figure 3)。遺伝子レベルの発現量推定における中央値パーセント誤差 (MPE: median percent error) はRSEMで最小であり、RefSeqデータセットで3.1%であった。アイソフォームレベルでも、RSEMはCufflinksに勝るか同等の精度を示した。参照ゲノムを与えない設定 (トランスクリプトFASTAのみ) でも、RSEMは参照ゲノムに依存する他手法と遜色ない精度を維持した。特に、CufflinksはEnsemblセットで異常に高い発現量推定値を生成し、短いトランスクリプトの処理に問題がある可能性が示唆された。IsoEMはRSEMと類似の性能を示したが、SEデータにおいてはRSEMがわずかに高い精度を示した。
多重マッピングリードの活用による精度向上: ユニークマッピングリードのみを使用する手法と比較して、多重マッピングリードをEMアルゴリズムで分配したRSEMの推定値は、真値との相関が顕著に高かった。特に、アイソフォームリッチな遺伝子 (スプライスバリアント数が多い遺伝子) でその差が大きかった。例えば、RefSeqシミュレーションデータにおいて、多重マッピングリードを考慮しない場合と比較して、RSEMは遺伝子レベルのMPEを大幅に改善した (Table 1)。多重マッピングリードを考慮しない場合、遺伝子レベルMPEは10.1%であったのに対し、RSEMでは3.1%に改善した。これは、リードマッピングの不確実性を統計的に取り込むことが定量精度の鍵であることを定量的に裏付けるものである。
実データでの再現性: MAQCのHuman Brain Reference RNAデータセットを用いて、RSEMの推定値をqRT-PCR測定値と比較した。その結果、RSEMのlog FPKMはqRT-PCRと高い相関 (Pearson r ≈ 0.86) を示し、参照ゲノム必須のツールと同等の精度であった (Table 2)。HBRサンプルではPearson r=0.69、UHRサンプルではPearson r=0.74を示した。ただし、MAQCデータセットは、単一アイソフォーム遺伝子や比較的ユニークな配列を持つ遺伝子に偏りがあり、アイソフォーム定量や多重マッピングリード処理の優位性を明確に区別するには限界があった。MAQCデータセットの遺伝子セットは、平均アイソフォーム数が1.1と、全遺伝子の平均1.7と比較して有意に低かった (p < 10⁻¹¹⁵)。
リード仕様の最適化に関する知見: 固定予算下で遺伝子レベルの定量精度を最大化する条件は「多数の短いSEリード」であることが示された。例えば、同じシーケンススループット (総塩基数 1400 MB) で比較した場合、2000万 PE 35bpリードよりも4000万 SE 35bpリードの方が、遺伝子レベルのMPEが低かった (マウスRefSeqデータでそれぞれ3.0% vs 2.3%) (Table 4)。これは、PE化やリード長伸長よりも総リード数増加の効果が大きいことを示唆する。一方、遺伝子内のアイソフォーム相対頻度の推定にはPEリードが有利であり、特にスプライスバリアント数が多い遺伝子 (ヒトEnsemblアノテーションなど) で改善が顕著であった。これは、実験設計を発現解析目的に応じて切り替えるべきという実践的指針を提供する。
De novoトランスクリプトーム対応: Trinityアセンブリに直接RSEMを適用したマウス実験では、参照ゲノムベースのパイプラインと比較して遺伝子レベルの発現量相関がr ≈ 0.95と高く、参照ゲノムが未整備な種でも信頼性のある定量が可能であることを示した。この実験では、n=20 million readsのシミュレーションデータが用いられた。これは、非モデル生物におけるRNA-Seq解析にRSEMが特に有用であることを裏付ける。
信頼区間と可視化: RSEMは95% CIとPMEを併記することで、点推定の不確実性を定量化できる。低発現アイソフォームではCI幅が広がる傾向が可視化された。シミュレーションデータを用いた検証では、95% CIが真の発現量を包含する割合は遺伝子レベルで94.3%、アイソフォームレベルで93.1%と、合理的な精度を示した (Table 6)。また、確率重み付きBAM出力により、IGVなどで多重マッピングリードの重みを直接確認できる点が、下流のレビューを容易にする (Figure 2)。
クオリティスコアの定量精度への影響: クオリティスコアを考慮したモデルと考慮しないモデル (リード配列のみを使用するモデル) の間で、定量精度に大きな差は認められなかった (Table 5)。Illuminaのようなエラープロファイルを持つリードでは、クオリティスコアを考慮したモデルが定量精度を大きく改善しないことが示された。これは、効果的なシーケンスエラーモデルがリード配列のみから学習可能であることを示唆する。
考察/結論
RSEMは、RNA-Seqにおける多重マッピングリード処理という根本課題に対し、EMアルゴリズムと完全な確率モデルで応える定量ツールとして開発され、広範に採用された。本研究で初めて、RSEMが参照ゲノムを必要とせずde novoトランスクリプトームで動作するため、非モデル生物、がんサンプルの融合トランスクリプト解析、新規アイソフォーム発見後の定量などで特に有用であることを示した。
先行研究との違い: これまでの多くの定量ツールが参照ゲノムへのアライメントを前提としていたのと異なり、RSEMは参照トランスクリプトセットのみで動作する点が画期的である。また、既存のEMベースのツールと比較して、RSEMはPEリード、可変長リード、フラグメント長分布、およびクオリティスコアをより詳細にモデル化しており、95% CIとPMEの計算機能も提供する点で優位性を持つ。IsoEMはRSEMと類似の性能を示すが、RSEMはPEデータにおいてわずかに高い精度を示し、ポリ(A)テール処理の改善がその一因であると考えられる。
新規性: 本研究で初めて、RNA-Seq実験の設計において、遺伝子レベルの定量には多数の短いSEリードが、遺伝子内のアイソフォーム相対頻度推定にはPEリードが最適であるという実践的な指針を新規に提供した。例えば、同じシーケンススループット (1400 MB) で比較した場合、4000万 SE 35bpリードは2000万 PE 35bpリードよりも遺伝子レベルMPEが優れていた (2.3% vs 3.0%)。また、クオリティスコア情報が定量精度に与える影響を評価し、Illuminaのようなエラープロファイルを持つリードでは、クオリティスコアを考慮したモデルがリード配列のみを使用するモデルと比較して、定量精度を大きく改善しないことを示した点も新規な知見である。
臨床応用: RSEMは、参照ゲノム情報が不足している非モデル生物の研究や、がん研究における新規融合トランスクリプトの定量など、臨床応用につながる幅広い分野で活用可能である。特に、95% CIの提供は、発現量の不確実性を定量化し、差次的発現解析の解釈に臨床的意義のある情報をもたらす。これにより、低発現遺伝子や多重マッピングの多い領域における推定値の信頼性を評価することが可能となる。
残された課題: 今後の検討課題として、アライメントベースであるRSEMの計算コスト (Salmonやkallistoなどのpseudoalignment系より遅い) の改善が挙げられる。また、低発現遺伝子におけるアイソフォームレベル推定の不確実性が残存する点、およびクオリティスコアモデルのキャリブレーションの必要性も今後の研究で取り組むべき課題である。さらに、シーケンス特異的なリード開始位置のバイアスやトランスクリプト特異的なリード分布など、追加のバイアスをモデルに組み込むことも今後の方向性として挙げられる。
方法
RSEMは、主に2つのステップで動作する。まず、参照トランスクリプトセットの準備と前処理を行う。これは、FASTA (Fast Alignment) 形式のトランスクリプト配列ファイル、またはGTF (Gene Transfer Format) 形式の遺伝子アノテーションファイルとゲノム配列ファイル (FASTA形式) のいずれかを入力として受け入れる。参照ゲノムは必須ではないため、de novoアセンブラが出力したトランスクリプトFASTAファイル単独でも動作する。RSEMは、ポリ(A) mRNA解析のために、参照トランスクリプトにポリ(A)テール配列を付加し、より正確なリードのアライメントを可能にする。
次に、RNA-Seqリードの参照トランスクリプトへのアライメントと、発現量の推定を行う。デフォルトでは、Bowtieアライメントプログラム (Langmead et al. GenomeBiol 2009) を使用し、RNA-Seq定量に特化したパラメータでリードをアライメントする。アライナーは、リードのすべての有効なアライメントを報告するように設定され、リードの短いプレフィックス内のミスマッチのみを考慮することで、RSEMがより詳細なモデルに基づいてアライメントの尤度を決定できるようにする。ユーザーはSAM (Sequence Alignment/Map) 形式でアライメントを提供することも可能である。
統計モデルは、リードの起源トランスクリプト、トランスクリプト内開始位置、リード方向、クオリティスコアを潜在変数として含む有向グラフィカルモデル (Figure 4) に基づいている。フラグメント長分布 (PEの場合) と位置特異的エラーモデルを組み合わせて尤度を計算する。多重マッピングリードは、EMアルゴリズムのEステップで各候補トランスクリプトに確率的に分配され、Mステップで最尤推定値 (ML estimate) が収束するまで反復計算される。収束基準は、値が10⁻⁷以上のすべてのθiについて相対変化が10⁻³未満となることである。また、Gibbsサンプラーにより事後分布から95% CIとPMEも算出される。EMアルゴリズムの高速化のため、多数 (デフォルトで200以上) のアライメントを持つリードはフィルタリングされる。
RSEMは、遺伝子およびアイソフォーム単位での期待カウント (expected count)、TPM (transcripts per million) 相当のフラクション推定値、各リードがどのトランスクリプトに重み付きで割り付けされたかを示す確率重み付きBAM (Binary Alignment/Map) ファイル、およびIGV (Integrative Genomics Viewer) などで閲覧可能なリード深度ウィグルピースを生成する。期待カウントは、Robinson et al. Bioinformatics 2010 やDESeqなどの差次的発現解析ツールで使用できる。さらに、RSEMには、設定したリード数、リード長、エラーモデルで仮想RNA-Seqデータを生成できるシミュレーターも内蔵されている。
本研究では、マウス RefSeq および Ensembl アノテーションセットを用いて、2000万 SE/PE リード (35bpまたは70bp) のシミュレーションデータセットを生成した。また、MAQC (MicroArray Quality Control) のHuman Brain Reference RNAデータセット (SRA accession SRX018974など) を用いて、実データでの性能を評価した。統計解析にはPearson correlation coefficientを用いた。MAQCデータセットは、ヒト脳組織 (HBR) と混合組織 (UHR) の2種類のRNAサンプルから構成され、qRT-PCRによる「ゴールドスタンダード」な発現量データが利用可能であった。RSEMの実行時間とメモリ使用量も評価され、8コア2.93 GHz Linuxサーバー (32 GB RAM) を使用して、2000万フラグメントのマウスRefSeqデータセットで測定された。