2019 4/12 dockerリンク追加
2021 3/25 condaインストール追記
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/releases
#bioconda(link)
mamba install -c bioconda mash -y
> mash
$ mash
Mash version 2.3
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 mixture of
sequences.
sketch Create sketches (reduced representations for fast operations).
taxscreen Create Kraken-style taxonomic report based on mash screen.
triangle Estimate a lower-triangular distance matrix.
> mash sketch
$ mash sketch
Version: 2.3
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...
-I <path> ID field for sketch of reads (instead of first sequence ID).
-C <path> Comment for a sketch of reads (instead of first sequence comment).
-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.3
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]
-C Show comment fields with reference/query names (denoted with ':').
(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.3
Usage:
mash screen [options] <queries>.msh <mixture> [<mixture>] ...
Description:
Determine how well query sequences are contained within a mixture of
sequences. The queries must be formatted as a single Mash sketch file (.msh),
created with the `mash sketch` command. The <mixture> files can be contigs or
reads, in fasta or fastq, gzipped or not, and "-" can be given for <mixture>
to read from standard input. The <mixture> 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
mixture.
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.3
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.
追記
dockerイメージ
https://hub.docker.com/r/staphb/mash/
ラン
テストデータをダウンロードする。
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
例１
disコマンドを走らせればすぐに２配列を比較できる。
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
例２
結果は例１と同じだが、最初に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
例３
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
例４
例１のように、distコマンドで1-3を比較するとgenome1とgenome2、genome1とgenome3の比較になる（例３との違いを確認する）。
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
例５
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
例６
コンタミが予想される時は、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.
例７
例えばメタゲノムのアセンブリなどの 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ファイルが１つできるので、これを全FASTAファイルと比較する。
mash dist -p 10 out.msh *fasta > result
１行ずつ 全結果が出力されるので、行列に変換して例えばヒートマップ等で可視化する。
MASHのパラメータ例
引用
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.
