macでインフォマティクス

macでインフォマティクス

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

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

 

 BLASTが1990年に初めてpublishされたとき、公開されたアーカイブには5000万塩基以下の塩基配列しか存在しなかった[論文より ref.2]。現在では、1つのシーケンシング機器1回の実行で1兆塩基を超えるシーケンス生成が可能である[ref.3]。この規模のデータを管理し、整理する新しい方法が必要である。これに対処するために、2つのシーケンス間の近似距離を計算する一般的な問題を検討し、大きなシーケンス(またはシーケンスセット)を圧縮Sketch表現に縮小するMinHash技法[ref.4]を使用する汎用ツールキットであるMashについて説明する。数千倍も小さいSketchのみを使用することで、元のシーケンスの類似性を有界誤差で素早く推定することができる。重要なことに、この計算の誤差はSketchのサイズにのみ依存し、ゲノムサイズとは無関係である。したがって、わずか数百の値を含むSketchを使用して、任意に大きなデータセットの類似性を近似することができる。これは、大規模なゲノムデータ管理と、長らく読まれている単一分子シークエンシング技術の重要なアプリケーションである。潜在的なアプリケーションには、おおよそのグローバル距離が許容可能な問題が含まれる。種のラベルを割り当て、大きなガイドツリーを構築し、誤って追跡されたサンプルを特定し、ゲノムデータベースを検索することができる。

 MinHash技法は、ほぼ重複するWebページや画像の検出に広く使用されているローカリティセンシティブハッシング[ref.5]の一形態であるが、10年以上前の初期アプリケーションにもかかわらずゲノミクスでの使用は限られている[ref.8]。最近、MinHashはゲノムアセンブリ[ref.9]、16S rDNA遺伝子クラスタリング[ref.10]、およびメタゲノムのクラスタリング[ref.12]等の関連問題に適用されている。この確率的アプローチのメモリとCPUの要求が非常に低いため、MinHashはゲノミクスのデータ集約的な問題に適している。これを容易にするために、著者らは柔軟な構造、操作、およびゲノムデータからのMinHashスケッチの比較のためのMashを開発した。過去のMinHashのアプリケーションをベースに、データベースの検索時にチャンスマッチを区別する新しい重要性テストを導き、MinHashスケッチから2つのシーケンス間の突然変異率を推定する新しい距離メトリクスMash距離を導き出す。同様の「アライメントフリー」法は、バイオインフォマティクスで長い歴史を持っている[ref.13,14]。しかしながら、word数に基づく従来の方法は、数ヌクレオチドの短いwordに依存しており、これは密接に関連する配列を区別することができず、解釈が難しい距離測定値を生成する[ref.15-18]。あるいは、ストリングマッチングに基づく方法は、非常に正確な突然変異距離の推定値を生成することができるが、すべてのペアを比較してシーケンス全体を処理しなくてはならない[ref.19-22]。対照的に、Mash距離はサイズ縮小されたスケッチだけから迅速に計算できるが、平均ヌクレオチドアイデンティティ(ANI)[ref.23]などのアライメントベースの測定値と強く相関する結果が得られる。したがって、マッシュは、マッチングに基づくアプローチの高い特異性と統計的アプローチの次元削減とを組み合わせ、多くの大きなゲノムとメタゲノムの間の正確な全ペア比較を可能にする。

 Mashはシーケンス比較のための2つの基本的な機能を提供する:Sketchとdist。Sketch関数は、シーケンスまたはシーケンスのコレクションをMinHash Sketchに変換する(論文 図1)。 dist関数は2つのスケッチを比較し、Jaccardインデックス(すなわち共有されるk-merの割合)、P value、単純な進化モデルによる配列突然変異の割合を推定するMash distanceを返す[ref.22] 「Method」を参照)。 Mashは長さkのサブストリングまたはk-mersの比較にのみ依存するため、入力は全ゲノム、メタゲノム、ヌクレオチド配列、アミノ酸配列、または生シーケンシングの読み込みである可能性がある。各入力は、既知のアルファベットから取られたk-merの集合として単純に扱われ、多くのアプリケーションを可能にする。ここでは、3つの具体的な使用例を検討する。(1)NCBI RefSeqゲノムデータベース全体のスケッチとクラスタリング。 (2)リアルタイムでスケッチされたRefSeqデータベースに対してアセンブリされた、およびアセンブリされていないゲノムを検索すること。 (3)アセンブリされたリードセットとされていないリードセットの両方を用いてメタゲノミックサンプル間の距離を計算すること、となる。

  

