hmmsearch寻找相似序列

2021-07-23 本文已影响0人生信师姐

HMMER是基于隐马尔可夫模型，用于生物序列分析工作的一个非常强大的软件包，它的一般用途是识别同源蛋白或核苷酸序列和进行序列比对。与BLAST、FASTA等序列比对和数据库搜索工具相比，HMMER更准确。

hmmsearch is used to search one or more profiles against a sequence database. For each profile in hmmfile, use that query profile to search the target database of sequences in seqdb, and output ranked lists of the sequences with the most significant matches to the profile.To build profiles from multiple alignments, see hmmbuild.

Either the query hmmfile or the target seqdb may be ’-’ (a dash character), in which case the query profile or target database input will be read from a stdin pipe instead of from a file. Only one input source can come through stdin, not both. An exception is that if the hmmfile contains more than one profile query, then seqdb cannot come from stdin, because we can’t rewind the streaming target database to search it with another profile.

The output format is designed to be human-readable, but is often so voluminous that reading it is impractical, and parsing it is a pain. The --tblout and --domtblout options save output in simple tabular formats that are concise and easier to parse. The -o option allows redirecting the main output, including throwing it away in /dev/null.

一、使用说明

Usage: hmmsearch [options] <hmmfile> <seqdb>
用法：hmmsearch 参数 hmm文件序列数据库

1. Options directing output 输出定向选项

-o <f> : direct output to file <f>, not stdout  
输出到文件中，而不是到标准输出

-A <f> : save multiple alignment of all hits to file <f>  
将所有命中的多序列比对输出到文件

--tblout <f> : save parseable table of per-sequence hits to file <f>  
保存每个序列的命中结果的解析表到文件中

--domtblout <f> : save parseable table of per-domain hits to file <f>  
保存每个结构域的命中结果的解析表到文件中

--pfamtblout <f> : save table of hits and domains to file, in Pfam format <f>
保存命中和结构域的表格到文件，Pfam形式

--acc : prefer accessions over names in output
在输出文件中将登录号覆盖名字

--noali : don't output alignments, so output is smaller
不要输出比对，这样输出会变得更小

--notextw : unlimit ASCII text output line width
不限制ASCII文本输出行的宽度

--textw <n> : set max width of ASCII text output lines [120] (n>=120)
设置ASCII文本输出行的最大宽度

2. Options controlling reporting thresholds 控制报告阈值的参数：

-E <x> : report sequences <= this E-value threshold in output [10.0] (x>0)
在结果中报告E值小于这个阈值的序列
-T <x> : report sequences >= this score threshold in output
在结果中报告得分大于这个阈值的序列
--domE <x> : report domains <= this E-value threshold in output [10.0] (x>0)
在结果中报告E值小于这个阈值的domain
--domT <x> : report domains >= this score cutoff in output
在结果中报告大于这个分数的domain

3. Options controlling inclusion (significance) thresholds 控制包含（显著性）阈值的参数：

--incE <x> : consider sequences <= this E-value threshold as significant
将<=此e值阈值的序列视为有意义的

--incT <x> : consider sequences >= this score threshold as significant
将得分大于这个阈值的序列设为显著的

--incdomE <x> : consider domains <= this E-value threshold as significant
将大于等于这个E值阈值的domain认为是显著的

--incdomT <x> : consider domains >= this score threshold as significant
将domain得分大于这个阈值的认为是显著的

3. Options controlling model-specific thresholding 控制模型特异性阈值的参数

--cut_ga : use profile's GA gathering cutoffs to set all thresholding
使用文件的GA gathering cutoffs来设置所有的阈值

--cut_nc : use profile's NC noise cutoffs to set all thresholding
使用文件的NC noise cutoffs来设置所有阈值

--cut_tc : use profile's TC trusted cutoffs to set all thresholding
使用文件的TC trusted cutoffs来设置所有阈值

4.Options controlling acceleration heuristics 控制加速启发式？搜索的参数

--max : Turn all heuristic filters off (less speed, more power)
关闭所有的启发式过滤器
--F1 <x> : Stage 1 (MSV) threshold: promote hits w/ P <= F1 [0.02]
--F2 <x> : Stage 2 (Vit) threshold: promote hits w/ P <= F2 [1e-3]
--F3 <x> : Stage 3 (Fwd) threshold: promote hits w/ P <= F3 [1e-5]
--nobias : turn off composition bias filter

5. Other expert options 其他的一些参数

