2018/12/21 ドラフトアセンブリ追記
2019 6/1 index追記、7/17追記、7/24 誤字修正、8/3 help更新
2020 1/19 追記、7/21 preset parameter追加
2021 1/17, 1/20 例追加、7/3 構成変更、10/9 新しい論文引用
2023/02/13 help更新, 07/06 追記
2024/02/15 分かりにくい説明を整頓
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 12/27
Minimap2 v2.24 released with alignment improvement in corner cases. https://t.co/gMlkKeUQfL
— Heng Li (@lh3lh3) December 26, 2021
2021 8/13
Minimap2 v2.22 released with a short preprint associated with it, more formally describing recent improvements since v2.19. Thank @chirgjain, @ArangRhie and all users for their invaluable feedback. https://t.co/POkFfpOYkW
— Heng Li (@lh3lh3) August 12, 2021
2021/5/27
Minimap2 v2.19 released with better and more contiguous alignment over long INDELs and in highly repetitive regions, improvements backported from unimap. These represent the most significant algorithmic change since v2.1. Use with caution. https://t.co/YiZU0QsUsm
— Heng Li (@lh3lh3) 2021年5月27日
2019/3/4
Talking about recent CPUs – the "avx" branch of minimap2 can use avx2 and avx512 for *non-splice* base alignment. You can use "make avx2=1" or "avx512=1" to compile. The speedup is marginal, so the change is not in master.
— Heng Li (@lh3lh3) 2019年3月4日
2018/5/10
Minimap2 published online at https://t.co/pYjXyN4hLJ. Preprint (not the latest): https://t.co/qzGFmMZY3A. LaTeX source: https://t.co/LhjeOcnNFM.
— Heng Li (@lh3lh3) 2018年5月10日
FAQ
https://github.com/lh3/minimap2/blob/master/FAQ.md
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 - lh3/minimap2: A versatile pairwise aligner for genomic and spliced nucleotide sequences
git clone https://github.com/lh3/minimap2
cd minimap2 && make
#conda
mamba install -c bioconda minimap2 -y
> ./minimap2 #v2.24-r1122
$ minimap2 -h
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[,NUM] chaining/alignment bandwidth and long-join bandwidth [500,20000]
-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 (larger value for lower divergence) [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 CLR/Nanopore vs reference mapping
- map-hifi - PacBio HiFi reads 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/splice:hq - long-read/Pacbio-CCS spliced alignment
- sr - genomic short-read mapping
ラン
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フォーマット出力になる。最近のsuper accuracy ONTリードでも動作するはずと述べられている(link)。
2、Pacific biosciencesのロングリード
#CLR
minimap2 -t 8 -ax map-pb ref.fa pacbio.fq > aln.sam
#HiFi(CCS)
minimap2 -t 8 -ax map-hifi 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、イルミナロングリード
minimap2 -ax map-iclr ref.fa iclr-reads.fq > aln.sam
5、ONTのRNA drect reads
minimap2 -ax splice -k14 -uf SIRV_E2.fa SIRV_ont-drna.fa > aln.sam
6、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
7、4Gより大きいラージゲノム
4Gbより長い参照塩基がある場合、minimap2はすべての配列を見ることができず、正しいSAMヘッダを出力しなくなってしまう(samtools viewでエラーが起きる)。2つの解決策があり、第一に、-Iオプションを例えば-I8gに増やして、より多くの参照塩基をインデックスすることができる(利用可能なメモリが十分なメモリがある場合)。十分なメモリが利用できない場合、--split-prefixを使う。この方法は一時的にディスク領域を使う。メモリ使用量は少なくなるが、処理速度は遅くなる。
メモリが潤沢なとき。-Iオプションの使用
#indexing
minimap2 -ax sr -I 8g -d large_genome.mmi large_genome.fna
#mapping
minimap2 -ax sr -t 40 large_genome.mmi pair1.fq.gz pair2.fq.gz > mapped.sam
メモリが潤沢に利用できない環境。--split-prefixオプションの使用
minimap2 -ax sr -t 8 --split-prefix=tmp large_genome.fna pair1.fq.gz pair2.fq.gz > mapped.sam
preset parameter
manualより
その他
ONTのロングリードをマッピングしてPAF出力
minimap2 -t 8 -x map-ont ref.fa nanopore.fq > mapping.paf
pafはstartとendの位置情報は持つが、1bp解像度のアライメント情報は持たない(miniasm PAF)。ゲノム間の相同性をHarr plotで解析したりするのに使う。
long RNAのマッピングについては、下記ページが詳しい。
リードグループを追加(*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
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等でコンセンサスコールを行い、精度を上げる。
minimap2からワンライナーでバリアントコール
The Genome Factory
http://thegenomefactory.blogspot.com/2018/10/a-unix-one-liner-to-call-bacterial.html
minimap2のpaftoolsを使うと、長い配列を比較して変異コールすることができる。minimap2のコンパニオンツールとしての性質を持つ。
all versus allリードマッピングにfpaのフィルタリングを取り込む。
引用
Minimap2: pairwise alignment for nucleotide sequences.
Li H
Bioinformatics. 2018 May 10.
publishされたため、preprintリンクからpubmedリンクに修正
2021 10/8
v2.22のパフォーマンスアップデートについての論文
https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btab705/6384570
参考ページ
*1
https://github.com/lh3/minimap2/issues/14
https://github.com/lh3/minimap2/commit/bbb37d95f207108acff6784b257171fdecfeb797
追記
関連
2021 12/11
With the recent (last year) advances in short-read aligners, it would be nice to see an independent benchmarking. I ran some of the tools in my small benchmarking, I include results and brief conclusions below. Could serve as inspiration for a more thorough benchmark. https://t.co/PHxBFLeA3a pic.twitter.com/iYBAI7vCiX
— Kristoffer Sahlin (@krsahlin) December 10, 2021
結論。正確さと実行時間の指標を評価する際に、%アラインドリードを見ることは重要である。難しいリード(繰り返しのある領域や、エラーや変動のためにレスキューが必要な領域)をスキップすることで、プログラムをより速く、より正確にすることができる(ツイートより)。
精度は最高でもマッピング率が低ければ実用的とは言えませんし、精度と感度の調和平均のような指標で表されるバランスが大事になります。
2022/05/26
"現在、ロングリードのマッピングに特化した数十のツールが存在し、そのヒューリスティックは時間とともにかなり進化している。この分野の急速な発展は、データの頻繁な改良と同期しており、文献や実装を追いかけることは容易ではない。そこで、本論文では、ロングリードのマッピング手法の概要とアルゴリズムの詳細について述べる。"
2022/07/07
T2T2アセンブリの論文。ONTリードからテロメア配列を組み立てるためにminimap2が使われています。テロメアリピートを持つロングリードを取り出し、minimap2によるAll-versus-allマッピングを最小クエリーカバレッジ95%を有する重複リードのみ行い、iGraphネットワークでキュレート、全ゲノムアセンブリと繋げて伸ばすという戦略が取られています。勉強になりました。
2022/07/30
M1 mac
Installing minimap2 on an m1 mac. The new apple silicon is aarch64 architecture and supports the neon instruction set. Took me longer than it should to put it together, given that it's in the README.
— Michael Barton (@bioinformatics) January 14, 2022
make aarch64=1 arm_neon=1
(1月のツイート)
2023/10/28