マニュアル

http://mash.readthedocs.io/en/latest/

  

インストール

mac os 10.12に導入した。

本体 Github

https://github.com/marbl/mash

バイナリリリース

https://github.com/marbl/Mash/releases

osx向けのバイナリ(version2)をダウンロードし、パスの通ったディレクリに移動しておく。

> mash

 mash

 

Mash version 2.0

 

Type 'mash --license' for license and copyright information.

 

Usage:

 

  mash <command> [options] [arguments ...]

 

Commands:

 

  bounds  Print a table of Mash error bounds.

 

  dist    Estimate the distance of query sequences to references.

 

  info    Display information about sketch files.

 

  paste   Create a single sketch file from multiple sketch files.

 

  screen  Determine whether query sequences are within a larger pool of

          sequences.

 

  sketch  Create sketches (reduced representations for fast operations).

 

> mash sketch

$ mash sketch

 

Version: 2.0

 

Usage:

 

  mash sketch [options] <input> [<input>] ...

 

Description:

 

  Create a sketch file, which is a reduced representation of a sequence or set

  of sequences (based on min-hashes) that can be used for fast distance

  estimations. Inputs can be fasta or fastq files (gzipped or not), and "-" can

  be given to read from standard input. Input files can also be files of file

  names (see -l). For output, one sketch file will be generated, but it can have

  multiple sketches within it, divided by sequences or files (see -i). By

  default, the output file name will be the first input file with a '.msh'

  extension, or 'stdin.msh' if standard input is used (see -o).

 

Options:

 

  Option     Description (range) [default]

 

  -h         Help

 

  -p <int>   Parallelism. This many threads will be spawned for processing. [1]

 

...Input...

 

  -l         List input. Lines in each <input> specify paths to sequence files,

             one per line.

 

...Output...

 

  -o <path>  Output prefix (first input file used if unspecified). The suffix

             '.msh' will be appended.

 

...Sketching...

 

  -k <int>   K-mer size. Hashes will be based on strings of this many

             nucleotides. Canonical nucleotides are used by default (see

             Alphabet options below). (1-32) [21]

 

  -s <int>   Sketch size. Each sketch will have at most this many non-redundant

             min-hashes. [1000]

 

  -i         Sketch individual sequences, rather than whole files, e.g. for

             multi-fastas of single-chromosome genomes or pair-wise gene

             comparisons.

 

  -S <int>   Seed to provide to the hash function. (0-4294967296) [42]

 

  -w <num>   Probability threshold for warning about low k-mer size. (0-1)

             [0.01]

 

  -r         Input is a read set. See Reads options below. Incompatible with -i.

 

...Sketching (reads)...

 

  -b <size>  Use a Bloom filter of this size (raw bytes or with K/M/G/T) to

             filter out unique k-mers. This is useful if exact filtering with -m

             uses too much memory. However, some unique k-mers may pass

             erroneously, and copies cannot be counted beyond 2. Implies -r.

 

  -m <int>   Minimum copies of each k-mer required to pass noise filter for

             reads. Implies -r. [1]

 

  -c <num>   Target coverage. Sketching will conclude if this coverage is

             reached before the end of the input file (estimated by average

             k-mer multiplicity). Implies -r.

 

  -g <size>  Genome size (raw bases or with K/M/G/T). If specified, will be used

             for p-value calculation instead of an estimated size from k-mer

             content. Implies -r.

 

