- 著者: Anthony M. Bolger, Marc Lohse, Bjoern Usadel
- Corresponding author: Bjoern Usadel (RWTH Aachen / Forschungszentrum Jülich, Germany)
- 雑誌: Bioinformatics
- 発行年: 2014
- Epub日: 2014-04-01
- Article種別: Original Article
- PMID: 24695404
背景
次世代シーケンシング (NGS) 技術の急速な発展に伴い、生成されるシーケンスデータの量は爆発的に増加している。しかし、IlluminaプラットフォームをはじめとするNGSデータには、低品質な塩基やライブラリ調製時に混入するアダプタ配列などの技術的配列が含まれることが多く、これらが下流のバイオインフォマティクス解析の精度を著しく低下させることが知られている。特に、DNA断片がリード長よりも短い場合に発生する「アダプタリードスルー」現象では、リードの3’末端に部分的なアダプタ配列が残存し、これを正確に同定して除去することは容易ではない。
当時、Cutadapt (Martin 2011)、FASTX-Toolkit、Scythe/Sickle、AdapterRemoval (Lindgreen 2012)、EA-Utils (Aronesty 2013)、Reaperなどのリード前処理ツールが既に多数存在していた。しかし、これらの既存ツールは、柔軟な処理ステップの組み合わせ、ペアエンドデータの「ペアアウェアネス (ペア関係の維持)」の完全な保持、および高速な処理性能を同時に満たすには至っていなかった。特に、フォワードおよびリバースのFASTQファイル間の対応関係を維持しつつ、片方のリードのみがフィルタリングされた「シングルトンリード」を適切に分離し、かつアダプタリードスルーを高感度に検出する機構は極めて限定的であった。
複数の単一目的ツールをシェルパイプラインで組み合わせるアプローチは、Newick (Junier and Zdobnov 2010) のような他のドメインでは有効であったが、ペアアウェアネスの必要性から、リードペア間の接続が失われるという課題があった。この問題を解決するには、リードペアを再調整し、シングルトンリードを別途保存する追加ステップが必要となり、処理ステップがリードペアを単一のユニットとして評価できないという問題が残されていた。中間ファイルの生成を伴う一連のツール実行も、データサイズを考慮すると非効率的であった。このように、既存のツール群では、アダプタリードスルーの検出感度と特異性の両立、および低品質な長リードに対する適応的な品質フィルタリングの手法が未確立であり、バイオインフォマティクスパイプラインにおける大きな知識ギャップおよび技術的不足が存在していた。
目的
本研究の目的は、Illumina NGSリードデータの前処理において、以下の要件を同時に満たす新しいオープンソースツールTrimmomaticを開発することである。
- ユーザー定義の任意の順序で複数の処理ステップ（アダプタ除去、品質フィルタリング、スライディングウィンドウ等）を組み合わせ可能な高い柔軟性の実現。
- ペアエンドデータにおけるリード間の対応関係を自動的に保持し、フィルタリングされたシングルトンリードを適切に分離して保存するペアアウェア機能の実装。
- 大規模データセットに対応可能なマルチスレッド対応の高速かつメモリ効率の良いJava実装。
- アダプタリードスルーを含む技術配列を高感度に検出する新しいアルゴリズムの開発。
これらの機能を通じて、参照ベースアライメント (Li et al. Bioinformatics 2009) およびde novoアセンブリのいずれの下流解析においても、既存ツールを上回る出力品質を達成し、特に低品質データや長いリード長において顕著な改善をもたらすことを実証することを目的とした。
結果
Reference-based alignment におけるアライメント率の向上: Dataset 1 (HiSeq 2×100bpリード、n=11008190 replicates) を用いたBowtie 2によるアライメントにおいて、Trimmomaticによる前処理はリードのアライメント率を大幅に向上させた。特に、ミスマッチを許容しない厳密なアライメント設定において、未処理データのユニークアライメント数が6401927リードであったのに対し、Trimmomaticのadapters+MIモードで処理されたリードは8748401リードと、アライメントされたリード数が顕著に増加した (Table 1)。BWAを用いた検証でも同様の傾向が観察され、adapters+MIモードが厳密なアライメント設定で9056403リードと最高値を示した (Table 1)。これは、Trimmomaticの品質フィルタリング、特にMIモードが、アライメントの厳密性が高い場合に極めて有効であることを示している (Fig 3)。
低品質・長リードデータでの劇的な改善: Dataset 2 (MiSeq 2×250bpリード、n=801192 replicates) は、リバースリードの品質が約120塩基以降で大きく低下する低品質データセットであった。このような挑戦的なデータセットにおいて、Trimmomaticによる前処理は劇的な改善をもたらした。BWAの厳密モードでは、未処理データはわずか11592リード（約1.4%）しか完全にアライメントされなかったが、Trimmomatic adapters+MIで処理されたデータでは634779リード（約79.2%）がアライメントされ、アライメント率が50-fold以上向上した (Table 1)。寛容なアライメントモードでも、未処理の60010リードから639740リードへと10-fold以上の向上を実現した。アライメント率の向上は統計的に極めて有意であり (p<0.001)、TrimmomaticのMIモードが低品質データに対する高い堅牢性を持つことが実証された。
De novo assembly の品質改善: Trimmomaticによる前処理は、de novoアセンブリの品質も大幅に向上させた。Dataset 1では、contig N50サイズが未処理データの60370 bpから95389 bpへと1.58-fold増加し、最大contig長も約28%向上した。Dataset 2ではさらに顕著な改善が見られ、N50サイズが100662 bpから177880 bpへと1.77-fold増加し、最大contig長も55%増加した。また、Dataset 1の未処理データからアセンブルされたゲノムには34 bpのアダプタ配列の完全一致が混入していたが、Trimmomaticで処理されたデータからのアセンブリではアダプタの混入は一切検出されなかった。
他ツールとの性能および実行速度比較: Dataset 1を用いたBowtie 2の厳密なアライメント設定において、TrimmomaticのMIモードは8748401リードをアライメントし、他ツール（AdapterRemoval: 8103596リード、Scythe/Sickle: 8060612リード、Cutadapt: 8086428リード）を明確に上回る性能を示した (Table 2)。実行時間に関して、EA-Utils (EA-Utils) が最速（シリアル9.3秒/パラレル8.0秒）であったが、Trimmomaticもそれに僅差で続き（シリアル33.7秒/パラレル9.6秒）、AdapterRemoval（960.2秒）、Scythe/Sickle（529.3秒）、Cutadapt（342.5秒）などの既存ツールと比較して桁違いに高速であった (Table 2)。これらの実行時間は、n=3 replicatesの試行の中央値として算出された。
考察/結論
先行研究との違い: 既存の多くの前処理ツールは、ペアエンドデータの「ペアアウェアネス」を完全に保持しつつ、柔軟な処理ステップ構成と高速性能を両立できていなかった。Trimmomaticは、AdapterRemoval (Lindgreen 2012) が類似のアプローチを実装していたペアアウェアアダプタ除去を高速化しつつ、Maximum Informationという適応的トリミングパラダイムを新規に導入した点で、これまでの固定閾値ベースのツールとは対照的である。
新規性: 本研究で初めて、リード長、カバレッジ、エラー確率を組み合わせたMaximum Information品質フィルタリングアルゴリズムを開発し、その有効性を実証した。この適応的なトリミングアプローチは、リードのどの部分を保持すべきかを知覚的に判断し、特に低品質な長リードにおいて、既存の固定閾値ベースの品質フィルタリングよりも大幅な改善をもたらすことを新規に示した。また、わずか1塩基のアダプタ痕跡でも低い偽陽性率で検出可能なPalindrome modeも新規性のある機能である。
臨床応用: Trimmomaticの知見は、NGSデータを用いる広範な研究および臨床応用において重要な意義を持つ。がんゲノム解析における低頻度変異検出やRNA-seqによる遺伝子発現定量、de novoアセンブリなど、あらゆる下流のバイオインフォマティクスパイプラインにおいて、上流のデータ品質を改善することで、解析結果の信頼性と精度を向上させる。これにより、疾患診断、治療標的の同定、個別化医療の推進など、臨床現場でのゲノム解析の精度向上に直結する。
残された課題: 今後の検討課題として、Trimmomaticは現在のところIllumina以外のプラットフォーム（PacBio、Nanoporeなどのロングリードシーケンス技術）には最適化されていない点が挙げられる。また、バーコード由来のデマルチプレキシングは対象外である。さらに、Maximum Informationモードのパラメータs（厳密性）とtarget length tの最適値は、特定の応用やデータセットに応じて調整が必要であり、これらのパラメータの自動最適化メカズムの開発が今後の研究方向性として考えられる。
方法
TrimmomaticはJava 1.5以上で動作するクロスプラットフォームのツールとして、GPL (General Public License) v3ライセンスの下で実装された。本ツールは、アダプタ除去と品質フィルタリングにおいて、以下のアルゴリズムを導入した。
1. 技術配列（アダプタ）除去アルゴリズム
- Simple mode: 各リードを5’末端から3’末端に向かって独立して走査し、ユーザーが指定した技術配列との近似マッチを検出する。16塩基のシードを用いた「シードアンドエクステンド」アプローチを採用し、4ビット/塩基で符号化された配列に対してビットワイズXOR (exclusive OR) 演算とpopcount操作を行うことで、高速にミスマッチ数を計算する。アライメントスコアがユーザー定義の閾値を超えた場合、その領域以降を切除する。
- Palindrome mode: ペアエンドデータに特化し、DNA断片がリード長より短い「アダプタリードスルー」の検出に最適化されている。フォワードリードとリバースリードにそれぞれのアダプタ配列を前置し、両者をグローバルアライメントする。両リードの有効配列が逆相補関係にあり、かつ残りの部分がそれぞれの反対側のアダプタに一致するという条件を同時に検出することで、わずか1塩基のアダプタ痕跡でも極めて低い偽陽性率で検出可能である。
2. 品質フィルタリングアルゴリズム
- Sliding Window: 5’末端からリードを走査し、指定されたウィンドウサイズ内の平均品質スコアが閾値を下回った位置で3’側を切除する。
- Maximum Information (MI): リード長、カバレッジ、エラー確率の3つの要因を組み合わせた独自のスコアリング関数を用いて、最適なトリミング位置を決定する。ロジスティック曲線を用いて長さ閾値スコアをモデル化し、保持されるシーケンス長に基づく線形カバレッジスコア、および各塩基の正確である確率の積に基づくエラー率スコアを算出する。ユーザー定義の厳密性パラメータs（0から1）でカバレッジとエラー率のバランスを調整し、総合スコアを最大化する位置でトリミングを行う。
3. 検証プロトコル 本ツールの検証には、Escherichia coli K-12 (MG1655株) のゲノムDNAから調製されたNGSデータセットが用いられた。参照ベースアライメントにはBowtie 2 (Langmead et al. NatMethods 2012) およびBWA (Li et al. Bioinformatics 2009) が、de novoアセンブリにはVelvetが用いられた。比較対象ツールは、Cutadapt、FASTX (FASTX-Toolkit)、Reaper、EA-Utils、Scythe/Sickle、およびAdapterRemovalであった。評価には、公開されているE. coliデータセットであるSRX131047 (HiSeq 2000, 2×100bpリード) とSRR519926 (MiSeq, 2×250bpリード) が使用された。統計解析には、アライメント率の比較にカイ二乗検定が、N50値の比較にMann-Whitney U testが用いられた。