macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

elPrep 5を使ったバリアントコール

 

 GATK Best Practices for variant callingに完全対応したelPrep5 (紹介) には、大きく分けて2つのモードが用意されています。1つ目は完全にRAM内で動作する(フィルタ)モードで、これは中間ファイルを全く書き出さず完全にRAM内で計算を進めるため、非常に高速です(GATK4で順番に進める場合の11−15倍高速)。2つ目は、ユーザーが提供したsam/bamファイルから実行するsfmモードで、このsfmモードでは、前もって中間ファイルであるsamやbamを作っておく必要がありますが、染色体領域ごとにデータを分割して進めるため、完全にRAM内で動作するフィルタモードよりもメモリ使用量を抑えることができます(GATK4で順番に進める場合の6-7.5倍高速)。利用可能な計算資源の量に対応して使い分けできるようになっているという事になりますGIthubの説明によると、sfmモードは、ヒトのWGS30xで128GBのメモリ使用量でランでき、ヒトWES30xでは24GBメモリ使用できるようです( WES30xのメモリモードだと96GGB使用)。elPrep5の論文ではm5.24xlargeインスタンスリンク)というかなりリッチな環境を使ってベンチマークなどを行なっていますが、このsfmモードなら、より少ない計算資源でも安定してランできそうです。

 

インストール

論文が出たからだと思いますが、ソースコードが公開され、condaで導入できるバージョンもv5になっています。

#anaconda (link)
mamba install -c bioconda elprep -y

 

 

実行方法

1、BQSR(参考サイト)を行うには、補正部位から除外する既知バリアントサイトのリスト (vcf or bed)も提供する必要がある。このvcf またはbed形式のファイルは、予めelprep vcf-to-elsites コマンドを走らせてelprepの形式に変換しておかないと使用できない。

ここでは時間の関係でdbSNP common SNV(生殖細胞由来、少なくとも1つの主要な集団においてマイナーアレル頻度(MAF)が>>0.01、少なくとも2人の血縁関係のない人がマイナーアレルを持っている)を使う。elprep vcf-to-elsitesコマンドは、入力のVCF、出力ファイルの順番に指定する。

elprep vcf-to-elsites dbSNP_common_all.vcf.gz  dbSNP_common_all.elsites

このコマンドは1回だけ実行すれば良いのですが、VCFが巨大だとかなりメモリを使います(最近のバージョンのdbSNP allだと300GB以上)。注意して下さい。

 

2、BQSRする際に必要なリファレンスfastaファイルも、elPrepの形式に変換しておく必要がある。elprep fasta-to-elfastaコマンドを使う。

elprep fasta-to-elfasta ucsc.GRCh38.fasta GRCh38.elfasta

 

3、elPrep5をランする。ここではメモリ使用量を抑えられるsfmモードを使う。sfmモードでは、ユーザーの方でアラインされたリードのsam かbam形式のファイルを提供する必要がある(論文中ではbwa memを使用)。

elprep sfm input.bam output.bam \
--mark-duplicates --mark-optical-duplicates NA12878.output.metrics\
--sorting-order coordinate \
--bqsr output.recal --known-sites dbSNP_common_all.elsites \
--reference GRCh38.elfasta --haplotypecaller output.vcf.gz \
--target-regions targetedregions.bed \
--filter-mapping-quality 0 \
--nr-of-threads 12 \

必要なら"--filter-mapping-quality 0"オプションも使います。 exome解析では、"--target-regions"オプションでbed形式ファイルを指定します。"filter-non-overlapping-reads”オプションを立てると、captureされた領域のbedファイルにマッチするリードだけがbamに保存され、バリアントコールの対象になります。バリアントコールは行わずにキャプチャ領域のリードだけ取り出したい時にも使えます。

注意

出力されるVCFはほぼすべてのポジションを含んでいます。バリアントコールに使うには、この後でフィルタリングする必要があります。

 

elPrep 5の論文では、NA12878を例にelPrep 5の使い方についても説明しています。また、GATKのオリジナルの Best Practicesのフローと比較して、リソース使用量を議論しています。簡単にまとめると、elPrepを使うことでディスク使用量、計算時間などを大幅に抑えられ、サーバーのレンタルコストを節約できるという結果です。論文にアクセスしてみてください。

引用

Multithreaded variant calling in elPrep 5

Charlotte Herzeel  ,Pascal Costanza ,Dries Decap,Jan Fostier,Roel Wuyts,Wilfried Verachtert

PLOS ONE, Published: February 4, 2021

 

参考


関連