...Sketching (alphabet)...

 

  -n         Preserve strand (by default, strand is ignored by using canonical

             DNA k-mers, which are alphabetical minima of forward-reverse

             pairs). Implied if an alphabet is specified with -a or -z.

 

  -a         Use amino acid alphabet (A-Z, except BJOUXZ). Implies -n, -k 9.

 

  -z <text>  Alphabet to base hashes on (case ignored by default; see -Z).

             K-mers with other characters will be ignored. Implies -n.

 

  -Z         Preserve case in k-mers and alphabet (case is ignored by default).

             Sequence letters whose case is not in the current alphabet will be

             skipped when sketching.

 

 > mash dist

$ mash dist

 

Version: 2.0

 

Usage:

 

  mash dist [options] <reference> <query> [<query>] ...

 

Description:

 

  Estimate the distance of each query sequence to the reference. Both the

  reference and queries can be fasta or fastq, gzipped or not, or Mash sketch

  files (.msh) with matching k-mer sizes. Query files can also be files of file

  names (see -l). Whole files are compared by default (see -i). The output

  fields are [reference-ID, query-ID, distance, p-value, shared-hashes].

 

Options:

 

  Option     Description (range) [default]

 

  -h         Help

 

  -p <int>   Parallelism. This many threads will be spawned for processing. [1]

 

...Input...

 

  -l         List input. Lines in each <query> specify paths to sequence files,

             one per line. The reference file is not affected.

 

...Output...

 

  -t         Table output (will not report p-values, but fields will be blank if

             they do not meet the p-value threshold).

 

  -v <num>   Maximum p-value to report. (0-1) [1.0]

 

  -d <num>   Maximum distance to report. (0-1) [1.0]

 

...Sketching...

 

  -k <int>   K-mer size. Hashes will be based on strings of this many

             nucleotides. Canonical nucleotides are used by default (see

             Alphabet options below). (1-32) [21]

 

  -s <int>   Sketch size. Each sketch will have at most this many non-redundant

             min-hashes. [1000]

 

  -i         Sketch individual sequences, rather than whole files, e.g. for

             multi-fastas of single-chromosome genomes or pair-wise gene

             comparisons.

 

  -S <int>   Seed to provide to the hash function. (0-4294967296) [42]

 

  -w <num>   Probability threshold for warning about low k-mer size. (0-1)

             [0.01]

 

  -r         Input is a read set. See Reads options below. Incompatible with -i.

 

...Sketching (reads)...

 

  -b <size>  Use a Bloom filter of this size (raw bytes or with K/M/G/T) to

             filter out unique k-mers. This is useful if exact filtering with -m

             uses too much memory. However, some unique k-mers may pass

             erroneously, and copies cannot be counted beyond 2. Implies -r.

 

  -m <int>   Minimum copies of each k-mer required to pass noise filter for

             reads. Implies -r. [1]

 

  -c <num>   Target coverage. Sketching will conclude if this coverage is

             reached before the end of the input file (estimated by average

             k-mer multiplicity). Implies -r.

 

  -g <size>  Genome size (raw bases or with K/M/G/T). If specified, will be used

             for p-value calculation instead of an estimated size from k-mer

             content. Implies -r.

 

...Sketching (alphabet)...

 

  -n         Preserve strand (by default, strand is ignored by using canonical

             DNA k-mers, which are alphabetical minima of forward-reverse

             pairs). Implied if an alphabet is specified with -a or -z.

 

  -a         Use amino acid alphabet (A-Z, except BJOUXZ). Implies -n, -k 9.

 

  -z <text>  Alphabet to base hashes on (case ignored by default; see -Z).

             K-mers with other characters will be ignored. Implies -n.

 

  -Z         Preserve case in k-mers and alphabet (case is ignored by default).

             Sequence letters whose case is not in the current alphabet will be

             skipped when sketching.

