elPrep 4 - macでインフォマティクス

2020, 2021 9/17 追記

elPrep 4はelPrep [ref.1]の大幅に拡張された再実装であり、DNAシーケンシングパイプラインでのバリアントコールのシーケンスアライメント/マップファイル（SAM / BAM）[ref.2]を準備するためのマルチスレッドツールである。パイプラインでどの準備ステップが使用されるかはアプリケーションによって異なるが、一般的に、統計解析のために何らかの方法でアラインメントされたリードデータを準備する。また、マッピングされていないリードまたは関心のあるゲノム領域に基づいてリードを除外するステップを含む。座標順のリードソート、optical duplicateまたはPCR duplicateのリードのマーク、basic quality scoreの再較正の計算および適用など。たとえば、GATKベストプラクティス[ref.3]では、最も広く使用されているバリアントコール元の1つであるGATK [ref.4]を使用して、バリアントコール用のデータを準備するための4ステップパイプラインと2つのバリエーションを定義している。

　elPrepは、SAMtools [ref.5]、Picard、GATK 4 [ref.4]などのSAM / BAMファイルを処理する他のツールとは異なり、データを1回だけパススルーすることでシーケンスパイプラインを実行できる。このソフトウェアアーキテクチャは、実行中にデータを可能な限りメモリ内に保持することによってファイルI/Oの繰り返しを回避し、異なる準備ステップの計算をマージし、実行を並列化しながら不要な同期を回避するように設計されている。

　elPrep 4は、Googleが開発したオープンソースのプログラミング言語GoでelPrep [ref.1]を完全に再設計および再実装したものである。 Goは、メモリの安全性、並列ガベージコレクション、型推論、およびマルチコアを利用した並行処理のサポートを特徴とする静的型付けされコンパイルされた言語であり、elPrepのパフォーマンスをさらに向上させる新しいソフトウェア最適化戦略にアクセスできる。 elPrepの最初の実装はCommon Lispで書かれていた。これはオプションの型宣言、コードインライン化、スタックベースのメモリ割り当て、そしてマルチスレッド機能のおかげで低レベルのパフォーマンス最適化をよくサポートする言語である。

　elPrepのようなシーケンシングアプリケーションに特有の1つの局面は、それが数百ギガバイトのデータを処理する必要があり、メモリ管理に多大な圧力をかけているということである[ref.6]。ほとんどのCommon Lispの実装は現在、stop-and-copy、stop-the-worldガベージコレクタを使用している。これは、elPrepのマルチスレッド実行がプログラムを頻繁に一時停止するため、あまりにも妨げになるため、オフにする必要があった。ガベージコレクションがなければ、メモリを再利用し、不要なメモリ割り当てをできるだけ回避し、プログラミングの複雑さを増し、elPrepを維持するという厳格なプログラミングスタイルを採用する必要があった。 Goには、この問題を解決する並行並列ガベージコレクタが付属している[ref.6]。 Goに切り替えることによるその他の利点には、移植性のある無料のコンパイラ、および型推論、デフォルトのUTF8、コンパイラによるエスケープ分析などの現代的な言語機能が含まれる。

　新しいelPrep 4フレームワークでは、新しい機能をより簡単に追加したり、GATKベストプラクティス[ref.3]で説明されている準備手順をすべて実装したりすることもできる。 2つの重要な貢献として、elPrepフレームワークでの効率的な並列実行のために最適化された、optical duplicatesのマーキングとbasic quality scoreの再較正のアルゴリズムがある。これらのアルゴリズムは、PicardとGATK 4のオリジナルのアルゴリズムと比較して、中間ファイルの使用を避け、データに対する複数の繰り返しループを避け、そして並列で、一方、PicardとGATK 4でのそれぞれの実装と比較して同じ結果が得られる。

　elPrepのGATKベストプラクティスの4ステップパイプラインをベンチマークし、それをGATK 3.8とGATK 4の両方のランタイムと比較することで、elPrep 4がシーケンスパイプラインの実行にかかる実行時間とリソースコストを劇的に削減することを示す。また、アマゾンウェブサービス（AWS）でスケーリング実験を行い、全エクソームと全ゲノムの両方のデータを処理するために、elPrep 4とGATK 4を実行した場合のコストを比較する。

elPrep 4では、以下の新機能が導入された。

BQSR：BQSRを実行するためのオプション（-bqsr）を追加した。
optical duplicatesマーキング：optical duplicatesマーキングを実行するためのオプション（-mark-optical-duplicates）を追加した。
メトリクス：elPrepは、マップされていないリード、セカンダリのリード、リードの重複、basic quality scoreなどの統計情報を含むメトリクスファイルを生成するようになった。 GATK 4. elPrepメトリクスファイルのフォーマットはPicard / GATK 4のものと同じで、視覚化のためにMultiQC [8]と互換性がある。
BAM解析：elPrep 4は、以前はBAM解析のためにSAMtoolsを呼び出す依存があったが、現在はGoの組み込みgzip圧縮ライブラリを使用してBAM解析自体を実装している。実行時の観点から圧縮がより効率的になった。
VCF解析：elPrep 4はVCF解析を提供する。
BEDファイルで指定されたゲノム領域に基づいてリードをフィルタリングする。これはSAMtools / Picard / GATKの-Lオプションに似たオプションである。
統合分割 - フィルタ - マージ（sfm）モード：elPrepには2つの実行モードがある。すなわち、完全にRAM内で動作するモードと、処理のためにゲノム領域を使用してデータを分割するモードである。これは以前はPython スクリプトを使って実装されていたが、これらはGoで実装されたsfmサブコマンドに置き換えられ、elPrepのインストールと使用がより簡単になった。
これらの新機能に加えて、実行時間とメモリ使用量の両方を減少させるさまざまなパフォーマンスの改善がelPrep 4に実装されている。