--nonull2 : turn off biased composition score corrections
关闭偏向组成分数校正？？
-Z <x> : set # of comparisons done, for E-value calculation
--domZ <x> : set # of significant seqs, for domain E-value calculation
--seed <n> : set RNG seed to <n> (if 0: one-time arbitrary seed) [42]
--tformat <s> : assert target <seqfile> is in format <s>: no autodetection
--cpu <n> : number of parallel CPU workers to use for multithreads
用于多线程的并行CPU工作程序的数量

用法如下：

hmmsearch  \
    --domtblout Ensembl.txt   \
    PF01535.hmm \
    Zea_mays.AGPv3.24.pep.all.fa

结果如下图：

image.png

结果按照E-values值从小到大排序，形式与blast类似。

target name是每个目标序列的名称；
query name是查询序列的名称；
score是比对得分，分值越高说明越相似；
E-value目标序列的期望值（统计意义）；最重要的是E-value值，值越小，越可信，相当于一个统计量。

target length输入序列总长度；query length 输入的查询序列也就是蛋白结构域的总长度；筛选的时候不应该筛选这两个长度，没有意义；应该看后面的两个 $\color{red}{from \ to}$ 看看这两条序列从哪里到哪里比对上了；一般我们会筛选第7列的E-value值；

In protein search programs, the --domtblout option produces the domain hits table.

There is one line for each domain. There may be more than one domain per sequence.

The domain table has 22 whitespacedelimited fields followed by a free text target sequence description, as follows:

（1）target name: The name of the target sequence or profile.

（2）target accession: Accession of the target sequence or profile, or - if none is available.

（3）tlen: Length of the target sequence or profile, in residues. This (together with the query length) is useful for interpreting where the domain coordinates (in subsequent columns) lie in the sequence.

（4）query name: Name of the query sequence or profile.

（5）accession: Accession of the target sequence or profile, or - if none is available.

（6）qlen: Length of the query sequence or profile, in residues.

（7） $\color{red}{E-value}$ : E-value of the overall sequence/profile comparison (including all domains).

（8）score: Bit score of the overall sequence/profile comparison (including all domains), inclusive of a null2 bias composition correction to the score.

（9）bias: The biased composition score correction that was applied to the bit score.
（10） $\color{red}{\#}$ : This domain’s number (1..ndom).

（11） $\color{red}{of }$ :The total number of domains reported in the sequence, ndom.

（12）c-Evalue: The “conditional E-value”, a permissive measure of how reliable this particular domain may be. The conditional Evalue is calculated on a smaller search space than the independent E-value. The conditional E-value uses the number of targets that pass the reporting thresholds. The null hypothesis test posed by the conditional E-value is as follows. Suppose that we believe that there is already sufficient evidence (from other domains) to identify the set of reported sequences as homologs of our query; now, how many additional domains would we expect to find with at least this particular domain’s bit score, if the rest of those reported sequences were random nonhomologous sequence (i.e. outside the other domain(s) that were sufficient to identified them as homologs in the first place)?

（13）i-Evalue: The “independent E-value”, the E-value that the sequence/profile comparison would have received if this were the only domain envelope found in it, excluding any others. This is a stringent measure of how reliable this particular domain may be. The independent E-value uses the total number of targets in the target database.

（14）score: The bit score for this domain.

（15）bias: The biased composition (null2) score correction that was applied to the domain bit score.

（16）from (hmm coord): The start of the MEA alignment of this domain with respect to the profile, numbered 1..N for a profile of N consensus positions.

（17）to (hmm coord): The end of the MEA alignment of this domain with respect to the profile, numbered 1..N for a profile of N consensus positions.

（18） $\color{red}{from\ (ali\ coord)}$ : The start of the MEA alignment of this domain with respect to the sequence, numbered 1..L for a sequence of L residues.

（19） $\color{red}{to\ (ali\ coord)}$ : The end of the MEA alignment of this domain with respect to the sequence, numbered 1..L for a sequence of L residues.

（20）from (env coord): The start of the domain envelope on the sequence, numbered 1..L for a sequence of L residues. The envelope defines a subsequence for which their is substantial probability mass supporting a homologous domain, whether or not a single discrete alignment can be identified. The envelope may extend beyond the endpoints of the MEA alignment, and in fact often does, for weakly scoring domains.

（21）to (env coord): The end of the domain envelope on the sequence, numbered 1..L for a sequence of L residues.

（22）acc: The mean posterior probability of aligned residues in the MEA alignment; a measure of how reliable the overall alignment is (from 0 to 1, with 1.00 indicating a completely reliable alignment according to the model).

（23）description of target: The remainder of the line is the target’s description line, as free text. As with the target hits table (above), this table is columnated neatly for human readability, but you should not write parsers that rely on this columnation; parse based on space-delimited fields instead.