macでインフォマティクス

macでインフォマティクス

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

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

mapping Pacbio Nanopore long read 高速なツール format 2018 Bioinformatics polish Nanopore direct RNA-seq

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

2021 8/13

2021/5/27

2019/3/4

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

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）。

 

２、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

 

３、ペアエンドショートリード

#ペアエンド(消えていたので修正)
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

 

４、イルミナロングリード

minimap2 -ax map-iclr ref.fa iclr-reads.fq > aln.sam  

 

５、ONTのRNA drect reads

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

 

６、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

 

７、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

f:id:kazumaxneo:20200721142907p:plain

manualより

 

 

その他

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

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

 pafはstartとendの位置情報は持つが、１bp解像度のアライメント情報は持たない（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等でコンセンサスコールを行い、精度を上げる。

 

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

 

 参考ページ

CSC Services for Research

 

*1

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

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

 

追記


関連

 

2021 12/11

結論。正確さと実行時間の指標を評価する際に、％アラインドリードを見ることは重要である。難しいリード（繰り返しのある領域や、エラーや変動のためにレスキューが必要な領域）をスキップすることで、プログラムをより速く、より正確にすることができる(ツイートより)。

精度は最高でもマッピング率が低ければ実用的とは言えませんし、精度と感度の調和平均のような指標で表されるバランスが大事になります。

 

2022/05/26

"現在、ロングリードのマッピングに特化した数十のツールが存在し、そのヒューリスティックは時間とともにかなり進化している。この分野の急速な発展は、データの頻繁な改良と同期しており、文献や実装を追いかけることは容易ではない。そこで、本論文では、ロングリードのマッピング手法の概要とアルゴリズムの詳細について述べる。"


2022/07/07

T2T2アセンブリの論文。ONTリードからテロメア配列を組み立てるためにminimap2が使われています。テロメアリピートを持つロングリードを取り出し、minimap2によるAll-versus-allマッピングを最小クエリーカバレッジ95％を有する重複リードのみ行い、iGraphネットワークでキュレート、全ゲノムアセンブリと繋げて伸ばすという戦略が取られています。勉強になりました。


2022/07/30

M1 mac

(１月のツイート)

 

2023/10/28

https://twitter.com/lh3lh3/status/1716295082259697846