macでインフォマティクス

macでインフォマティクス

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

高速なロング/ショートリードアライナー minimap2

2018/12/21 ドラフトアセンブリ追記

2019 6/1 index追記

 

 

  Single Molecule Real-Time(SMRT)シークエンシング技術とOxford Nanopore technologies(ONT)は、10kbp以上の長さのリードを約15%のエラー率で生成する。そのようなデータのためにいくつかのアライナーが開発されている(論文より Chaisson and Tesler、2012; Li、2013; Liu et al、2016; Sovic et al、2016; Liu et al、2017; Lin and Hsu、2017; Sedlazeck et al、2017)。それらの大半は、秒あたりマッピングされる塩基数の点で、主流のショートリードアライナー(Langmead and Salzberg、2012; Li、2013)の5倍の遅さであった。ショートリードアライメントのボトルネックであるリピート領域をより効果的にスキップすることができるため、10kbの配列は100bpのリードよりもマッピングが容易でなければならないという考えには、スピードアップの余地が大きいと推測された。著者らは、BWA-MEM(Li、2016)より約50倍速い近似マッピングを達成することによって、推測を確認した。 Suzuki and Kasahara(2018)(BMC)は、ベースレベルのアラインメントを生成するための高速かつ新規なアルゴリズムを使用してこの論文の著者の作業を拡張した。それがきっかけとなって、minimap2の開発が促された。

 SMRTおよびONTの両方が、スプライシングされたmRNA(RNA-seq)の配列決定に適用されている。伝統的なmRNAアライナーは機能するが(Wu and Watanabe、2005; Iwata and Gotoh、2012)、長時間のノイズの多いシーケンスリードには最適化されておらず、専用のロングリードアライナーよりも数十倍遅い。最初にゲノムDNAのみをアライメントさせるためのminimap2を開発する際には、わずかな改変がmRNAをマップするためのベースアルゴリズムを可能にすることが分かった。 Minimap2は、長いノイズの多いロングリード用に特別に設計された最初のRNA-seqアライナーになる。また、元のアルゴリズムを拡張して、いくつかの主流のショート・リード・マッパーよりも速い速度でショート・リードをマッピングする。

 この論文では、minimap2アルゴリズムとさまざまなタイプの入力シーケンスへのそのアプリケーションについて説明する。 いくつかのシミュレートされたデータセットとリアルデータセットでminimap2の性能と精度を評価し、minimap2の多様性を実証する。

 

 BWA(pubmed)でよく知られるHeng Liさん(wiki)(オーサー1人で有名な論文いくつも出しているすごい人です)が、自身のBlogでminimap2の論文のアクセプトを報告している。非常に興味深い内容になっている。どうやら"BWA2"も開発中らしい。

Minimap2 and the future of BWA

〜10kbのノイズの多いシーケンスの場合、minimap2はBLASR、BWA-MEM、NGMLR、GMAPなどより数十倍高速とされる。100bp以下のショートリードでは、例外はあるが(上記リンク参照)、BWA-MEMとBowtie2の3倍速いとされ、研究者からの評価もたいへん高い。 

twitter

 

DNAのショートリード、ロングリードのmappingだけでなく、long RNAのスプリットアライメントも可能で、RNAマッピングにも使えます。また、デフォルトはPAF出力で、ゲノム間の相同性をゲノムワイドに調べる用途にも使えます(LASTやMAMMERに似た使い方)。

 PAF format (D-GENIES manualより)

http://dgenies.toulouse.inra.fr/documentation/formats#index-file

 

 

minimap2 cookbook

 

インストール

Github

GitHub - lh3/minimap2: A versatile pairwise aligner for genomic and spliced nucleotide sequences

git clone https://github.com/lh3/minimap2 
cd minimap2 && make

> ./minimap2

$ ./minimap2 

Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]

Options:

  Indexing:

    -H           use homopolymer-compressed k-mer

    -k INT       k-mer size (no larger than 28) [15]

    -w INT       minizer window size [10]

    -I NUM       split index for every ~NUM input bases [4G]

    -d FILE      dump index to FILE

  Mapping:

    -f FLOAT     filter out top FLOAT fraction of repetitive minimizers [0.0002]

    -g NUM       stop chain enlongation if there are no minimizers in INT-bp [5000]

    -G NUM       max intron length (effective with -xsplice; changing -r) [200k]

    -F NUM       max fragment length (effective with -xsr or in the fragment mode) [800]

    -r NUM       bandwidth used in chaining and DP-based alignment [500]

    -n INT       minimal number of minimizers on a chain [3]

    -m INT       minimal chaining score (matching bases minus log gap penalty) [40]

    -X           skip self and dual mappings (for the all-vs-all mode)

    -p FLOAT     min secondary-to-primary score ratio [0.8]

    -N INT       retain at most INT secondary alignments [5]

  Alignment:

    -A INT       matching score [2]

    -B INT       mismatch penalty [4]

    -O INT[,INT] gap open penalty [4,24]

    -E INT[,INT] gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [2,1]

    -z INT[,INT] Z-drop score and inversion Z-drop score [400,200]

    -s INT       minimal peak DP alignment score [80]

    -u CHAR      how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]

  Input/Output:

    -a           output in the SAM format (PAF by default)

    -Q           don't output base quality in SAM

    -L           write CIGAR with >65535 ops at the CG tag

    -R STR       SAM read group line in a format like '@RG\tID:foo\tSM:bar'

    -c           output CIGAR in PAF

    --cs[=STR]   output the cs tag; STR is 'short' (if absent) or 'long' [none]

    --MD         output the MD tag

    -Y           use soft clipping for supplementary alignments

    -t INT       number of threads [3]

    -K NUM       minibatch size for mapping [500M]

    --version    show version number

  Preset:

    -x STR       preset (always applied before other options) []

                 map-pb: -Hk19 (PacBio vs reference mapping)

                 map-ont: -k15 (Oxford Nanopore vs reference mapping)

                 asm5: -k19 -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 (asm to ref mapping; break at 5% div.)

                 asm10: -k19 -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 (asm to ref mapping; break at 10% div.)

                 ava-pb: -Hk19 -Xw5 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)

                 ava-ont: -k15 -Xw5 -m100 -g10000 -r2000 --max-chain-skip 25 (ONT read overlap)

                 splice: long-read spliced alignment (see minimap2.1 for details)

                 sr: short single-end reads without splicing (see minimap2.1 for details)

 