mash screen

$ mash screen

 

Version: 2.0

 

Usage:

 

  mash screen [options] <queries>.msh <pool> [<pool>] ...

 

Description:

 

  Determine how well query sequences are contained within a pool of sequences.

  The queries must be formatted as a single Mash sketch file (.msh), created

  with the `mash sketch` command. The <pool> files can be contigs or reads, in

  fasta or fastq, gzipped or not, and "-" can be given for <pool> to read from

  standard input. The <pool> sequences are assumed to be nucleotides, and will

  be 6-frame translated if the <queries> are amino acids. The output fields are

  [identity, shared-hashes, median-multiplicity, p-value, query-ID,

  query-comment], where median-multiplicity is computed for shared hashes, based

  on the number of observations of those hashes within the pool.

 

Options:

 

  Option    Description (range) [default]

 

  -h        Help

 

  -p <int>  Parallelism. This many threads will be spawned for processing. [1]

 

  -w        Winner-takes-all strategy for identity estimates. After counting

            hashes for each query, hashes that appear in multiple queries will

            be removed from all except the one with the best identity (ties

            broken by larger query), and other identities will be reduced. This

            removes output redundancy, providing a rough compositional outline.

 

...Output...

 

  -i <num>  Minimum identity to report. Inclusive unless set to zero, in which

            case only identities greater than zero (i.e. with at least one

            shared hash) will be reported. Set to -1 to output everything.

            (-1-1) [0]

 

  -v <num>  Maximum p-value to report. (0-1) [1.0]

 

mash info

 $ mash info

 

Version: 2.0

 

Usage:

 

  mash info [options] <sketch>

 

Description:

 

  Display information about sketch files.

 

Options:

 

  Option  Description (range) [default]

 

  -h      Help

 

  -H      Only show header info. Do not list each sketch. Incompatible with -d,

          -t and -c.

 

  -t      Tabular output (rather than padded), with no header. Incompatible with

          -d, -H and -c.

 

  -c      Show hash count histograms for each sketch. Incompatible with -d, -H

          and -t.

 

  -d      Dump sketches in JSON format. Incompatible with -H, -t, and -c.

 

ラン

テストデータをダウンロードする。

wget https://gembox.cbcb.umd.edu/mash/genome1.fna 
wget https://gembox.cbcb.umd.edu/mash/genome2.fna
wget https://gembox.cbcb.umd.edu/mash/genome3.fna

 

例1

disコマンドを走らせればすぐに2配列を比較できる。

mash dist genome1.fna genome2.fna

解析が終わると、STDOUTにReference-ID、Query-ID、Mash-distance、P-value、Matching-hashesが出力される。

genome1.fna  genome2.fna  0.0222766  0  456/1000

 

例2

結果は例1と同じだが、最初にsketchを走らせてsketcchファイルを作ることで、繰り返し時のランタイムを短縮する。

mash sketch genome1.fna 
mash sketch genome2.fna

sketchコマンドを実行するとgenome1.fna.msh、genome2.fna.mshというk-mer indexのバイナリファイルができる。中身はinfoコマンドで確認できる。

mash info genome1.fna.msh

 mash info genome1.fna.msh 

Header:

  Hash function (seed):          MurmurHash3_x64_128 (42)

  K-mer size:                    21 (64-bit hashes)

  Alphabet:                      ACGT (canonical)

  Target min-hashes per sketch:  1000

  Sketches:                      1

 

Sketches:

  [Hashes]  [Length]  [ID]         [Comment]

 

  1000      4639675   genome1.fna  -

 

distコマンドの引数に.mshファイルを指定してランする。 

mash dist genome1.fna.msh genome2.fna.msh

 

例3

genome3をgenome1、2と比較する。はじめにgenome1とgenome2のsketchを実行し、できた.mshファイルをgenome3と比較する。

mash sketch -o reference genome1.fna genome2.fna
mash info reference.msh
mash dist reference.msh genome3.fna

