macでインフォマティクス

macでインフォマティクス

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

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

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

2019 6/1 index追記

2019 71/17追記

2019 7/124 誤字修正

2019 8/3 help更新

2020 1/19 追記

2020 7/21 preset parameter追加

2021 1/17, 1/20 例追加

2021 7/3 構成変更

 

  

  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の多様性を実証する。

Minimap2は、ほとんどのフルゲノムアライナーで使用されている典型的なseed-chain-align手順に従っています。参照配列のミニマム化(Roberts et al., 2004)を収集し、ミニマム化のハッシュをキー、ミニマム化のコピーの位置のリストを値とするハッシュテーブルにインデックスします。次に、各クエリ配列に対して、minimap2は、クエリの最小化子をシードとし、参照配列との完全一致(すなわちアンカー)を見つけ、コリニアなアンカーのセットをチェーンとして識別します。ベースレベルのアライメントが要求された場合、minimap2はダイナミックプログラミング(DP)を適用して、チェーンの端から延長し、チェーン内の隣接するアンカー間の領域を閉じます。

minimap2は、minimap(Li, 2016)に類似したインデックス化とシード化のアルゴリズムを使用しており、より正確なチェイニング、ベースレベルアライメントの生成機能、スプライスされたアライメントのサポートなど、前身の機能をさらに強化している。

 

Minimap2は、ほとんどのフルゲノムアライナーで使用されている典型的なseed-chain-align手順に従っている。参照配列のminimizers (Roberts et al., 2004)を収集し、minimizers のハッシュをキー、ミニマム化のコピーの位置のリストを値とするハッシュテーブルにインデックスする。次に、各クエリ配列に対して、minimap2は、クエリのminimizers をシードとし、参照配列との完全一致(すなわちアンカー)を見つけ、コリニアなアンカーのセットをチェーンとして識別する。ベースレベルのアライメントが要求された場合、minimap2はダイナミックプログラミング(DP)を適用して、チェーンの端から延長し、チェーン内の隣接するアンカー間の領域を閉じる。

minimap2は、minimap(Li, 2016)に類似したインデックス化とシード化のアルゴリズムを使用しており、より正確なチェイニング、ベースレベルアライメントの生成機能、スプライスされたアライメントのサポートなど、前身の機能をさらに強化している。

 

ブログ

Minimap2 and the future of BWA

ブログより

”しかし、DNAnexusのAndrew Carroll氏は最近、彼の手による2つのNovaSeqの実行において、minimap2がbwa-memよりも遅いことを示してくれました。その理由の一つは、NovaSeqの2回の実行ではエラー率がやや高く、minimap2の高価なヒューリスティックがより頻繁に起動されるためだと思います。さらに、bwa-memは短いマッチに対してより敏感なので、Hi-Cアライメントではbwa-memの方がminimap2よりも優れているだろうとも考えています。

中略

minimap2は競争力のあるショートリードマッパーであると考えており、私の研究プロジェクトでは頻繁に使用します。しかし、様々な品質のショートリードに対してminimap2の性能がbwa-memほど安定していないことを考えると、少なくとも私がminimap2を改良する方法を見つけるまでは、生産的な用途ではbwa-memの方がまだ良いでしょう”。

 

2021 8/13

2021/5/27

2019/3/4

2018/5/10

 

 

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

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

Options:

  Indexing:

    -H           use homopolymer-compressed k-mer (preferrable for PacBio)

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

    -w INT       minimizer 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)

    -o FILE      output alignments to FILE [stdout]

    -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

    --eqx        write =/X CIGAR operators

    -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; see minimap2.1 for details) []

                 - map-pb/map-ont: PacBio/Nanopore vs reference mapping

                 - ava-pb/ava-ont: PacBio/Nanopore read overlap

                 - asm5/asm10/asm20: asm-to-ref mapping, for ~0.1/1/5% sequence divergence

                 - splice: long-read spliced alignment

                 - sr: genomic short-read mapping

 

See `man ./minimap2.1' for detailed description of these and other advanced 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するようにします。"-O BAM"はなくても機能します。


#sortコマンドが高速なsambambaに切り替える(フィルタリング機能がある、view, sortしているのでフィルタリングが不要ならsamtools sortよりかえって遅くなる可能性がある)
#sambamba sortはbamを受け付けないのでsambamba viewでbamにしてからsortする。
minimap2 -R "@RG\tID:X\tLB:Y\tSM:sample1\tPL:ILLUMINA" -t 8 -ax sr \
ref.fa pair1.fq.gz pair2.fq.gz \
| sambamba view -f bam -S /dev/stdin |sambamba sort -t 20 -o sorted.bam /dev/stdin

 

ーK”でチャンクサイズ(データのサブセットのサイズ)を変える。デフォルトは500Mとなっているが50Mの間違いと思われる。I/Oとの兼ね合いで少なくすると高速化する場合もある(少なくすればメモリ使用量は確実に減る)。

minimap2 -t3 -ax sr -K100M ref.fasta pair_*.fq.gz > out.sam

 

 

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.fasta

#any2fasta(紹介)を使うこともできる。
any2fasta assembly.gfa > assembly.fasta

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

 

nanopore RNA drect reads

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

 

preset  parameter

f:id:kazumaxneo:20200721142907p:plain

manualより

 

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

 

2019 11/24  追記

minimap2からワンライナーでバリアントコール

The Genome Factory

http://thegenomefactory.blogspot.com/2018/10/a-unix-one-liner-to-call-bacterial.html

 

2020 1/19 追記

minimap2のpaftoolsを使うと、長い配列を比較して変異コールすることができます。

 

all versus allリードマッピングにfpaのフィルタリングを取り込む。

引用

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

 

追記


関連