2019 2/17 テスト環境の誤り修正
Mashツール[ref.1]のリリース以来、MinHashのようなデータスケッチは比較ゲノミクスにおいて有益になっている。それらは大規模データベースからのゲノムのクラスター化[ref.1]、特定のシーケンス内容を持つデータセットの検索[ref.2]、ゲノムアセンブラでの重複ステップの加速[ref.3、4]、シーケンスリードのマッピング[ref.5]、類似性しきい値の検索に使用される。また種レベルの違いを特徴付けるためにも使用される[ref.6]。 MinHashはもともと似たようなWebページを見つけるために開発されたものだが[ref.7]、ここではリファレンスゲノムやシーケンスデータセットのような大きなゲノム配列コレクションをまとめるために使われている。コレクションは代表的なk-merの集合に縮小され、最終的に整数のリストとして格納される。要約は元のデータよりはるかに小さいが、結合のサイズや2つのゲノムのk-merコンテンツの共通部分など、関連する集合濃度を推定するために使用できる。これらの濃度から、Jaccard係数（J）または「MASH distance」を得ることができる。これは、平均ヌクレオチド同一性（ANI）の代用に利用できる[ref.1]。これらは、配列をクラスター化すること、およびそうでなければ巨大なゲノム最近傍問題を解決することを可能にする。
MinHashはバイオインフォマティクスの他のコアメソッドに関連している。MinHashの特別な場合と考えることができるミニマイザーは、メタゲノミクス分類[ref.8]とアライメントとアセンブリ[ref.9]で広く使われている。より一般的には、MinHashは一種の局所性依存ハッシュ（LSH）と見なすことができる。これは、類似の入力を同じ値にマップするように設計されたハッシュ関数を含む。LSHは、相同性検索[ref.10]やメタゲノム分類[ref.11]など、バイオインフォマティクスにも使用されている。
MinHashの効用に拍車をかけて、他のグループは検索とデータマイニングからの新しいアイデアを使った代替案を提案した。 BinDash [ref.12]は、小さいメモリフットプリントでMashと比較して高い精度と速度を達成するために、bビットの1順列ローリングMinHashを使用する。HyperMin-Hash [ref.13]とSuperMinHash [ref.14]の研究では、他の理論的な改善が提案されている。
いくつかの研究はMinHashの欠点を指摘した。KoslickiとZabetiは、集合が非常に異なるサイズであると、MinHash濃度推定値が損なわれると主張している[ref.15]。これは珍しいシナリオではない。非常に異なる長さの2つのゲノム間の距離を見つけるとき、または短い配列（たとえば細菌ゲノム）と大規模なコレクション（たとえば、広範囲にわたるメタゲノム解析データセット）の間の類似性を見つけるとき。
ここでは、入力セットのサイズが非常に異なる場合やスケッチデータ構造が非常に小さい場合など、さまざまなシナリオで優れた正確性と速度を示すMinHashの代替としてHyperLogLog（HLL）スケッチ[ref.16]を使用する。 HLLはバイオインフォマティクスの他の分野、例えば、ゲノムまたはデータコレクション中の異なるk-mer数を数えること［ref.17、18、19］。集合和集合と交点[ref.20]、Jを推定するのに必要な要素、その他の類似度を評価するための基数推定における最近の理論的改良をさらに使用する。
GPLv3ライセンスの下でフリーでオープンソースのHashをDashingソフトウェアツール（https://github.com/ dnbaker / dashing）に実装した。 DashingはMash [ref.1]（紹介）、BinDash [ref.12]（紹介）、Sourmash [ref.21]（紹介）のような同様のツールで利用可能な機能をサポートしている。ダッシュは、FASTAファイル（アセンブリされたゲノム用）またはFASTQファイル（データセットのシーケンス用）を含む、入力シーケンスセットのスケッチ（ダッシュスケッチ）を作成することができる。 Dashingには、スケッチ前にシーケンスエラーを含む可能性があるk-merを削除するためのスケッチベースの機能がある。dashing dist関数は、例えばRefSeqデータベースからのすべての完全なゲノムなど、大規模なコレクションでデータセットの全ペア間の距離比較を実行する。 Dashingのスケッチ機能は非常に高速であるため、Dashingはスケッチと全ペア距離計算の両方を同じコマンドで実行でき、ステップ間でスケッチをディスクに保存する必要がない。Dashingは並列化されており、効率的に100スレッドに拡張できることを示している。Dashingはまた、現代の汎用コンピュータプロセッサ上でSIMD（Single Instruction Multiple Data）命令を使用して、HLL計算に固有のよりきめの細かい並列処理を利用する。
Dashingに関するツイート
Revised preprint by @dnb_hopkins on Dashing: Fast and Accurate Genomic Distances with HyperLogLog. https://t.co/bf2nttTnEr. v0.1.2 has portability improvements and bug fixes. Changes from 1st preprint are: (a) new Methods figure explaining HLL & cardinality estimation... (1/3)
— Ben Langmead (@BenLangmead) February 5, 2019
インストール
ubuntu18.04でテストした。
依存
- Dashing is written in C++14, which means that it requires a relatively new compiler. Dashing is tested under gcc
{5.4-9}, but fails for gcc4. For OSX, we recommend using Homebrew to install gcc-8（"brew install gcc@8"）. On Linux, we recommend package managers. (For instance, our Travis-CI Ubuntu example upgrades to a sufficiently new GCC usingsudo update-alternatives.
本体 Github
git clone https://github.com/dnbaker/dashing
cd dashing && make update dashing
> dashing
# ./dashing --help
Usage: ./dashing <subcommand> [options...]. Use ./dashing <subcommand> for more options. [Subcommands: sketch, dist, setdist, hll, printmat.]
> dashing dist
# ./dashing dist
Usage: dist <opts> [genomes if not provided from a file with -F]
Flags:
-h/-? Usage
-k Set kmer size [31]
-W Cache sketches/use cached sketches
-p Set number of threads [1]
-b Emit distances in binary (default: human-readable, upper-triangular)
-U Emit distances in PHYLIP upper triangular format(default: human-readable, upper-triangular)
-s add a spacer of the format <int>x<int>,<int>x<int>,..., where the first integer corresponds to the space between bases repeated the second integer number of times
-w Set window size [max(size of spaced kmer, [parameter])]
-S Set sketch size [10, for 2**10 bytes each]
-H Treat provided paths as pre-made sketches.
-C Do not canonicalize. [Default: canonicalize]
-P Set prefix for sketch file locations [empty]
-x Set suffix in sketch file names [empty]
-o Output for genome size estimates [stdout]
-I Use Ertl's Improved Estimator
-E Use Ertl's Original Estimator
-J Use Ertl's JMLE Estimator [default:Uses Ertl-MLE]
-O Output for genome distance matrix [stdout]
-e Emit in scientific notation
-F Get paths to genomes from file rather than positional arguments
-M Emit Mash distance (default: jaccard index)
-T postprocess binary format to human-readable TSV (not upper triangular)
-Z Emit genome sizes (default: jaccard index)
-N Autodetect fastq or fasta data by filename (.fq or .fastq within filename).
-y Filter all input data by count-min sketch.
-q Set count-min number of hashes. Default: [4]
-c Set minimum count for kmers to pass count-min filtering.
-t Set count-min sketch size (log2). Default: ceil(log2(max_filesize)) + 2
-R Set seed for seeds for count-min sketches
> dashing sketch
# dashing sketch
No paths. See usage.
Usage: sketch <opts> [genomes if not provided from a file with -F]
Flags:
-h/-?: Emit usage
Sketch options --
-k Set kmer size [31]
-s add a spacer of the format <int>x<int>,<int>x<int>,..., where the first integer corresponds to the space between bases repeated the second integer number of times
-w Set window size [max(size of spaced kmer, [parameter])]
-S Set sketch size [10, for 2**10 bytes each]
-C Do not canonicalize. [Default: canonicalize]
Run options --
-p Set number of threads [1]
-P Set prefix for sketch file locations [empty]
-x Set suffix in sketch file names [empty]
-F Get paths to genomes from file rather than positional arguments
-c Cache sketches/use cached sketches (save sketches to disk in directory of the file [default] or in folder specified by -P
-z Write gzip compressed. (Or zstd-compressed, if compiled with zlibWrapper.
Estimation methods --
-E Use Flajolet with inclusion/exclusion quantitation method for hll. [Default: Ertl MLE]
-I Use Ertl Improved estimator [Default: Ertl MLE]
-J Use Ertl JMLE
Filtering Options --
Default: consume all kmers. Alternate options:
-f Autodetect fastq or fasta data by filename (.fq or .fastq within filename).
-B Filter all input data by count-min sketch.
Options for count-min filtering --
-H Set count-min number of hashes. Default: [4]
-q Set count-min sketch size (log2). Default: ceil(log2(max_filesize)) + 2
-n Provide minimum expected count for fastq data. If unspecified, all kmers are passed.
-R Set seed for seeds for count-min sketches
----
> dashing setdist
# dashing setdist
No paths. See usage.
Usage: setdist <opts> [genomes if not provided from a file with -F]
Flags:
-h/-? Usage
-k Set kmer size [31]
-W Cache sketches/use cached sketches
-p Set number of threads [1]
-b Emit distances in binary (default: human-readable, upper-triangular)
-U Emit distances in PHYLIP upper triangular format(default: human-readable, upper-triangular)
-s add a spacer of the format <int>x<int>,<int>x<int>,..., where the first integer corresponds to the space between bases repeated the second integer number of times
-w Set window size [max(size of spaced kmer, [parameter])]
-S Set sketch size [10, for 2**10 bytes each]
-H Treat provided paths as pre-made sketches.
-C Do not canonicalize. [Default: canonicalize]
-P Set prefix for sketch file locations [empty]
-x Set suffix in sketch file names [empty]
-o Output for genome size estimates [stdout]
-I Use Ertl's Improved Estimator
-E Use Ertl's Original Estimator
-J Use Ertl's JMLE Estimator [default:Uses Ertl-MLE]
-O Output for genome distance matrix [stdout]
-e Emit in scientific notation
-F Get paths to genomes from file rather than positional arguments
-M Emit Mash distance (default: jaccard index)
-T postprocess binary format to human-readable TSV (not upper triangular)
-Z Emit genome sizes (default: jaccard index)
-N Autodetect fastq or fasta data by filename (.fq or .fastq within filename).
-y Filter all input data by count-min sketch.
-q Set count-min number of hashes. Default: [4]
-c Set minimum count for kmers to pass count-min filtering.
-t Set count-min sketch size (log2). Default: ceil(log2(max_filesize)) + 2
-R Set seed for seeds for count-min sketches
> dashing hll
# dashing hll
[E:int bns::hll_main(int, char**):872] Usage: hll <opts> <paths>
Flags:
-k: kmer length (Default: 31. Max: 31)
-w: window size (Default: -1) Must be -1 (ignored) or >= kmer length.
-s: spacing (default: none). format: <value>x<times>,<value>x<times>,...
Omitting x<times> indicates 1 occurrence of spacing <value>
-S: sketch size (default: 24). (Allocates 2 << [param] bytes of memory per HyperLogLog.
-p: number of threads.
-F: Path to file which contains one path per line
> dashing printmat
# dashing printmat
printmat printmat <path to binary file> [- to read from stdin]
-o Specify output file (default: stdout)
-s Emit in scientific notation
実行方法
distコマンド
#genomeのfastaを指定
dashing dist -k 31 -p 13 -O distance_matrix.txt -o size_estimates.txt genome1.fna.gz genome2.fna genome3.fasta
#またはワイルドカードで指定
dashing dist -k 31 -p 13 -O distance_matrix.txt -o size_estimates.txt genome*gz
#list（genome.fastaのfull path）で指定
dashing dist -k 31 -p 13 -O distance_matrix.txt -o size_estimates.txt -F list.txt
sketchコマンドはdistコマンドとほぼ同じだが、sketchのみ計算する。
テストラン
RefSeqのゲノムをダウンロードしてテストする。ncbi-genome-download（紹介）かpyani（紹介）のスクリプトを使うとtaxidsや分類学の系統を指定して自動でダウンロードできる。
#list作成
find genome_dir/*fa -type f > list
#listファイルを指定してdashing dist実行
dashing dist -k 31 -p 13 -O distance_matrix.txt -o size_estimates.txt -F list.txt
正当なANI計算のツール（*1）でAll vs AllのANI計算に4日ほどかかったおよそ1500バクテリアゲノムを使ってdistを実行してみると、3秒で終わった（*2）。ANI近似値だがおよそ10^5のオーダーも早く終わったことになる。
NCBI GenomeのBrowse by Organism（リンク）からgenomeレポートを表示して確認すると、bacteria genomeの総数は189502あったが（190216時点）、completeはこのうち10%以下だった。このcomplete genome全てを使っても、Dashingを使えば、数分〜数十分で計算が終わることになる。今後も増え続けるスモールゲノムの大規模比較にも対応できるスケール性能を持ってると言えますね（Preprint中で議論されてます）。
他のコマンドはGithubの説明を読んでください。
2021/11/12
Dashing2
Dashing: Fast and Accurate Genomic Distances with HyperLogLog
Daniel N Baker, Ben Langmead
bioRxiv preprint first posted online Dec. 20, 2018
*1
xeon E5 v4 2680 x2のマシンで56 threds使用。CPU使用使用率は４日間ほぼ100%に張り付いていた。
*2
xeon E5 v4 2680x2のマシンで56 threads指定時。16threadsだと5.8秒かかった。
