macでインフォマティクス

macでインフォマティクス

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

MinHashを利用した長い配列(ゲノムやロングリード)のアライナー MashMap

タイトル修正

 

 ハイスループットDNAシーケンサーによって生成されたリードをリファレンスゲノムにマッピングすることは、根本的かつ広く研究されている課題である[Preprintより ref.16,24]。この問題は、BWA [ref.15]やBowtie [ref.12]のようなマッピングアルゴリズムと広く使用されているショートリードシーケンスに対して特によく研究されている。 Pacific BiosciencesとOxford Nanoporeの1分子シーケンサーの人気が高まっており、シーケンスリードの長さ(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紹介

MinHashを使い高速にゲノムを比較する MASH

 

インストール

mac os 10.13でテストした。

Github

https://github.com/marbl/MashMap

#linuxmac 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/

> 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

 

-r <value>, --ref <value>

    an input reference file (fasta/fastq)[.gz]

 

--refList <value>, --rl <value>

    a file containing list of reference files, one per line

 

-q <value>, --query <value>

    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]

 

-t <value>, --threads <value>

    count of threads for parallel execution [default : 1]

 

-o <value>, --output <value>

    output file name [default : mashmap.out]

 

-k <value>, --kmer <value>

    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とリファレンスゲノムも非常に高速に行える。

出力は1リード1行のフォーマットになっている。クエリ名、リード長、start(0開始)、end、strand、ターゲット名、ターゲット長、ターゲットのstart、ターゲットのend、nucleotide identityとなっている。minimapの出力フォーマットに近い(PAF)。

 

 

複数のリファレンスにマッピング

mashmap --rl List.txt -q query.fa

リストファイルには、リファレンスのFASTAのパスを1行に1つずつ記載する。

 

パラメータに関しては、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 map is 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スクリプトを使う。

perl generateDotPlot png medium mashmap.out 

 

Preprintでは、ロングリードのhg38へのマッピングや、Refseq全ゲノム(6万以上)839Gbpへのマッピング実行時間等をBWA MEM(Li H et al 2013)やMinimap(Li H et al 2016)と比較している。同一条件でのヒトゲノムへのマッピング時間の比較から、BWA MEMよりMashMapが290倍早く、Refseq全ゲノムへのロングリード12万のマッピングも、ピークメモリ660GB、29CPU時間で終了している(メモリ使用率はPreprint 表2も参照(リンク))。 一方、BWAとminimapは、メモリ1TBのサーバーではリファレンスのインデックス作成ができなかった。

 

 

感想 

ONTの30万のロングリード(合計3.5GB)のリファレンスのA.thalianaゲノムへのマッピングタイムをは測ってみると、52秒だった(12スレッド使用)。同一条件でminimap2は6分19秒かかったので、7倍ほど高速という計算になる(マッピング精度もそれ以外のmetricsも調べてません)。

現状ではnucmer4(紹介)を高速化したツールのような扱いで、ロングリードのアライメントを行なっても、方法論からして当然ですが、アライメント内容は報告せず、アライメント境界のポジションとindentityを出力するだけです。今後、アライメントを全て報告し、samフォーマットに準じた出力が可能になれば、ロングリードのアライナーとしても活躍しそうですね。Githubのページに以下のようなコメントが記載されています。

"In future, we plan to add an optional alignment support to generate base-to-base alignments."

 

引用

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でも使われています。