在基因家族分析中，hmmsearch用來在很多候選序列中尋找具有某種基因家族結(jié)構域的蛋白。在尋找時，需要提供基因家族對應的隱馬模型。前段時間給大家介紹了pfam數(shù)據(jù)庫更新后hmm文件不能使用的問題，有粉絲給出了如下解決方案，下載的文件是壓縮格式，將文件添加.gz后綴后，解壓即可。hmmsearch之前介紹過，今天主要給大家介紹輸出結(jié)果。
part1. 參數(shù)說明
# hmmsearch :: search profile(s) against a sequence database
# HMMER 3.0 (March 2010); http://hmmer.org/
# Copyright (C) 2010 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query HMM file: query.hmm
# target sequence database: target.fasta
# sequence reporting threshold: E-value <= 1e-05
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
這里面主要是說明用了哪個模型、輸入序列以及設置的參數(shù)閾值，這里只設置了evalue一個參數(shù)，E-value <= 1e-05。
part2. 篩選結(jié)果
Query: target.test.aln [M=583]
Scores for complete sequences (score includes all domains):
--- full sequence --- --- best 1 domain --- -#dom-
E-value score bias E-value score bias exp N Sequence Description
------- ------ ----- ------- ------ ----- ---- -- -------- -----------
6.5e-165 551.5 0.1 8.6e-165 551.1 0.0 1.0 1 test.query01G2293 gene=test.query01G2293
1.8e-162 543.5 0.0 2.3e-162 543.1 0.0 1.1 1 test.query05G2589 gene=test.query05G2589
6.8e-162 541.6 0.0 9.1e-162 541.2 0.0 1.1 1 test.query03G0799 gene=test.query03G0799
4.1e-159 532.4 0.0 5.4e-159 532.0 0.0 1.0 1 test.query07G1076 gene=test.query07G1076
8.6e-158 528.0 0.0 1.1e-157 527.7 0.0 1.1 1 test.query07G1966 gene=test.query07G1966
使用hmmsearch篩選的時候目的很明確，就是確定哪些序列中包含目標結(jié)構域，這個統(tǒng)計表很直觀。統(tǒng)計結(jié)果中會統(tǒng)計全部序列和結(jié)構域比對的結(jié)果，每個都有對應的E-value、score和bias三個統(tǒng)計結(jié)果，最后是擁有該結(jié)構域的候選序列。
part3. 比對結(jié)果
這部分首先是輸出每條序列中篩選獲得的所有結(jié)構域比對統(tǒng)計結(jié)果，之后給出每個結(jié)構域與query序列詳細的序列比對結(jié)果。
- 統(tǒng)計結(jié)果
>> test.query01G2293 gene=test.query01G2293
# score bias c-Evalue i-Evalue hmmfrom hmm to alifrom ali to envfrom env to acc
--- ------ ----- --------- --------- ------- ------- ------- ------- ------- ------- ----
1 ! 10.7 0.0 0.00082 0.4 212 263 .. 9 61 .. 2 75 .. 0.88
2 ! 41.6 0.1 3.5e-13 1.7e-10 310 473 .. 97 261 .. 90 266 .. 0.84
3 ! 21.4 0.0 4.6e-07 0.00022 507 582 .. 286 361 .. 282 364 .. 0.90
上面的結(jié)果顯示test.query01G2293序列中存在3個結(jié)構域。
第一列的！表示該結(jié)構域全部或者部分存在于比對全序列及結(jié)構中，如果不存在，表示為？。
對于c-Evalue和i-Evalue，c表示conditional，計算時假定全序列都是同源序列進行比對，i表示independent E-value，無前面的假設。part2比對結(jié)果中的E-value為i-Evalue。
后面幾列是比對區(qū)間，hmm為模型文件中的比對區(qū)間，ali為輸入序列中的比對區(qū)間，env為輸入序列中結(jié)構域的區(qū)間。
acc為準確度
- 比對結(jié)果
之前沒有注意過這個結(jié)果，發(fā)現(xiàn)這里面的東西還挺多，詳細說一下。
Alignments for each domain:
== domain 1 score: 10.7 bits; conditional E-value: 0.00082
target.test.aln 212 lrelLvecaeavsednlelaqellarlselsSpegn.pseRlaayfaeaLear 263
+lL++ca+a++ ++l+ a ll + l+ + + Rl++yfa+aL r
test.query01G2293 9 ALRLLLSCAKAIEDGDLKSADALLHIILVLADERPYlYESRLVKYFADALVRR 61
5689***********************9999988887889**********987 PP
== domain 2 score: 41.6 bits; conditional E-value: 3.5e-13
target.test.aln 310 ilealenedrvHiiDfdisqGlQwpsLlqeLasrsnkspslriTgigepesgvrseeglretgdrLaqfaeelnvpfe..f 388
i +al ++ r+H iDf i + + s l++L + s+ + +r+ i +p + + e ++ + L++ a +lnv++e
test.query01G2293 97 IDDALMGNRRLHLIDFSIPYENFEGSVLRTLPTFSGDPLPVRVSYILPPFLK-KYVESFSQM-EFLTKDAMNLNVKLEaeL 175
558999999**************************************99998.445555544.679999999998886115 PP
target.test.aln 389 navvs.kledirlesLkv...kegEtvaVnlvfqlhklldesVtesardelLrlikslnPkvvvlaeqeantndasFltrf 465
++v l +++ +L++ +e+E v+V f+l kll ++ +a+++ L ++k++nP +v++ + n+ + Flt +
test.query01G2293 176 KVVYAnSLAEVDEYKLDFkrrREDEMVVVYYKFKLDKLLTDG---KAMERELVRLKEINPTIVIMLDFYSNHTHSNFLTCL 253
56666677778888887633369*****************96...5555556679************************** PP
target.test.aln 466 iealeyYs 473
++ +yYs
test.query01G2293 254 EHSFQYYS 261
*******9 PP
==后面的內(nèi)容為每個結(jié)構域的統(tǒng)計結(jié)果，下方是序列的比對結(jié)果。
第一行為數(shù)據(jù)庫中的序列，.表示在輸入序列中存在插入，大寫的氨基酸表示該位置高度保守(emission probability比較高，隱馬模型中術語，具體原理及算法有興趣的可以查閱相關資料)。第三行為輸入的待測序列，-表示相對于數(shù)據(jù)庫序列，輸入序列中存在缺失。
最后一行為后驗概率值，區(qū)間從0-9，0表示0-5%，1表示5-15%，以此類推，*表示95-100%。
兩條序列中間為一致性序列(consensus sequence)，字母的大小寫與數(shù)據(jù)庫中序列保持一致，+表示保守型置換(conservative substitution)。這里并不是說一對相同的氨基酸只要出現(xiàn)就被定義為保守置換(示例如下)。
17 sllsssrk.........................llsdalaaqakeksl.....................vqhP......lr 45
+ +s++ + +++d+la+qa+eksl v+ P +
79 DTFSPTDDtdvsdtvlkyisqvlleeemeekpcMFHDSLALQAAEKSLyevlgesypprdqapvcvdpsVESPdncsfgTS 159
556666666666666666666666677778888********************8888885555544444444476666633 PP
第二位和倒數(shù)第二位都是l-->T，第二個為保守性置換。這是因為在計算時，只有該位點emission probability大于背景頻率才會定義成保守性置換。
part4. summary結(jié)果
Internal pipeline statistics summary:
-------------------------------------
Query model(s): 1 (583 nodes)
Target sequences: 40960 (16088871 residues)
Passed MSV filter: 2227 (0.0543701); expected 819.2 (0.02)
Passed bias filter: 1064 (0.0259766); expected 819.2 (0.02)
Passed Vit filter: 179 (0.00437012); expected 41.0 (0.001)
Passed Fwd filter: 90 (0.00219727); expected 0.4 (1e-05)
Initial search space (Z): 40960 [actual number of targets]
Domain search space (domZ): 84 [number of targets reported over threshold]
# CPU time: 2.95u 0.08s 00:00:03.03 Elapsed: 00:00:00.51
# Mc/sec: 18355.80
對于整個篩選過程，先后經(jīng)過MSV、bias、Vit以及Fwd過濾，具體原理及算法可以查看官方文檔[1]。summary中Domain search space (domZ)為最終篩選后獲得的序列，有84條。
