macでインフォマティクス

macでインフォマティクス

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

ショートリードのマッピングを行う Whisper

 

 リファレンスゲノムへのリードのマッピングは、シークエンシングデータ解析パイプラインの最初のステップである。シーケンシングコストが削減していることから、合理的な時間内に増大する量の生成データを処理することができるアルゴリズムに対する必要性が高まっている。
 著者らはリードをソートし、それからリファレンスゲノムとそのreverse complementのsuffix arraysに対してそれらをマッピングするという考えに基づいて、正確で高性能なマッピングツールであるWhisperを紹介する。 タスクとデータの並列処理を採用し、一時データをディスクに格納すると、妥当なメモリ要件で優れた時間効率が得られる。 Whisperは、大規模なNGSリードコレクション、特にIlluminaの一般的なWGSカバレッジのデータ処理で優れた性能を示す。 実データによる実験では、よく知られているBWA-MEMおよびBowtie 2ツールで必要とされる時間の約15%で、同等の精度でこのソリューションが機能することが示されている。

 

 

本体 Github

リリースよりダウンロードする

https://github.com/refresh-bio/Whisper/releases

#linux v1.1
wget https://github.com/refresh-bio/Whisper/releases/download/v1.1/whisper
chmod +x whisper
sudo mv whisper /usr/local/bin/

#linux
wget https://github.com/refresh-bio/Whisper/releases/download/v1.1/whisper-index
chmod +x whisper-index
sudo mv whisper-index /usr/local/bin/

whisper

$ whisper

Whisper v. 1.1 (2018-07-10)

Usage:

   whisper [options] <index_name> @<files> 

   whisper [options] <index_name> file_se 

   whisper [options] <index_name> file_pe1 file_pe2

Parameters:

  index_name   - name of the index (as created by asm_pp)

  files        - name of the file containing list of FASTQ files with seq. reads

  file_se      - FASTQ file (single-end)

  file_pe[1|2] - FASTQ files (paired-end)

Options:

  -b <value> - no. of temporary files (minimum: 100, default: 384)

  -d[fr/ff/rf] - mapping orientation (default: -dfr (forward - reverse)

  -dist_paired <value> - max. distance for paired read (default: 1000)

  -e <value> - max. no of errors (default: 4, max: 5% of read length)

  -e-paired <value> - max. fraction of errors in paired read (default: 0.06)

  -enable-boundary-clipping <value> - enable clipping at boundaries when a lot of mismatches appears (default: 0)

  -filter <value> - store only mappings for given chromosome (default: )

  -gap-open <value> - score for gap open (default: -6)

  -gap-extend <value> - score for gap extend (default: -1)

  -gzipped-SAM-level <value> - gzip compression level of SAM/BAM, 0 - no compression (default: 0)

  -hit-merging-threshold <value> - minimal distance between different mappings (default: 12)

  -high-confidence-sigmas <value> - (default: 4)

  -hit-merging-wrt-first <value> - calculate distance in marged group w.r.t. first (default: 1)

  -m[f/s/a] - mode: first stratum/second stratum/all strata (default: first stratum)

  -mask-lqb <value> - mask bases of quality lower than value (default: 0)

  -out <name> - name of the output file (default: whisper)

  -penalty-saturation <value> - no. of sigmas for max. penalty in matching pairs (default: 7)

  -rg <read_group> - complete read group header line, (example: '@RG\tID:foo\tSM:bar')

                     '\t' character will be converted to a TAB in the output SAM/BAM,

                      while the read group ID will be attached to every read 

  -r[s|p] - single or paired-end reads <value> (default: single)

  -score-discretization-threshold (default: 0.5)

  -score-match <value> - score for matching symbol (default: 1)

  -score-clipping <value> score for clipping (default: -10)

  -score-mismatch <value> - score for mismatching symbol (default: -5)

  -sens <value> - turn on/off sensitive mode (default: 1)

  -sens-factor <value> - sensitivity factor (default: 3)

  -stdout - use stdout to store the output

  -store-BAM - turn on saving in BAM

  -t <value> - no. of threads (0-adjust to hardware) (default: 0)

  -temp <name> - prefix for temporary files (default: ./whisper_temp_)

  -x - load complete suffix arrays in main memory (default: 0)

Examples:

  whisper human @files

  whisper -temp temp/ human reads1.fq reads2.fq

  whisper -out result.sam -temp temp/ -t 12 human reads1.fq reads2.fq

kazu@kazu:~$ 

 

> whisper-index

$ whisper-index 

Whisper index construction v. 1.1 (2018-07-10)

Usage: 

   whisper-index <index_name> <ref_seq_file_name> <dest_dir> <temp_dir>

   whisper-index <index_name> <@ref_seq_files_name> <dest_dir> <temp_dir>

 

 

実行方法

たとえば、圧縮されていないFASTQ形式のサンプル読み取りのサイズが100GBで、処理が16個のCPUスレッドによって行われる場合、Whisperは600のファイルを開く。そのため、同時に開けるファイル数を上げる必要がある。linuxなら

ulimit -n 600

fastqサイズや使用スレッド数が増えると増加する。 必要な数は

num_files = num_bins (384 by default) + num_threads + 2 * total_size_of_FASTQ_in_GB

で算出できる。

 

1、indexing

出力と作業ディレクトリを指定する。ここでは作業ディレクトリは/tmpか/var/tmpあたりにする。

mdkir output_dir
whisper-index genome_index ref.fa output_dir /tmp

 

2、mapping

ペアエンドfastqとindexを指定する。

whisper -out output -t 8 output_dir/genome_index \
pair_1.fq.gz pair_2.fq.gz

output.samが出力される。

 

引用

Whisper: read sorting allows robust mapping of DNA sequencing data
Sebastian Deorowicz Agnieszka Debudaj-Grabysz Adam Gudyś Szymon Grabowski
Bioinformatics, Volume 35, Issue 12, June 2019, Pages 2043–2050