インストール

miniconda3.4.0.5環境でテストした（ホストOS、ubuntu16.0.4）。

依存

elPrep (since version 3.0) is implemented in Go.

elPrep have been tested with the following Linux distributions:

Ubuntu 14.04.5 LTS
Manjaro Linux
Red Hat Enterprise Linux 6.4 and 6.5
Amazon Linux 2

本体　GIthub

go get -u github.com/exascience/elprep
export PATH=$PATH:~/go/bin

#bioconda (link)
conda install -c bioconda -y elprep

> elprep --help-extended

# elprep --help-extended

elprep version 4.1.5 compiled with go1.12.7 - see http://github.com/exascience/elprep for more information.

Available commands: filter, split, merge, merge-optical-duplicates-metrics, sfm, vcf-to-elsites, bed-to-elsites, fasta-to-elfasta

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]

[--mark-optical-duplicates file]

[--optical-duplicates-pixel-distance nr]

[--remove-duplicates]

[--remove-optional-fields [all | list]]

[--keep-optional-fields [none | list]]

[--sorting-order [keep | unknown | unsorted | queryname | coordinate]]

[--clean-sam]

[--bqsr recal-file]

[--bqsr-reference elfasta]

[--quantize-levels nr]

[--sqq list]

[--max-cycle nr]

[--known-sites list]

[--nr-of-threads nr]

[--timed]

[--log-path path]

[--mark-optical-duplicates-intermediate file]

[--bqsr-tables-only table-file]

[--bqsr-apply path]

[--recal-file file]

split parameters:

elprep split (sam-file | /path/to/input/) /path/to/output/

[--output-prefix name]

[--output-type [sam | bam]]

[--single-end]

[--nr-of-threads nr]

[--timed]

[--log-path path]

[--contig-group-size nr]

merge parameters:

elprep merge /path/to/input sam-output-file

[--single-end]

[--nr-of-threads n]

[--timed]

[--log-path path]

sfm parameters:

elprep sfm 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]

[--mark-optical-duplicates file]

[--optical-duplicates-pixel-distance nr]

[--remove-duplicates]

[--remove-optional-fields [all | list]]

[--keep-optional-fields [none | list]]

[--sorting-order [keep | unknown | unsorted | queryname | coordinate]]

[--clean-sam]

[--bqsr]

[--bqsr-reference elfasta]

[--quantize-levels nr]

[--sqq list]

[--max-cycle nr]

[--known-sites list]

[--nr-of-threads nr]

[--timed]

[--log-path path]

[--intermediate-files-output-prefix name]

[--intermediate-files-output-type [sam | bam]]

[--tmp-path path]

[--single-end]

[--contig-group-size nr]

merge-optical-duplicates-metrics parameters:

elprep merge-optical-duplicates-metrics sam-input-file sam-output-file metrics-file /path/to/intermediate/metrics

[--remove-duplicates]

[--nr-of-threads nr]

[--timed]

[--log-path path]

vcf-to-elsites parameters:

elprep vcf-to-elsites vcf-file elsites-file

[--log-path path]

bed-to-elsites parameters:

elprep bed-to-elsites bed-file elsites-file

[--log-path path]

fasta-to-elfasta parameters:

elprep fasta-to-elfasta fasta-file elfasta-file

[--log-path path]

実行方法

version3のelprepは以前紹介しています。基本的な使い方はこちらも参考にして下さい。

elrepには大きく分けでfilter, split, merge, merge-optical-duplicates-metrics, sfm, vcf-to-elsites, bed-to-elsites, fasta-to-elfastaのモードがある。

下のコマンドを実行すると、論文中で記載されているのではGATKベストプラクティスのvariantコール前の煩雑な前処理パイプラインを１度に実行できる。

#BQSRする時はリファレンスが必要。fastaはあらかじめelprep用に変換して使用する。
elprep fasta-to-elfasta input.fasta output.elfasta

elprep sfm input.bam output.bam 
 ––mark–duplicates ––mark–optical–duplicates output.metrics
 ––sorting–order coordinate
 ––bqsr output.recal
 ––known–sites dbsnp_138.hg38.elsites
 ––bqsr–reference hg38.elfasta

入力と出力はそれぞれsuffix（sam/bam/cram）で判断する。入力はsamフォーマット、出力はbamフォーマットにするには、出力側のsuffixをxxx.bamとすればよい。

2021 9/20

BQSRする際に除外する既知バリアントサイトのリスト (vcf or bed) はelprepの形式に前もって変換しておく必要がある。

elprep vcf-to-elsites dbsnp_137.hg19.vcf dbsnp_137.hg19.elsites

さらにBQSRする際に必要なリファレンスfastaも、elprepの形式に前もって変換しておく必要がある。

elprep fasta-to-elfasta ucsc.hg19.fasta ucsc.hg19.elfasta

2020 7/24 追記

mappingしてelprepに渡しフィルタリング、coordinate sortしたbam出力する。unmapリードは除く。20スレッド指定。

#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 --filter-unmapped-reads-strict

#samtools indexing
samtools index filtered_sorted.bam
#sambamba indexing
sambamba index filtered_sorted.bam -t 10 -p

#2021 9/17 bedの領域以外は除く
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 --filter-unmapped-reads-strict \
 --filter-non-overlapping-reads capture.bed