jackhmmer(1)

📦
version 3.4 · 3 versions
, in debian trixie, nixos 26.05, nixpkgs unstable, ubuntu noble, ubuntu resolute, ubuntu stonking

jackhmmer(1) — HMMER Manual 💖

okay so bestie, welcome to jackhmmer, the HMMER tool that just keeps going and going and going (affectionate). this is part 1 of the page, so buckle up, we're not done yet 🎀

✨ NAME ✨

jackhmmer — iteratively search sequence(s) against a sequence database 🦋

✨ SYNOPSIS ✨

jackhmmer [options] seqfile seqdb

✨ DESCRIPTION ✨

so basically, jackhmmer takes every query sequence in seqfile and searches it against the target sequence(s) in seqdb — but it doesn't just do it once, it iterates. like, it keeps searching, then rebuilding, then searching again. obsessed.

the first iteration is literally identical to a phmmer search. no cap.

then for the next iteration, here's the glow-up: a multiple alignment gets assembled out of the query plus all the target sequences that satisfied the inclusion thresholds ✨, a profile gets constructed from that alignment (which is identical to using hmmbuild on the alignment), and then a profile search of seqdb happens (identical to an hmmsearch with the profile). it's giving main character energy for real.

now, the query seqfile can be '-' (a dash character), and in that case the query sequences are read from a stdin pipe instead of from a file 🫶

as for seqdb — bestie, listen — it needs to be a 'normal' sequence file. it cannot be read from a stdin stream, because jackhmmer needs to do multiple passes over the database. it cannot be a compressed (gzipped) file either, because gzipped files are treated essentially as stdin streams, calling an external decompression program. so like, don't do that 💀

the output format is designed to be human-readable, but it's often so voluminous that reading it is impractical and parsing it is a whole pain 😭. that's why the --tblout and --domtblout options save output in simple tabular formats that are concise and easier to parse. and the -o option lets you redirect the main output, including throwing it away in /dev/null. iconic.

✨ OPTIONS ✨

  • -h — help! prints a brief reminder of command line usage and all available options 🌸
  • -N <n> — set the maximum number of iterations to <n>. the default is 5. and if N=1, the result is equivalent to a phmmer search ✨

✨ OPTIONS CONTROLLING OUTPUT ✨

by default, output for each iteration appears on stdout in a somewhat human readable, somewhat parseable format. these options let you redirect that output or save additional kinds of output to files, including checkpoint files for each iteration 🎀

  • -o <f> — direct the human-readable output to a file <f>.
  • -A <f> — after the final iteration, save an annotated multiple alignment of all hits satisfying inclusion thresholds (also including the original query) to <f> in Stockholm format 🦋
  • --tblout <f> — after the final iteration, save a tabular summary of top sequence hits to <f> in a readily parseable, columnar, whitespace-delimited format 📊
  • --domtblout <f> — after the final iteration, save a tabular summary of top domain hits to <f> in a readily parseable, columnar, whitespace-delimited format 📊
  • --chkhmm prefix — at the start of each iteration, checkpoint the query HMM, saving it to a file named prefix-n.hmm where n is the iteration number (from 1..N) 💾
  • --chkali prefix — at the end of each iteration, checkpoint an alignment of all domains satisfying inclusion thresholds (e.g. what will become the query HMM for the next iteration), saving it to a file named prefix-n.sto in Stockholm format, where n is the iteration number (from 1..N) 💾
  • --acc — use accessions instead of names in the main output, where available for profiles and/or sequences.
  • --noali — omit the alignment section from the main output. this can greatly reduce the output volume 🙈 (ngl, sometimes you need this)
  • --notextw — unlimit the length of each line in the main output. the default is a limit of 120 characters per line, which helps in displaying the output cleanly on terminals and in editors, but can truncate target profile description lines 😭
  • --textw <n> — set the main output's line length limit to <n> characters per line. the default is 120.

✨ OPTIONS CONTROLLING SINGLE SEQUENCE SCORING (FIRST ITERATION) ✨