See `man ./minimap2.1' for detailed description of command-line options.

 

ラン

1、ONTのロングリードのマッピング

minimap2 -t 8 -ax map-ont ref.fa nanopore.fq > aln.sam
  • -a       output in the SAM format (PAF by default)
  • -t       number of threads [3]
  •  -x <STR>   preset (always applied before other options) []
    map-pb: -Hk19 (PacBio vs reference mapping)
    map-ont: -k15 (Oxford Nanopore vs reference mapping)
    asm5: -k19 -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 (asm to ref mapping; break at 5% div.)
    asm10: -k19 -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 (asm to ref mapping; break at 10% div.)
    ava-pb: -Hk19 -Xw5 -m100 -g10000 --max-chain-skip 25 (PacBio read overlap)
    ava-ont: -k15 -Xw5 -m100 -g10000 -r2000 --max-chain-skip 25 (ONT read overlap)
    splice: long-read spliced alignment (see minimap2.1 for details)
    sr: short single-end reads without splicing (see minimap2.1 for details)

flagなしだとpafフォーマット出力になる。

 

2、Pacific biosciencesのロングリードのマッピング

minimap2 -t 8 -ax map-pb ref.fa pacbio.fq > aln.sam

 

3、ペアエンドショートリードのマッピング

#ペアエンド(消えていたので修正)
minimap2 -t 8 -ax sr ref.fa pair1.fq.gz pair2.fq.gz > aln.sam

#interleaveのペアエンド
minimap2 -t 8 -ax sr ref.fa interleaved.fq.gz > aln.sam

 

追記

4、ONTのロングリードをマッピングしてPAF出力。

minimap2 -t 8 -x map-ont ref.fa nanopore.fq > mapping.paf

 pafはstartとendの位置情報は持つが、1bp解像度のアライメント情報は持たない(miniasm PAF)。ゲノム間の相同性をHarr plotで解析したりするのに使う。

 long RNAマッピングについては、下記ページが詳しい。

CSC - minimap2 - Software

 

リードグループを追加(*1)。サンプル名はsample1とする。

minimap2 -R "@RG\tID:X\tLB:Y\tSM:sample1\tPL:ILLUMINA" -t 8 -ax sr ref.fa pair1.fq pair2.fq > aln.sam

#ヘッダーを確認
samtools view -H aln.sam

 

パイプしてbam出力。スレッドは8、メモリは2G / thread。

minimap2 -R "@RG\tID:X\tLB:Y\tSM:sample1\tPL:ILLUMINA" -t 8 -ax sr \
ref.fa pair1.fq.gz pair2.fq.gz \
|samtools sort -@ 8 -m 2G -O BAM - > sorted.bam \
&& samtools index -@ 8 sorted.bam

CRAM出力なら-O BAMを-O CRAMに変えて下さい。&& でindexを実行することでbam出力が成功した時にみindexするようにします。

 

ONTロングリードのドラフトアセンブリ

minimap2 -x ava-ont long-read.fq long-read.fq > ont.paf
miniasm -f long-read.fq ont.paf > assembly.gfa
awk '/^S/{print ">"$2"\n"$3}' assembly.gfa | fold > assembly.fa

このあとracon等でコンセンサスコールを行い、精度を上げる。

 

nanopore RNA drect reads

minimap2 -ax splice -k14 -uf SIRV_E2.fa SIRV_ont-drna.fa > aln.sam

 

 

2019 6/1 追記

indexを作成してからmapping

ヒトゲノムのラージゲノムのmampingを繰り返し行うなら、これによって時間短縮できる。 ただし、k-mer長などのパラメータはインデックス作成後には変更できないので注意する。

minimap2 -ax sr -d GRCh38_latest_genomic.mmi GRCh38_latest_genomic.fna
minimap2 -ax sr -t 40 GRCh38_latest_genomic.mmi pair1.fq.gz pair2.fq.gz > mapped.sam

 

 

引用

Minimap2: pairwise alignment for nucleotide sequences.

Li H

Bioinformatics. 2018 May 10.
publishされたため、preprintリンクからpubmedリンクに修正 

 

 参考ページ

CSC Services for Research

 

*1

https://github.com/lh3/minimap2/issues/14

https://github.com/lh3/minimap2/commit/bbb37d95f207108acff6784b257171fdecfeb797