genome1.fna genome3.fna 0 0 1000/1000

genome2.fna genome3.fna 0.0222766 0 456/1000

 

例4

例1のように、distコマンドで1-3を比較するとgenome1とgenome2、genome1とgenome3の比較になる(例3との違いを確認する)。

mash dist genome1.fna genome2.fna genome3.fna

genome1.fna genome2.fna 0.0222766 0 456/1000

genome1.fna genome3.fna 0 0 1000/1000

 

例5

RefSeqのデータベースと比較する。著者があらかじめ.mshファイルを準備してくれているので、これをダウンロードする。

#chromosome
wget https://gembox.cbcb.umd.edu/mash/refseq.genomes.k21s1000.msh
#plasmid
wget https://gembox.cbcb.umd.edu/mash/refseq.plasmid.k21s1000.msh

 

 ペアエンドシーケンスデータはコンカテネートしておく。gzに圧縮してもO.K。

cat pair*fq > merge.fq

NGSのデータを使う時は、低頻度k-merを除外するオプションを指定する。-m 2なら頻度2以上のk-merのみ利用される。

mash sketch -m 2 -p 4 merge.fq
  • -m    Minimum copies of each k-mer required to pass noise filter for reads. Implies -r. [1]
  • -p     Parallelism. This many threads will be spawned for processing. [1]

 Refseqのデータベースを指定してdistコマンドを走らせる。

mash dist refseq.genomes.k21.s1000.msh merge.fq.msh > distances.tab

sortしてtop10ヒットを確認する。

sort -gk3 distances.tab | head

 

例6

コンタミが予想される時は、screenコマンドを使う。sketchでlow k-merを捨てず、fastqをそのままscreenの引数に指定する。

mash screen refseq.genomes.k21s1000.msh ERR024951.fastq > screen.tab
sort -gr screen.tab | head

トップヒットに特定の種の様々な株ばかりヒットしてくる場合、winner take all戦略のフラグ(-w)を立てることで回避できる。

mash screen -w -p 4 refseq.genomes.k21s1000.msh ERR024951.fastq > screen.tab 
sort -gr screen.tab | head
  • -w   Winner-takes-all strategy for identity estimates. After counting hashes for each query, hashes that appear in multiple queries will be removed from all except the one with the best identity (ties broken by larger query), and other identities will be reduced. This removes output redundancy, providing a rough compositional outline.  

 

例7

例えばメタゲノムのアセンブリなどの multi-fastaを比較する。はじめにfastaを分割する(FASTAの分割)。

#ヘッダー数を調べる
grep
-n ">" metagenome_assembly.fasta |wc -l
3621と出たとする。

#ヘッダ数で分割。
partition.sh in=metagenome_assembly.fasta out=contig%.fasta ways=3621

MASHで計算する。論文中(MethodsのMetagenomic clustering)では以下のパラメータが使用されている。 ワイルドカードを使ってFASTAを指定している。

mash sketch -s 400 -k 16 -o out contig*.fasta
  • -s     Sketch size. Each sketch will have at most this many non-redundant              min-hashes. [1000]
  • -k     K-mer size. Hashes will be based on strings of this many              nucleotides. Canonical nucleotides are used by default (see              Alphabet options below). (1-32) [21]
  • -o     Output prefix (first input file used if unspecified). The suffix '.msh' will be appended.

 out.mshファイルが1つできるので、これを全FASTAファイルと比較する。

 mash dist -p 10 out.msh *fasta > result

1行ずつ 全結果が出力されるので、行列に変換して例えばヒートマップ等で可視化する。

 

MASHのパラメータ例

kmer-dbのsupplementary data 

引用

Mash: fast genome and metagenome distance estimation using MinHash

Brian D. Ondov, Todd J. Treangen, Páll Melsted, Adam B. Mallonee, Nicholas H. Bergman, Sergey Koren, and Adam M. Phillippy

Genome Biol. 2016; 17: 132.