2019 7/28 コマンド修正
2019 8/3 リンク追加
一般にDNA配列解析はマッピングとそれに続く分析からなる（論文 図1）。マッピング段階では、BWA [論文より ref.1]のようなアライメントツールを介して既知の参照ゲノムにマッピングされる。その後、マッピングされたリードは、GATK [ref.2]のような変異検出のための解析ツールによって処理される。さまざまなアライメントツールと分析ツールが存在し、それぞれ独自のユースケースがある。
マッピングおよび解析ツールは、SAMファイル、またはその圧縮されたBAM / CRAM）[ref.4,5]の標準化されたファイルフォーマットを介してやり取りする。 SAMファイルは、クエリファイルのテンプレート名やセグメントシーケンス、アライメントツールによって生成された情報など、シーケンサによって生成されたリードに関する情報を格納するタブで区切られたファイルになっている。リファレンスゲノム、およびこれらのポジション[ref.6]にリードがどのくらい上手くマップされているかを記述するCIGAR文字列が含まれる。 SAMフォーマットは、オプションおよびツール固有の情報を格納できる非常に柔軟なフォーマットである。
実際には、さまざまなアライメントツールがわずかに異なる出力を生成し、異なる解析ツールは、わずかに異なるSAM構造に依存して動作する。例えば、分析ツールの中にはオプションの情報が必要な場合や、マップされていないリードを削除する場合や、リードが特定の順序で格納されている場合にのみ処理する場合など、リードをフィルタリングする必要があるものがある。このため実際には、アライメントツールと解析ツールの間に、SAMファイルを解析ツールで受け入れられるフォーマットに書き換えるためのステップが数多くある（論文 図2）。例えば、GATKベストプラクティス[ref.7]とbcbio-nextgenプロジェクト[ref.8]は、さまざまなアラインメントツールと分析ツールをうまく組み合わせるためにSAMToolsを呼び出す必要があると勧告している。
SAMtools [ref.3]とPicard（http://picard.sourceforge.net/）は、間違いなくSAMファイルを操作するために最も広く使われているツールである。これらは、リードのソートやフィルタリング、オプション情報の追加、マッピング位置に基づくポリメラーゼ連鎖反応（PCR）duplicationのマーキングなどのコマンドラインツールである。パイプラインスクリプトはこれらのコマンドのいくつかを次々と呼び出し、中間のSAMファイルを作成し、最終的に分析ツールへの入力として渡されるSAMファイルになる。
準備ステップに費やされる計算時間は無視できない。たとえば、Janssen Pharmaceuticaで使用されている5ステップの準備パイプラインを、全面的なBAMファイルで実行すると、標準の24コアサーバーで約1時間40分かかる。典型的な臨床実験では数百のBAMファイルを処理するため、準備ステップに費やされる計算時間は数百時間に達し、待機時間が長くなってしまう。クラウドの計算ノード上で必要なBAMファイルはexomeファイルの10〜20倍の大きさであるため、全ゲノムデータの計算上の課題はさらに厳しくなる。このランタイムとコストを大幅に削減するために、ソフトウェアを再設計する機会があることを示している。
複数のコマンドラインツールを順番に呼び出して準備パイプラインを作成する標準的な方法は、パフォーマンスの観点から次の欠点がある。
- BAM / CRAMの圧縮/解凍を含む、各コマンドラインの呼び出しで新しいSAMファイルが生成されるため、ステップ間に繰り返しファイルのI / Oが発生する。
- 各中間処理ツールは、特定の計算を実行するためにそのデータを表すSAMファイル全体を反復するため、同じインクリメンタルに変更されたデータに複数のトラバーサルが存在する。
- 並列化のチャンスは、各ツールによって同期ポイントが導入されるたびに制限される。
著者らは、どの準備段階から独立した準備パイプラインの実行でも、SAMファイルを1回だけ通過するソフトウェアアーキテクチャを提案する。このアーキテクチャは、elPrepという具象的なツールの形で実装されている。 elPrepを使用すると、ユーザーは適用するすべての準備ステップをリストする単一のコマンドのみを発行することになる。ソフトウェアは内部的に異なるステップの実行をマージして並列化する。これにより、中間ファイル用のストレージの命名と編成、Unixパイプなどの概念の理解、適用可能かどうかのエンドユーザーでの必要性が無くなり、ランタイムと関連するコストが削減される。
elprepの特徴
- efficient multi-threaded execution
- operates completely in-memory, no intermediate files are generated
- 100% equivalent output to output produced by SAMtools and Picard for overlapping functionality
- compatible with existing tools such as GATK, SAMtools, and Picard
- modular implementation
２の特徴の通りin memoryで動作するので、入力のsamファイルの6x程度のメモリが必要となる（メモリが不足すればシステムドライブにswapされ、結局I/Oが頻繁に起きてしまう）。ただしメモリ要求量はchromosomeごとなどでソートされているかでも変わってくる（See Github）。ほかの特徴として、elprepは拡張子でファイルを判断する（.sam/.bam/.cram）。入力がsamでも、出力を.bamにすれば bamファイルを出力することができる（samtoolsが必要）。unixのパイプにも対応しており、パースしてコマンドをつなぐことができる (この場合samとして渡す)。動作が高速なのは、複数の処理を同時に行えることも関係している。例えばduplication readsを除きながらmapping quality がゼロのリードを除く作業を一回のコマンドで実行できる。
上記の３だが、具体的には以下のツールの出力と互換性がある（記載のバージョンでチェックされている）。
- gatk 3.1-1, 3.7, 3.8
- gatk 2.7.2
- gatk 1.6
- samtools-1.0, samtools-1.1, samtools-1.2, samtools-1.5
- samtools-0.1.19
- picard-tools-2.9.2
- picard-tools-1.113
インストール
上記Githubには、テストデータとしてNA12878のエキソームシーケンスのchromosome 22のsamも含まれている。
#GOで導入できる。GO言語がなければbrewで導入しておく。
go get github.com/exascience/elprep #ダウンロード
go build github.com/exascience/elprep #ビルド elprepの実行ファイルがカレントにできる
#bioconda
conda install -c bioconda -y elprep
> ./elprep
$ elprep
elprep version 3.05 - see http://github.com/exascience/elprep for more information.
2018/03/31 09:56:10 Incorrect number of parameters.
Available commands: filter, split, merge
Filter parameters:
elprep filter sam-file sam-output-file
[--replace-reference-sequences sam-file]
[--filter-unmapped-reads]
[--filter-unmapped-reads-strict]
[--filter-mapping-quality mapping-quality]
[--filter-non-exact-mapping-reads]
[--filter-non-exact-mapping-reads-strict]
[--filter-non-overlapping-reads bed-file]
[--replace-read-group read-group-string]
[--mark-duplicates]
[--remove-duplicates]
[--remove-optional-fields [all | list]]
[--keep-optional-fields [none | list]]
[--sorting-order [keep | unknown | unsorted | queryname | coordinate]]
[--clean-sam]
[--nr-of-threads nr]
[--timed]
[--reference-t fai-file]
[--reference-T fasta-file]
Split parameters:
elprep split (sam-file | /path/to/input/) /path/to/output/
[--output-prefix name]
[--output-type [sam | bam | cram]]
[--single-end]
[--nr-of-threads nr]
[--reference-t fai-file]
[--reference-T fasta-file]
Merge parameters:
elprep merge /path/to/input sam-output-file
[--single-end]
[--nr-of-threads n]
[--reference-t fai-file]
[--reference-T fasta-file]
コマンド一覧にあるように、大きく分けてFilter、Split、Mergeのコマンドがある。Filterが一番内容が充実している。パスの通ったディレクトリに移動しておく。
ラン
filter
基本コマンドが以下になる。
elprep filter input.sam output.sam
入力にbamを使った時は内部でsamtoolsが動き（samtools view -h）、インメモリでsamに変換してから作業が行われる。
上の基本の型に、下に載せたオプションを追加することで様々なフィルタリングを同時に行うことができる。
- --replace-reference-sequences ヘッダの変更。希望のdictファイルを指定する。
- --filter-unmapped-reads アンマップリードを除く。
- --filter-unmapped-reads-strict アンマップリードをより徹底的に除く。
- --filter-mapping-quality mapping-quality マッピングクオリティが指定以下のリードを除く。
- --filter-non-exact-mapping-reads
- --filter-non-exact-mapping-reads-strict
- --filter-non-overlapping-reads bed-file
- --replace-read-group read-group-string
- --mark-duplicates
- --remove-duplicates
- --remove-optional-fields [all | list]
- --keep-optional-fields [none | list]
- --clean-sam
このコマンドは常に以下の順番で動作している。なので、例えば8のフラグと9のフラグを同時に立てても（9を先につけても、常に8 の処理が先に行われる）、正常にduplicatipn readsの除去は行うことができる。
- filter-unmapped-reads or filter-unmapped-reads-strict
- filter-mapping-quality
- filter-non-exact-mapping-reads or filter-non-exact-mapping-reads-strict
- filter-non-overlapping-reads
- clean-sam
- replace-reference-sequences
- replace-read-group
- mark-duplicates
- remove-duplicates
- remove-optional-fields
- keep-optional-fields
例えば以下の作業をワンライナーで行う。
- アンマップリードを完全に除く（--filter-unmapped-reads-strict）。
- duplicationのタグを付け（--mark-duplicates）duplicationリードを捨てる（--remove-duplicates）。
- mapping qualityがゼロのリード（bwaならmultiple mappingはゼロ）を捨てる（--filter-mapping-quality）
- samのcleaning（--clean-sam）。
elprep filter input.sam output.sam --mark-duplicates --remove-duplicates --filter-mapping-quality 0 --clean-sam
修正 "--filter-mapping-quality 0"ではmultiple mappingはフィルタリングされない。multiple mappingをフィルタリングするなら"--filter-mapping-quality 1"とする。
追記
すでにbamに変換したファイルも扱える。上と同じ作業を実行して、さらにcoordinateソートする。スレッドは12使う。
elprep filter input.bam output.bam --mark-duplicates --remove-duplicates --filter-mapping-quality 0 --clean-sam \
--nr-of-threads 12 --sorting-order coordinate
出力もbamになる。
mappingしてelprepに渡しフィルタリング、coordinate sortしたbam出力する。
#mapperはminimap2とする。
minimap2 -ax sr -R "@RG\tID:X\tLB:Y\tSM:sample1\tPL:ILLUMINA"\
-t 20 ref.fasta pair_1.fq.gz pair_2.fq.gz \
|elprep filter /dev/stdin filtered_sorted.bam \
--mark-duplicates --remove-duplicates \
--filter-mapping-quality 0 --clean-sam --nr-of-threads 20 \
--sorting-order coordinate \
&& samtools index filtered_sorted.bam
* パイプでの入力も対応しているが、入力は"-"でなく"/dev/stdin"とする。メモリが足りなければthread数を減らす。
split インプット (sam/bam/cram)をクロモソームごとに分割する。unmapリードは別出力する。
elprep split input.sam output/ --output-prefix "split-sam" --output-type sam --nr-of-threads 8 --single-end
ソートやindexファイルがなくても動作する。シングルエンドのinputは--single-endをつける。デフォルトでは入力をペアエンドとして扱う。
merge splitで分割した複数のインプット (sam/bam/cram)をマージする。
elprep merge input/ output.sam --nr-of-threads 8
シングルエンドのinputは--single-endをつける。デフォルトでは入力をペアエンドとして扱う。bamとcram出力にも対応している。
samやbamと同様にcramファイルも同様に扱えます。multi-socket serverでの使用や使い方についてのpythonスクリプトも用意されています。GithubのREADMEで確認してください。
https://github.com/ExaScience/elprep
追記
引用
elPrep: High-Performance Preparation of Sequence Alignment/Map Files for Variant Calling.
Herzeel C, Costanza P, Decap D, Fostier J, Reumers J.
PLoS One. 2015 Jul 16;10(7):e0132868.