by default, the first iteration uses a search model constructed from a single query sequence. this model is constructed using a standard 20x20 substitution matrix for residue probabilities, plus two additional parameters for position-independent gap open and gap extend probabilities. these options let you change the default single-sequence scoring parameters 🧚

  • --popen <x> — set the gap open probability for a single sequence query model to <x>. the default is 0.02. <x> must be >= 0 and < 0.5.
  • --pextend <x> — set the gap extend probability for a single sequence query model to <x>. the default is 0.4. <x> must be >= 0 and < 1.0.
  • --mx <s> — obtain residue alignment probabilities from the built-in substitution matrix named <s>. several standard matrices are built-in and do not need to be read from files 🔥 the matrix name <s> can be PAM30, PAM70, PAM120, PAM240, BLOSUM45, BLOSUM50, BLOSUM62, BLOSUM80, or BLOSUM90. only one of the --mx and --mxfile options may be used, bestie, pick a lane 💅
  • --mxfile mxfile — obtain residue alignment probabilities from the substitution matrix in file mxfile. the default score matrix is BLOSUM62 (this matrix is internal to HMMER and does not have to be available as a file). the format of a substitution matrix mxfile is the standard format accepted by BLAST, FASTA, and other sequence analysis software. see ftp.ncbi.nlm.nih.gov/blast/matrices/ for example files. (the only exception: matrices are required to be square, so for DNA, use files like NCBI's NUC.4.4, not NUC.4.2.)

✨ OPTIONS CONTROLLING REPORTING THRESHOLDS ✨

reporting thresholds control which hits are reported in output files (the main output, --tblout, and --domtblout). in each iteration, sequence hits and domain hits are ranked by statistical significance (E-value) and output is generated in two sections called per-target and per-domain output. in per-target output, by default, all sequence hits with an E-value <= 10 are reported. in the per-domain output, for each target that has passed per-target reporting thresholds, all domains satisfying per-domain reporting thresholds are reported. by default, these are domains with conditional E-values of <= 10. the following options let you change the default E-value reporting thresholds, or use bit score thresholds instead 🌈

  • -E <x> — report sequences with E-values <= <x> in per-sequence output. the default is 10.0.
  • -T <x> — use a bit score threshold for per-sequence output instead of an E-value threshold (any setting of -E is ignored). report sequences with a bit score of >= <x>. by default this option is unset.
  • -Z <x> — declare the total size of the database to be <x> sequences, for purposes of E-value calculation. normally E-values are calculated relative to the size of the database you actually searched (e.g. the number of sequences in target_seqdb). in some cases (for instance, if you've split your target sequence database into multiple files for parallelization of your search), you may know better what the actual size of your search space is 🫶
  • --domE <x> — report domains with conditional E-values <= <x> in per-domain output, in addition to the top-scoring domain per significant sequence hit. the default is 10.0.
  • --domT <x> — use a bit score threshold for per-domain output instead of an E-value threshold (any setting of --domE is ignored). report domains with a bit score of >= <x> in per-domain output, in addition to the top-scoring domain per significant sequence hit. by default this option is unset.
  • --domZ <x> — declare the number of significant sequences to be <x> sequences, for purposes of conditional E-value calculation for additional domain significance. normally conditional E-values are calculated relative to the number of sequences passing per-sequence reporting threshold.

✨ OPTIONS CONTROLLING INCLUSION THRESHOLDS ✨

inclusion thresholds control which hits are included in the multiple alignment and profile constructed for the next search iteration 👑 by default, a sequence must have a per-sequence E-value of <= 0.001 (see -E option) to be included, and any additional domains in it besides the top-scoring one must have a conditional E-value of <= 0.001 (see --domE option). the difference between reporting thresholds and inclusion thresholds is that inclusion thresholds control which hits actually get used in the next iteration (or the final output multiple alignment if the -A option is used), whereas reporting thresholds control what you see in output. reporting thresholds are generally more loose so you can see borderline hits in the top of the noise that might be of interest 👀

  • --incE <x> — include sequences with E-values <= <x> in subsequent iteration or final alignment output by -A. the default is 0.001.
  • --incT <x> — use a bit score threshold for per-sequence inclusion instead of an E-value threshold (any setting of --incE is ignored). include sequences with a bit score of >= <x>. by default this option is unset.
  • --incdomE <x> — include domains with conditional E-values <= <x> in subsequent iteration or final alignment output by -A, in addition to the top-scoring domain per significant sequence hit. the default is 0.001.
  • --incdomT <x> — use a bit score threshold for per-domain inclusion instead of an E-value threshold (any setting of --incdomE is ignored). include domains with a bit score of >= <x>. by default this option is unset.

okay that's part 1 bestie, and there's more coming so don't go anywhere 💖

🔥 OPTIONS CONTROLLING ACCELERATION HEURISTICS 🔥

okay so here's the tea about how HMMER3 moves so fast: it runs a three-step filter pipeline, like a lil assembly line of vibes. We got the MSV filter, then the Viterbi filter, then the Forward filter. First one is the speediest and the most "eh, close enough" — the last one is the full Forward scoring algorithm, which is the slowest but the most accurate. There's also a bias filter step chilling between MSV and Viterbi. Any target that survives all the steps of the acceleration pipeline then gets sent to postprocessing, which is where we do domain identification and scoring with the Forward/Backward algorithm. 💅

Basically the only free parameters you get to play with in HMMER's heuristic filters are the P-value thresholds — they control the expected fraction of nonhomologous sequences that make it through the filters. Crank the default thresholds higher and more nonhomologous sequence passes through, which ups your sensitivity but costs you speed. Flip it: lower P-value thresholds mean a smaller fraction passes, so you lose a lil sensitivity but gain speed. And setting a filter's P-value threshold to 1.0 means it'll pass all sequences — which effectively disables the filter entirely. periodt.

Lowkey important: changing filter thresholds only removes or includes targets from consideration. It does NOT change bit scores, E-values, or alignments — those are all decided solely in postprocessing. So don't come for the scores, bestie, they're not the vibe you're tweaking here. 🫶
  • --max — Maximum sensitivity mode. Turns off all the filters, bias filter included, and runs full Forward/Backward postprocessing on every single target. It bumps sensitivity up a lil bit, but the speed cost is huge. it's giving "I have time and a big computer" energy.
  • --F1 <x> — First filter threshold; sets the P-value threshold for the MSV filter step. The default is 0.02, which means roughly 2% of the highest scoring nonhomologous targets are expected to pass the filter. 💖
  • --F2 <x> — Second filter threshold; sets the P-value threshold for the Viterbi filter step. The default is 0.001.
  • --F3 <x> — Third filter threshold; sets the P-value threshold for the Forward filter step. The default is 1e-5.
  • --nobias — Turn off the bias filter. This increases sensitivity somewhat, but can come at a high cost in speed, especially if your query has biased residue composition — think a repetitive sequence region, or a membrane protein with big regions of hydrophobicity. Without the bias filter, way too many sequences may pass the filter when your query is biased, and then performance gets slower than expected because the computationally intensive Forward/Backward algorithms end up shouldering an abnormally heavy load. 😭 no cap, that's the whole trap.

🎀 OPTIONS CONTROLLING PROFILE CONSTRUCTION (LATER ITERATIONS) 🎀

okay so here's a thing that's just inherently true about jackhmmer: it always includes your original query sequence in the alignment result at every iteration, and consensus positions are always defined by that query sequence. Like, always. Which means a jackhmmer profile is always the same length as your original query, at every iteration. fr fr.

So that's why jackhmmer gives you less control over profile construction than hmmbuild does — it doesn't have the --fast, or --hand, or --symfrac options. The one and only profile construction option available in jackhmmer is --fragthresh, bestie:

  • --fragthresh <x> — We only wanna count terminal gaps as deletions if the aligned sequence is known to be full-length, not if it's a fragment (like, say, only part of it got sequenced). HMMER uses a simple lil rule to infer fragments: if the sequence length L is less than or equal to a fraction <x> times the alignment length in columns, then the sequence gets handled as a fragment. Default is 0.5. Setting --fragthresh 0 will define no (nonempty) sequence as a fragment — you might wanna do this if you know you've got a carefully curated alignment of full-length sequences. Setting --fragthresh 1 will define all sequences as fragments — you might wanna do this if you know your alignment is entirely composed of fragments, like translated short reads in metagenomic shotgun data. pick your fighter ✨

💜 OPTIONS CONTROLLING RELATIVE WEIGHTS 💜

Whenever a profile gets built from a multiple alignment, HMMER uses an ad hoc sequence weighting algorithm to downweight closely related sequences and upweight distantly related ones. The whole point is making models less biased by uneven phylogenetic representation. For example, two identical sequences would typically each receive half the weight that one sequence would — and that's exactly why jackhmmer isn't worried about always including your original query sequence in each iteration's alignment, even if it finds it again in the database you're searching. 😌 it's math, babe.

These options control which algorithm gets used:

  • --wpb — Use the Henikoff position-based sequence weighting scheme [Henikoff and Henikoff, J. Mol. Biol. 243:574, 1994]. This is the default. 👑
  • --wgsc — Use the Gerstein/Sonnhammer/Chothia weighting algorithm [Gerstein et al, J. Mol. Biol. 235:1067, 1994].
  • --wblosum — Use the same clustering scheme that was used to weight data in calculating BLOSUM substitution matrices [Henikoff and Henikoff, Proc. Natl. Acad. Sci 89:10915, 1992]. Sequences are single-linkage clustered at an identity threshold (default 0.62; see --wid) and within each cluster of c sequences, each sequence gets relative weight 1/c.
  • --wnone — No relative weights. All sequences get assigned uniform weight. everyone's equal, we love democracy. 🌈
  • --wid <x> — Sets the identity threshold used by single-linkage clustering when using --wblosum. Invalid with any other weighting scheme — don't try it. Default is 0.62.

🦋 OPTIONS CONTROLLING EFFECTIVE SEQUENCE NUMBER 🦋

After relative weights are determined, they get normalized to sum to a total effective sequence number, eff_nseq. This number may be the actual number of sequences in the alignment, but it's almost always smaller than that. The default entropy weighting method (--eent) reduces the effective sequence number to reduce the information content — that's relative entropy, or the average expected score on true homologs — per consensus position. The target relative entropy is controlled by a two-parameter function, where the two parameters are settable with --ere and --esigma. 🧚

  • --eent — Adjust effective sequence number to achieve a specific relative entropy per position (see --ere). This is the default.
  • --eclust — Set effective sequence number to the number of single-linkage clusters at a specific identity threshold (see --eid). This option is not recommended; it's for experiments evaluating how much better --eent is. it's giving "science experiment, not production."
  • --enone — Turn off effective sequence number determination and just use the actual number of sequences. One reason you might wanna do this is to try to maximize the relative entropy/position of your model, which may be useful for short models.
  • --eset <x> — Explicitly set the effective sequence number for all models to <x>.
  • --ere <x> — Set the minimum relative entropy/position target to <x>. Requires --eent. Default depends on the sequence alphabet; for protein sequences, it's 0.59 bits/position.
  • --esigma <x> — Sets the minimum relative entropy contributed by an entire model alignment, over its whole length. The effect here is that short models end up with higher relative entropy per position than --ere alone would give. Default is 45.0 bits.
  • --eid <x> — Sets the fractional pairwise identity cutoff used by single linkage clustering with the --eclust option. Default is 0.62.

⭐ OPTIONS CONTROLLING PRIORS ⭐

In profile construction, by default, weighted counts get converted to mean posterior probability parameter estimates using mixture Dirichlet priors. Default mixture Dirichlet prior parameters for protein models and for nucleic acid (RNA and DNA) models are built in. The following options let you override the default priors. ✨

  • --pnone — Don't use any priors. Probability parameters will simply be the observed frequencies, after relative sequence weighting. raw vibes only.
  • --plaplace — Use a Laplace +1 prior in place of the default mixture Dirichlet prior.

🌸 OPTIONS CONTROLLING E-VALUE CALIBRATION 🌸

Estimating the location parameters for the expected score distributions for MSV filter scores, Viterbi filter scores, and Forward scores requires three short random sequence simulations. yeah. three lil simulations. 😌

  • --EmL <n> — Sets the sequence length in simulation that estimates the location parameter mu for MSV filter E-values. Default is 200.
  • --EmN <n> — Sets the number of sequences in simulation that estimates the location parameter mu for MSV filter E-values. Default is 200.
  • --EvL <n> — Sets the sequence length in simulation that estimates the location parameter mu for Viterbi filter E-values. Default is 200.
  • --EvN <n> — Sets the number of sequences in simulation that estimates the location parameter mu for Viterbi filter E-values. Default is 200.
  • --EfL <n> — Sets the sequence length in simulation that estimates the location parameter tau for Forward E-values. Default is 100.
  • --EfN <n> — Sets the number of sequences in simulation that estimates the location parameter tau for Forward E-values. Default is 200.
  • --Eft <x> — Sets the tail mass fraction to fit in the simulation that estimates the location parameter tau for Forward evalues. Default is 0.04.

✨ OTHER OPTIONS ✨

okay bestie, this is the grab-bag section and it's lowkey iconic 💅 these are the little knobs that don't fit anywhere else but absolutely matter.

  • --nonull2 — turns OFF the null2 score corrections, which are the ones that fix up biased composition. so yeah, no more null2 for you 🚫
  • -Z <x> — tells the program "hey, the total number of targets in my searches is <x>" and it uses that number for the per-sequence E-value calculations, instead of just counting how many targets it actually saw. it's giving "trust me bro, this is the real total" 💁
  • --domZ <x> — same idea as -Z but for the per-domain conditional E-value calculations: you assert that the total number of targets in your searches is <x>, rather than using the number of targets that actually passed the reporting thresholds. so not "how many survived" but "how many there were" 🌸
  • --seed <n> — seeds the random number generator with <n>, an integer >= 0. here's the tea:
    • if <n> is >0, any stochastic simulations are reproducible — same command, same results, fr fr ✨
    • if <n> is 0, the RNG gets seeded arbitrarily and stochastic simulations will vary from run to run of the same command 😭
    • default seed is 42 (of course it is, iconic number, no cap) 🎀
  • --qformat <s> — asserts that the input query seqfile is in format <s>, which bypasses format autodetection. common choices for <s> include: fasta, embl, genbank. alignment formats also work; common choices include: stockholm, a2m, afa, psiblast, clustal, phylip. and here's the important part bestie: jackhmmer always uses a single sequence query to start its search, so when the input seqfile is an alignment, jackhmmer reads it one unaligned query sequence at a time, not as an alignment 👀. for more info, and for codes for some less common formats, see the main documentation. also the string <s> is case-insensitive (so fasta or FASTA both work) 💖
  • --tformat <s> — asserts that the input target sequence seqdb is in format <s>. same accepted choices as --qformat above, go look 👆
  • --cpu <n> — sets the number of parallel worker threads to <n>. on multicore machines the default is 2. you can also control this number by setting an environment variable, HMMER_NCPU. and heads up: there's also a master thread, so the actual number of threads that HMMER spawns is <n>+1 🧵
this option is NOT available if HMMER was compiled with POSIX threads support turned off. rip 💀
  • --stall — for debugging the MPI master/worker version: it pauses after start so the developer can attach debuggers to the running master and worker(s) processes. send a SIGCONT signal to release the pause. (under gdb: (gdb) signal SIGCONT) — only available if optional MPI support was enabled at compile-time 🛠️
  • --mpi — runs under MPI control with master/worker parallelization (using mpirun, for example, or equivalent). only available if optional MPI support was enabled at compile-time 🏃♀️

👀 SEE ALSO 👀

  • see hmmer(1) for a master man page with a list of all the individual man pages for programs in the HMMER package 📚
  • for complete documentation, see the user guide that came with your HMMER distribution (Userguide.pdf); or see the HMMER web page (http://hmmer.org/) 🌐

copyright (C) 2023 Howard Hughes Medical Institute. freely distributed under the BSD open source license 🫶

for additional information on copyright and licensing, see the file called COPYRIGHT in your HMMER source distribution, or see the HMMER web page (http://hmmer.org/).

🌟 AUTHOR 🌟

http://eddylab.org — that's the whole author line bestie, it's giving mysterious and I respect it 💜

and that's the whole page, we ate and left no crumbs 👑 this was HMMER 3.4, dated Aug 2023, and it's the jackhmmer(1) man page. go run something slay-worthy with it, ily, bye 🦋✨

📦 where this page comes from