2022/04/17 コマンド修正
MinHash Sketch（BBSketchの解説）を構築し、Jaccard指数で比較・検索するsourmashは、発表後もバージョンアップが続けられていて、現在では様々なコマンドが利用できるようになっています。そこで今日は、sourmashのグネチャファイルをDNAやタンパク質から作るsketchコマンドと、低メモリで検索可能なシグネチャデータベースを作るindexコマンドなどを紹介します。
document
https://sourmash.readthedocs.io/en/latest/tutorials.html
インストール
ubuntu18でmambaを使ってインストールした。
本体 Github
mamba install -c bioconda sourmash -y
実行方法
DNA配列（ゲノム配列など）を使う。塩基配列の場合、K値は31-merに設定されている。
sourmash sketch dna *fasta
ゲノムごとにシグネチャファイルが出力され、-o output.sigのように出力名を指定すると、単一のシグネチャファイルが出力される（どちらでもsourmash compareに使える）。
compareコマンドで比較し、plotコマンドでヒートマップを出力する。
#20スレッド指定。k値はシグネチャ作成時と同じ。
sourmash compare *.sig -o output -k 31 --csv matrix.csv -p 20
sourmash plot --pdf --labels output
タンパク質配列（protein.faa）を使う。タンパク質の場合、K値は10-mer長のアミノ酸に設定されている。
sourmash sketch protein *faa
sourmash compare *.sig -o output -k 10 --csv matrix.csv --protein
sourmash plot --pdf --labels output
シグネチャファイル (.sig) からSBTインデックス；低メモリで高速に検索・収集可能なインデックスを作成し、類似性の高い配列を高速検索することができる。
sourmash index database *sig
- --from-file a text file containing a list of files to load signatures from
database.sbt.zipが出力される。
SBTインデックスから高い一致度を示すシグネチャファイル (.sig) を探索する。
sourmash search query.fasta.sig database.sbt.zip
SBTインデックスは大規模なゲノムデータベース（ntなど）から問い合わせするために特に有用で、大規模な例としてCOBS indexと共に作られた661K_sourmash_index_scaled.sbt.zip（ファイルサイズ48GB）が有名（論文）。BBSketchも同じような方法でネットワーク越しに高速検索しているものと思われる。sourmash のより実践的な使い方は、引用したsourmashのプロトコルの論文が役立ちます。
ダウンロードした661K_sourmash_index_scaled.sbt.zip（MinHashシグネチャ； num signatures）を調べてみる。
sourmash signature fileinfo 661K_sourmash_index_scaled.sbt.zip
661K_sourmash_index_scaled.sbt.zipはMinHashシグネチャのデータベースなので、クエリのほうもパラメータ固定のsourmash sketchコマンドではなく、パラメータを指定可能なsourmash computeの方でシグネチャファイルを作成する必要がある（例；sourmash compute -k 31 -n 5000 --dna genome*.fasta）。
引用
Large-scale sequence comparisons with sourmash.
Pierce NT, Irber L, Reiter T, Brooks P, Brown CT
F1000Res. 2019 Jul 4;8:1006
関連
