新たに塩基配列が決定された生物の多様性は極めて高く、最新の配列データベースは非常に大規模であるため、配列アノテーションにおける感度とスピードという相反するニーズの間で緊張関係が生じている。プロファイル隠れマルコフモデル(pHMM)に基づくアライメントは最先端の感度を示しているが、最近のアルゴリズムの進歩により、BLASTに近い感度を持つ超高速アノテーションツールが開発されている。
pHMMsが提供する感度の大部分を維持しながら、ここで紹介するnailと呼ばれるツールは、確率質量の大部分を含むFB動的計画行列のセルの疎な部分集合を特定することにより、pHMM Forward/Backward(FB)アルゴリズムの発見的近似を実行し、MMseqs2のような高速アノテーション手法に匹敵する速度に達する。本手法は、高速かつ少量のメモリでpHMMスコアとE値の正確な近似を生成する。タンパク質のベンチマークにおいて、nailはMMseqs2とHMMER間のリコール差の大部分を回復し、実行時間はHMMER3より25倍速い(MMseqs2の高感度バリアントより2.5倍遅いだけ)。nailはオープンBSD-3-clauseライセンスでリリースされ、https://github.com/TravisWheelerLab/nailから利用できる。
インストール
依存
- To run the nail search pipeline, mmseqs search and hmmbuild must be available in your system path.
https://github.com/TravisWheelerLab/nail
#mmseqs2
mamba install -c conda-forge -c bioconda mmseqs2 -y
#Hmmer suite
mamba install -c bioconda -y hmmer
#本体
git clone https://github.com/TravisWheelerLab/nail
cd nail/
cargo build --release
export PATH=$(pwd)/target/release//:$PATH
> target/release/nail -h
Using MMseqs2 to find rough alignment seeds, perform bounded profile HMM sequence alignment
Usage: nail <COMMAND>
Commands:
search Run the entire nail pipeline: prep, seed, & align
prep Prepare a query (MSA) file and target (fasta) file for the seed step
seed Use MMseqs2 to create a set of alignment seeds for the align step
align Search with the query against the target, using alignment seeds
help Print this message or the help of the given subcommand(s)
Options:
-h, --help Print help
> nail search -h
Run the entire nail pipeline: prep, seed, & align
Usage: nail search [OPTIONS] <QUERY.[fasta:sto]> <TARGET.fasta>
Arguments:
<QUERY.[fasta:sto]> Query file
<TARGET.fasta> Target file
Options:
-q, --query-hmm <QUERY.hmm>
The path to a pre-built P7HMM file
-p, --prep <PATH>
Where MMseqs2 intermediate files are placed [default: ./prep/]
-E <F>
Only report hits with an E-value below this value [default: 10]
-T, --tab-output <path>
Where to place tabular output [default: results.tsv]
-O, --output <path>
Where to place alignment output
-Z <N>
Override the target database size (number of sequences) used for E-value calculation
-A <F>
Pruning parameter alpha [default: 12]
-B <F>
Pruning parameter beta [default: 20]
-G <N>
Pruning parameter gamma [default: 5]
--cloud-thresh <F>
The P-value threshold for promoting hits past cloud search [default: 0.001]
--forward-thresh <F>
The P-value threshold for promoting hits past forward [default: 0.0001]
--full-dp
Compute the full dynamic programming matrices during alignment
--mmseqs-k <K>
MMseqs2 prefilter: k-mer length (0: automatically set to optimum) [default: 0]
--mmseqs-k-score <K_SCORE>
MMseqs2 prefilter: k-mer threshold for generating similar k-mer lists [default: 80]
--mmseqs-min-ungapped_score <MIN_UNGAPPED_SCORE>
MMseqs2 prefilter: Accept only matches with ungapped alignment score above threshold [default: 15]
--mmseqs-max-seqs <MAX_SEQS>
MMseqs2 prefilter: Maximum results per query sequence allowed to pass the prefilter [default: 1000]
--mmseqs-pvalue-threshold <PVALUE_THRESHOLD>
MMseqs2 align: Include matches below this P-value as seeds [default: 0.01]
-t, --threads <n>
The number of threads to use [default: 8]
-h, --help
Print help (see more with '--help')
> nail prep -h
Prepare a query (MSA) file and target (fasta) file for the seed step
Usage: nail prep [OPTIONS] <QUERY.[fasta:sto]> <TARGET.fasta>
Arguments:
<QUERY.[fasta:sto]> Query file
<TARGET.fasta> Target file
Options:
--skip-hmmbuild Don't build a profile HMM with the input MSA
-p, --prep <PATH> Where MMseqs2 intermediate files are placed [default: ./prep/]
-t, --threads <n> The number of threads to use [default: 8]
-h, --help Print help
> nail seed -h
Use MMseqs2 to create a set of alignment seeds for the align step
Usage: nail seed [OPTIONS] <PATH>
Arguments:
<PATH> The location of files prepared with nail prep
Options:
-s, --seeds-path <SEEDS_PATH>
Where to place the seeds output file [default: seeds.json]
-q, --query-hmm <QUERY.hmm>
The path to a pre-built P7HMM file
--mmseqs-k <K>
MMseqs2 prefilter: k-mer length (0: automatically set to optimum) [default: 0]
--mmseqs-k-score <K_SCORE>
MMseqs2 prefilter: k-mer threshold for generating similar k-mer lists [default: 80]
--mmseqs-min-ungapped_score <MIN_UNGAPPED_SCORE>
MMseqs2 prefilter: Accept only matches with ungapped alignment score above threshold [default: 15]
--mmseqs-max-seqs <MAX_SEQS>
MMseqs2 prefilter: Maximum results per query sequence allowed to pass the prefilter [default: 1000]
--mmseqs-pvalue-threshold <PVALUE_THRESHOLD>
MMseqs2 align: Include matches below this P-value as seeds [default: 0.01]
-t, --threads <n>
The number of threads to use [default: 8]
-h, --help
Print help (see more with '--help')
> nail align -h
Search with the query against the target, using alignment seeds
Usage: nail align [OPTIONS] <QUERY.[fasta:hmm:sto]> <TARGET.fasta> <SEEDS.json>
Arguments:
<QUERY.[fasta:hmm:sto]> Query file
<TARGET.fasta> Target file
<SEEDS.json> Alignment seeds from running nail seed (or elsewhere)
Options:
-Z <N> Override the target database size (number of sequences) used for E-value calculation
-A <F> Pruning parameter alpha [default: 12]
-B <F> Pruning parameter beta [default: 20]
-G <N> Pruning parameter gamma [default: 5]
--cloud-thresh <F> The P-value threshold for promoting hits past cloud search [default: 0.001]
--forward-thresh <F> The P-value threshold for promoting hits past forward [default: 0.0001]
--full-dp Compute the full dynamic programming matrices during alignment
-E <F> Only report hits with an E-value below this value [default: 10]
-T, --tab-output <path> Where to place tabular output [default: results.tsv]
-O, --output <path> Where to place alignment output
-t, --threads <n> The number of threads to use [default: 8]
-h, --help Print help
実行方法
4つのサブコマンド:search、prep、seed、align、がある。nail searchは、prep、seed、alignの各ステップを連続して実行する。
searchサブコマンドの入力は、Stockholm形式のクエリMSA(Multiple Sequence Alignment)ファイルと、FASTA形式のターゲット配列データベースファイル。順番に指定する(注:Stockholm形式でなくシングルのfastaファイルも認識自体はする)。
#nail searchのラン
nail search query.sto target.fa > align.txt
出力
カレントにprepディレクトリができる。最終出力は表形式でresults.tsv に書き出され、アライメント出力は標準出力に書き出される。
results.tsvはタブ区切りで、
- ターゲット名
- クエリ名
- ターゲット開始
- ターゲット終了
- クエリ開始
- クエリ終了
- ビットスコア
- コンポジションバイアスコア
- E値
- セル率(スパースF/Bで計算されたセルの割合)
がプリントされている。
レポジトリより
-
各パイプラインステップを分離して実行することもサポートしている。これは特に以下のような場合に便利である:
prepまたはseedの結果をキャッシュし、次のステップでパラメータを実験したい場合。
MMseqs2以外のソースから生成されたアライメントシードで実験したい場合。
引用
nail: software for high-speed, high-sensitivity protein sequence annotation
Jack W. Roddy, David H. Rich, Travis J. Wheeler
bioRxiv, Posted January 30, 2024.
関連
https://kazumaxneo.hatenablog.com/entry/2018/09/22/220752