2018 タイトル修正
2019 6/21 インストール追記、論文追記
2020 10/19 help更新、文章修正
ハイスループットDNAシーケンサーによって生成されたリードをリファレンスゲノムにマッピングすることは、根本的かつ広く研究されている課題である[Preprintより ref.16,24]。この問題は、BWA [ref.15]やBowtie [ref.12]のようなマッピングアルゴリズムと広く使用されているショートリードシーケンスに対して特によく研究されている。 Pacific BiosciencesとOxford Nanoporeの１分子シーケンサーの人気が高まっており、シーケンスリードの長さ（10 kbp以上）が絶えず向上しているため、ロングリードのマッピングアルゴリズムに関心が高まっている。しかし、ロングリードは、現在のところ、非常に高いエラー率（最大15-20％）を伴う。高いエラー率にもかかわらず、デノボゲノムアセンブリ[ref.7,10]およびリアルタイム病原体同定[ref.2,22]を含む多くのアプリケーションにおいて、ロングリードが有利であることが判明している。
ナノポアデバイスからのシーケンスデータは、サンプルを導入した直後に入手できる。これにより、大規模なリファレンスデータベースに対してデータストリームをマッピングできる高速の計算方法と組み合わせると、リアルタイムにゲノム解析を行うことが可能になる。しかし、未加工のシーケンスをマッピングすることは、引き続き多くのアプリケーションにとってボトルネックになる。この問題は、Oxford NanoporeのPromethIONが1日あたり数テラの塩基配列を生成すると予想されるため、より悪化するだろう。並行して、リファレンスデータベースのサイズは絶えず拡大しており、non-redundant NCBI RefSeqデータベースはテラベース（塩基）のサイズを超えている。1分子シーケンスの高いエラー率は、計算上の複雑さをさらに増す。
マッピングの問題は、Smith-Watermanアライメントアルゴリズム[ref.27]の適切なバリエーションを設計することによって正確に解決することができる。しかし、ハイスループットシーケンサーから大きなリファレンスゲノムにマッピングすることは、計算上、禁止されている。 Seed-and-extendマッピングのヒューリスティックスは、アラインメントアルゴリズム実行前に、正確な短いwordまたは最大の共通部分文字列一致が発生する場所に検索を限定することによって、この制限に対処している[ref.1,8,5 ]。 BLASR [ref.5]とBWA-MEM [ref.13]は、正確な位置合わせに基づくロングリードマッパーに含まれている。しかし、高いシークエンシングエラー率により、正しいマッピングに変換されない反復シードがスケーラビリティを制限する。さらに、アライメントベースのマッピングアルゴリズムは、インデックス内の完全なリファレンスシーケンスを保存するため、テラベースのリファレンスデータベースに拡大することはできない。多くのゲノミクスアプリケーションでは、詳細な塩基間のアライメント情報は必要なく、代わりにアラインメント境界とアイデンティティサマリーのみを使用する。そのようなアプリケーションには、カバレッジデプス分析、メタゲノミックリードのアサイン、構造変化検出、selective sequencing[ref.18]などがある。これらの問題の効率的なアルゴリズムをナノポアシーケンスと組み合わせることで、患者、病原体、ガン、および微生物のリアルタイムゲノム解析が可能になる。
高速で近似的なマッピングのためのアルゴリズムの1つのクラスは、Webドキュメント間の類似点を見つけるために元々開発されたアイデアに依存している。 Broder [ref.4]は、2つのセット間のJaccard類似性係数の不偏推定値が、sketchと呼ばれるハッシュされた要素のサブセットを使用して効率的に計算できることを証明した 。Schleimer et alは Web文書間の局所的類似性をより迅速に推定する手段として、テキストの各連続ウィンドウから最小ハッシュ項目（minimizer[23]としても知られている）を選択するwinnowing algorithmを提案した[25]。これらのアイデアは、MinHash Alignment Process [ref.3]、minimap [ref.14]、BALAUR [ref.21]などのロングリードための新しいマッピングとアセンブリアルゴリズムを開発するために使用されてきた。今日まで、これらのアプローチの有効性は、empiricallyでのみ実証されている。
本稿では、十分な理論的保証と実証的検証を用いて、ラージリファレンスデータベースにロングリードをマッピングするための高速近似アルゴリズムを提案する。（一部略）。この理論は、PacBioおよびMin-IONデータセットを使用して検証され、PacBioメタゲノミックシーケンスリードをRefSeqデータベース全体にマッピングすることによって、スケーラビリティを実証する。著者らのアルゴリズムの速度と空間効率はリアルタイムマッピングを可能にし、minimapと比較して、本方法は、大きな反復ゲノムに対してより高い精度で高い感度を維持する。実装はgithub.com/MarBL/MashMapから入手できる。
MASHのラストオーサAdam M. Phillippyさんの解説。
MashMap: approximate long-read mapping using minimizers and MinHash – Genome Informatics Section
MASH紹介
インストール
ubuntu16.04のminiconda3-4.3.30環境でテストした。
https://github.com/marbl/MashMap
#<old> linuxとmac os向けバイナリがリリースに準備されている。mac版をダウンロードする。
wget https://github.com/marbl/MashMap/releases/download/v2.0/mashmap-OSX64-v2.0.tar.gz
tar -zxvf mashmap-OSX64-v2.0.tar.gz
chmod u+x mashmap-OSX64-v2.0/mashmap
#パスの通ったディレクトリに移動。
mv mashmap-OSX64-v2.0/mashmap /usr/local/bin/
#bioconda (link)#現在v3が最新
mamba install -c bioconda -y mashmap
> mashmap
$ mashmap -h
-----------------
Mashmap is an approximate long read or contig mapper based on Jaccard
similarity
-----------------
Example usage:
$ mashmap -r ref.fa -q seq.fq [OPTIONS]
$ mashmap --rl reference_files_list.txt -q seq.fq [OPTIONS]
Available options
-----------------
-h, --help
Print this help page
an input reference file (fasta/fastq)[.gz]
--refList <value>, --rl <value>
a file containing list of reference files, one per line
an input query file (fasta/fastq)[.gz]
--ql <value>, --queryList <value>
a file containing list of query files, one per line
-s <value>, --segLength <value>
mapping segment length [default : 5,000]
sequences shorter than segment length will be ignored
--noSplit
disable splitting of input sequences during mapping [enabled by default]
--perc_identity <value>, --pi <value>
threshold for identity [default : 85]
count of threads for parallel execution [default : 1]
output file name [default : mashmap.out]
kmer size <= 16 [default : 16]
-f <value>, --filter_mode <value>
filter modes in mashmap: 'map', 'one-to-one' or 'none' [default: map]
'map' computes best mappings for each query sequence
'one-to-one' computes best mappings for query as well as reference sequence
'none' disables filtering
ラン
mashmap -r reference.fna -q query.fa -o mashmap.out -t 8
- -s mapping segment length [default : 5,000]. Sequences shorter than segment length will be ignored.
- -t count of threads for parallel execution [default : 1]
- --pi Threshold for identity [default : 85]
mashmap.outができる。ロングリードだけでなく、 contigとリファレンスゲノムも非常に高速に行える。
出力は１リード１行のフォーマットになっている。クエリ名、リード長、start（0開始）、end、strand、ターゲット名、ターゲット長、ターゲットのstart、ターゲットのend、nucleotide identityとなっている。minimapの出力フォーマットに近い（PAF）。
複数のリファレンスにマッピング。
mashmap --rl List.txt -q query.fa
リストファイルには、リファレンスのFASTAのパスを１行に１つずつ記載する。
パラメータに関しては、Githubに以下のようなコメントが記載されている。
-
Identity threshold (--perc_identity, --pi) : By default, it is set to 85, implying mappings with 85% or more identity should be reported. For example, it can be set to 80% to account for more noisy long-read datasets or 95% for mapping human genome assembly to human reference.
-
Minimum segment length (-s, --segLength) : Default is 5,000 bp. Sequences below this length are ignored. Mashmap provides guarantees on reporting local alignments of length twice this value.
-
Filtering options (-f, --filter_mode) : Mashmap implements a plane-sweep based algorithm to perform the alignment filtering. Similar to delta-filter in nucmer, different filtering options are provided that are suitable for long read or assembly mapping. Option
-f mapis suitable for reporting the best mappings for long reads, whereas-f one-to-oneis suitable for reporting orthologous mappings among all computed assembly to genome mappings.
perl generateDotPlot png medium mashmap.out
Githubより
Preprintでは、ロングリードのhg38へのマッピングや、Refseq全ゲノム（６万以上）839Gbpへのマッピング実行時間等をBWA MEM（Li H et al 2013）やMinimap（Li H et al 2016）と比較している。同一条件でのヒトゲノムへのマッピング時間の比較から、BWA MEMよりMashMapが290倍早く、Refseq全ゲノムへのロングリード12万のマッピングも、ピークメモリ660GB、29CPU時間で終了している（メモリ使用率はPreprint 表２も参照（リンク））。 一方、BWAとminimapは、メモリ1TBのサーバーではリファレンスのインデックス作成ができなかった。
感想
ONTの30万のロングリード（合計3.5GB）のリファレンスのA.thalianaゲノムへのマッピングタイムをは測ってみると、52秒だった（12スレッド使用）。同一条件でminimap2は6分19秒かかったので、７倍ほど高速という計算になる（マッピング精度もそれ以外のmetricsも調べてません）。
現状では配列jの類似性を求めるための近似の実装であり、アライメントは行っていません（アライメント境界のポジションとindentityは出力される）。しかし、今後、samフォーマットに準じた出力が可能になれば、ロングリードのアライナーとしても活躍しそうですね。Githubのページに以下のようなコメントが記載されています。
"In future, we plan to add an optional alignment support to generate base-to-base alignments."
引用
A fast adaptive algorithm for computing whole-genome homology maps
Chirag Jain, Sergey Koren, Alexander Dilthey, Adam M Phillippy, Srinivas Aluru
Bioinformatics, Volume 34, Issue 17, 01 September 2018, Pages i748–i756,
A Fast Approximate Algorithm for Mapping Long Reads to Large Reference Databases
Jain Chirag , Dilthey Alexander , Koren Sergey , Aluru Srinivas , and Phillippy Adam M.
Published Online: 30 Apr 2018 https://doi.org/10.1089/cmb.2018.0036
https://www.liebertpub.com/doi/pdf/10.1089/cmb.2018.0036
A Fast Adaptive Algorithm for Computing Whole-Genome Homology Maps
Chirag Jain, Sergey Koren, Alexander Dilthey, Adam M. Phillippy, and Srinivas Aluru.
BioRxiv, 2018.
https://www.biorxiv.org/content/biorxiv/early/2017/03/24/103812.full.pdf
A fast approximate algorithm for mapping long reads to large reference databases
Chirag Jain, Alexander Dilthey, Sergey Koren, Srinivas Aluru, and Adam M. Phillippy.
In International Conference on Research in Computational Molecular Biology, Springer, Cham, 2017.
https://github.com/marbl/MashMap/issues/9
MashmapはD-GENIESでも使われています。
