nhmmscan is used to search nucleotide sequences against collections of nucleotide profiles. For each sequence in <seqfile>, use that query sequence to search the target database of profiles in <hmmfile>, and output ranked lists of the profiles with the most significant matches to the sequence.
The <seqfile> may contain more than one query sequence. It can be in FASTA format, or several other common sequence file formats (genbank, embl, and uniprot, among others), or in alignment file formats (stockholm, aligned fasta, and others). See the --qformat option for a complete list.
Table of hits
Save a simple tabular (space-delimited) file summarizing the per-target output, with one data line per homologous target model found.
Table of hits (dfam)
Save a tabular (space-delimited) file summarizing the per-hit output, similar to --tblout but more succinct.
List of per-position scores for each hit (--aliscoreout)
Save to file a list of per-position scores for each hit. This is useful, for example, in identifying regions of high score density for use in resolving overlapping hits from different models.
Reporting thresholds control which hits are reported in output files (the main output, --tblout, and --domtblout).
E-value (-E)
In the per-target output, report target profiles with an E-value of <= <x>. The default is 10.0, meaning that on average, about 10 false positives will be reported per query, so you can see the top of the noise and decide for yourself if it’s really noise.
Bit score (-T)
Instead of thresholding per-profile output on E-value, instead report target profiles with a bit score of >= <x>.
Inclusion thresholds are stricter than reporting thresholds. Inclusion thresholds control which hits are considered to be reliable enough to be included in an output alignment or a subsequent search round. In hmmscan, which does not have any alignment output (like hmmsearch or phmmer) nor any iterative search steps (like jackhmmer), inclusion thresholds have little effect. They only affect what domains get marked as significant (!) or questionable (?) in domain output.
E-value of per target inclusion threshold
Use an E-value of <= <x> as the per-target inclusion threshold. The default is 0.01, meaning that on average, about 1 false positive would be expected in every 100 searches with different query sequences.
Bit score of per target inclusion threshold
Instead of using E-values for setting the inclusion threshold, instead use a bit score of >= <x> as the per-target inclusion threshold. It would be unusual to use bit score thresholds with hmmscan, because you don’t expect a single score threshold to work for different profiles; different profiles have slightly different expected score distributions.
Curated profile databases may define specific bit score thresholds for each profile, superseding any thresholding based on statistical significance alone. To use these options, the profile must contain the appropriate (GA, TC, and/or NC) optional score threshold annotation; this is picked up by hmmbuild from Stockholm format alignment files. Each thresholding option has two scores: the per-sequence threshold <x1> and the per-domain threshold <x2> These act as if -T<x1> --incT<x1> --domT<x2> --incdomT<x2> has been applied specifically using each model’s curated thresholds.
--cut_ga
Use the GA (gathering) bit scores in the model to set per-sequence (GA1) and per-domain (GA2) reporting and inclusion thresholds. GA thresholds are generally considered to be the reliable curated thresholds defining family membership; for example, in Pfam, these thresholds define what gets included in Pfam Full alignments based on searches with Pfam Seed models.
--cut_nc
Use the NC (noise cutoff) bit score thresholds in the model to set per-sequence (NC1) and per-domain (NC2) reporting and inclusion thresholds. NC thresholds are generally considered to be the score of the highest-scoring known false positive.
--cut_tc
Use the NC (trusted cutoff) bit score thresholds in the model to set per-sequence (TC1) and per-domain (TC2) reporting and inclusion thresholds. TC thresholds are generally considered to be the score of the lowest-scoring known true positive that is above all known false positives.
MSV filter
The sequence is aligned to the profile using a specialized model that allows multiple high-scoring local ungapped segments to match. The optimal alignment score (Viterbi score) is calculated under this multi- segment model, hence the term MSV, for “multi-segment Viterbi”. This is HMMER’s main speed heuristic. The MSV score is comparable to BLAST’s sum score (optimal sum of ungapped alignment segments). Roughly speaking, MSV is comparable to skipping the heuristic word hit and hit extension steps of the BLAST acceleration algorithm.
The MSV filter is very, very fast. In addition to avoiding indel calculations in the dynamic programming table, it uses reduced precision scores scaled to 8-bit integers, enabling acceleration via 16-way parallel SIMD vector instructions.
The MSV score is a true log-odds likelihood ratio, so it obeys conjectures about the expected score distribution (Eddy, 2008) that allow immediate and accurate calculation of the statistical significance (P- value) of the MSV bit score.
By default, comparisons with a P-value of ≤ 0.02 pass this filter, meaning that about 2% of nonhomol- ogous sequences are expected to pass. You can use the --F1 option to change this threshold. For example, --F1 <0.05> would pass 5% of the comparisons, making a search more sensitive but slower. Setting the threshold to ≥ 1.0 (--F1 99 for example) assures that all comparisons will pass. Shutting off the MSV filter may be worthwhile if you want to make sure you don’t miss comparisons that have a lot of scattered insertions and deletions. Alternatively, the --max option causes the MSV filter step (and all other filter steps) to be bypassed.
The MSV bit score is calculated as a log-odds score using the null model for comparison. No correction for a biased composition or repetitive sequence is done at this stage. For comparisons involving biased sequences and/or profiles, more than 2% of comparisons will pass the MSV filter. At the end of search output, there is a line like:
Passed MSV filter: 107917 (0.020272); expected 106468.8 (0.02)
which tells you how many and what fraction of comparisons passed the MSV filter, versus how many (and what fraction) were expected.
Viterbi filter
The sequence is now aligned to the profile using a fast Viterbi algorithm for optimal gapped alignment.
This Viterbi implementation is specialized for speed. It is implemented in 8-way parallel SIMD vector instructions, using reduced precision scores that have been scaled to 16-bit integers. Only one row of the dynamic programming matrix is stored, so the routine only recovers the score, not the optimal alignment itself. The reduced representation has limited range; local alignment scores will not underflow, but high scoring comparisons can overflow and return infinity, in which case they automatically pass the filter.
The final Viterbi filter bit score is then computed using the appropriate null model log likelihood (by default the biased composition filter model score, or if the biased filter is off, just the null model score). If the P-value of this score passes the Viterbi filter threshold, the sequence passes on to the next step of the pipeline.
The --F2 <x> option controls the P-value threshold for passing the Viterbi filter score. The default is 0.001. The --max option bypasses all filters in the pipeline. At the end of a search output, you will see a line like:
Passed Vit filter: 2207 (0.00443803); expected 497.3 (0.001)
which tells you how many and what fraction of comparisons passed the Viterbi filter, versus how many were expected.
Forward filter/parser
The sequence is now aligned to the profile using the full Forward algorithm, which calculates the likelihood of the target sequence given the profile, summed over the ensemble of all possible alignments.
This is a specialized time- and memory-efficient Forward implementation called the “Forward parser”. It is implemented in 4-way parallel SIMD vector instructions, in full precision (32-bit floating point). It stores just enough information that, in combination with the results of the Backward parser (below), posterior probabilities of start and stop points of alignments (domains) can be calculated in the domain definition step (below), although the detailed alignments themselves cannot be.
The Forward filter bit score is calculated by correcting this score using the appropriate null model log likelihood (by default the biased composition filter model score, or if the biased filter is off, just the null model score). If the P-value of this bit score passes the Forward filter threshold, the sequence passes on to the next step of the pipeline.
The bias filter score has no further effect in the pipeline. It is only used in filter stages. It has no effect on final reported bit scores or P-values. Biased composition compensation for final bit scores is done by a more complex domain-specific algorithm, described below.
The --F3 <x> option controls the P-value threshold for passing the Forward filter score. The default is 1e-5. The --max option bypasses all filters in the pipeline. At the end of a search output, you will see a line like:
Passed Fwd filter: 1076 (0.00216371); expected 5.0 (1e-05)
which tells you how many and what fraction of comparisons passed the Forward filter, versus how many were expected.
Bias Filter Options
The --max option bypasses all filters in the pipeline, including the bias filter.
The --nobias option turns off (bypasses) the biased composition filter. The simple null model is used as a null hypothesis for MSV and in subsequent filter steps. The biased composition filter step compromises a small amount of sensitivity. Though it is good to have it on by default, you may want to shut it off if you know you will have no problem with biased composition hits.
Advanced Documentation
A more detailed look at the internals of the various filter pipelines was posted on the developer's blog. The information posted there may be useful to those who are struggling with poor-scoring sequences.
The next number, the bias, is a correction term for biased sequence composition that has been applied to the sequence bit score.1 For instance, for the top hit MYG PHYCA that scored 222.7 bits, the bias of 3.2 bits means that this sequence originally scored 225.9 bits, which was adjusted by the slight 3.2 bit biased- composition correction. The only time you really need to pay attention to the bias value is when it’s large, on the same order of magnitude as the sequence bit score. Sometimes (rarely) the bias correction isn’t aggressive enough, and allows a non-homolog to retain too much score. Conversely, the bias correction can be too aggressive sometimes, causing you to miss homologs. You can turn the biased-composition score correction off with the --nonull2 option (and if you’re doing that, you may also want to set --nobias, to turn off another biased composition step called the bias filter, which affects which sequences get scored at all).
nonull2
can be too aggressive sometimes, causing you to miss homologs. You can turn the biased-composition score correction off with the --nonull2 option (and if you’re doing that, you may also want to set --nobias, to turn off another biased composition step called the bias filter, which affects which sequences get scored at all).
domZ
Assert that the total number of targets in your searches is <x>, for the purposes of per-domain conditional E-value calculations, rather than the number of targets that passed the reporting thresholds.
Z
Assert that the total number of targets in your searches is <x>, for the purposes of per-sequence E-value calculations, rather than the actual number of targets seen.
Window length tail mass (--w_beta)
The upper bound, W, on the length at which nhmmer expects to find an instance of the model is set such that the fraction of all sequences generated by the model with length >= W is less than <x>. The default is 1e-7.
Model instance length upper bound (--w length)
Override the model instance length upper bound, W, which is otherwise controlled by --w beta. It should be larger than the model length. The value of W is used deep in the acceleration pipeline, and modest changes are not expected to impact results (though larger values of W do lead to longer run time).
Seed the random number generator with <n>, an integer >= 0. If <n> is nonzero, any stochastic simulations will be reproducible; the same command will give the same results. If <n> is 0, the random number generator is seeded arbitrarily, and stochastic simulations will vary from run to run of the same command.
This Galaxy tool relies on HMMER3 from http://hmmer.janelia.org/ Internally the software is cited as:
# hmmscan :: search sequence(s) against a profile database # HMMER 3.1 (February 2013); http://hmmer.org/ # Copyright (C) 2011 Howard Hughes Medical Institute. # Freely distributed under the GNU General Public License (GPLv3). # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
The wrappers were written by Eric Rasche and is licensed under Apache2. The documentation is copied from the HMMER3 documentation.