寫在前面
做基因功能的人都會特別注意基因上有什么功能結構域够挂，通常我們認為，結構域決定了這個基因的功能藕夫。隨著高通量測序技術的發(fā)展孽糖，我們完全可以通過一級序列來預測該基因的結構域枯冈，pfam和smart數(shù)據(jù)庫是一個不錯的選擇，下面我將結合我的經(jīng)驗办悟，來完成對多序列文件的結構域注釋尘奏。雖然網(wǎng)頁工具也能做這部分內容(https://www.ebi.ac.uk/Tools/hmmer/search/hmmscan)，但是如果你的蛋白質大小大于5000個氨基酸病蛉，網(wǎng)頁工具就沒法完成這個事情了炫加，畢竟太吃內存了。所以學好代碼是多么的重要铡恕，可以自己隨心所欲地搭環(huán)境琢感，跑程序。
1.工具
HMMER
2.下載
#下載安裝工具并添加到環(huán)境變量
wget http://eddylab.org/software/hmmer/hmmer.tar.gz
#下載pfam數(shù)據(jù)庫
ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam33.1/Pfam-A.hmm.gz
gzip -d Pfam-A.hmm.gz
# 得到 PFAM 數(shù)據(jù)庫的 HMM 文件探熔。 HMM 文件是文本文件驹针，需要將其變成二進制格式，以加快運算速度诀艰，同時進行壓縮柬甥，并建立成索引數(shù)據(jù)庫。
hmmpress Pfam-A.hmm
3.使用hmmscan 進行Pfam注釋
Pfam 數(shù)據(jù)庫中每個編號代表一個蛋白質家族其垄。Pfam 分 A 和 B 兩個數(shù)據(jù)庫苛蒲，其中 A 數(shù)據(jù)庫是經(jīng)過手工校正的高質量數(shù)據(jù)庫， B 數(shù)據(jù)庫雖然質量低些绿满，依然可以用來尋找蛋白質家族的保守位點臂外。Pfam 30.0以后的版本，Pfam-B就去除掉了喇颁，所以我們只關心PfamA就好了漏健。
# 因為我要查看結構域的大小，方便我去人工刪除不滿足大小的基因(如只包含結構域片段的基因)
hmmscan -o out.txt --tblout out.tbl --domtblout out.dom --noali -E 1e-5 /pfam/Pfam-A.hmm multifile.fasta
附上參數(shù)解讀
$ hmmscan [-options]
?
-h
顯示幫助信息
-o FILE
將結果輸出到指定的文件中橘霎。默認是輸出到標準輸出蔫浆。
--tblout FILE
將蛋白質家族的結果以表格形式輸出到指定的文件中。默認不輸出該文件姐叁。
--domtblout FILE
將蛋白結構域的比對結果以表格形式輸出到指定的文件中瓦盛。默認不輸出該文件。該表格中包含query序列起始結束位點與目標序列起始結束位點的匹配信息外潜。
--acc
在輸出結果中包含 PF 的編號原环，默認是蛋白質家族的名稱。
--noali
在輸出結果中不包含比對信息处窥。輸出文件的大小則會更小嘱吗。
-E FLOAT default:10.0
設定 E_value 閾值，推薦設置為 1e-5 碧库。
-T FLOAT
設定 Score 閾值柜与。
--domE FLOAT default:10.0
設定 E_value 閾值。該參數(shù)和 -E 參數(shù)類似嵌灰，不過是 domain 比對設定的值弄匕。
--cpu
多線程運行的CPU。默認應該是大于1的沽瞭，表示支持多線程運行迁匠。但其實估計一般一個hmmscan程序利用150%個CPU。并且若進行并行化調用hmmscan驹溃，當并行數(shù)高于4的時候城丧，會報錯：Fatal exception (source file esl_threads.c, line 129)。這時豌鹤，設置--cpu的值為1即可亡哄。
參考鏈接：
1.http://www.chenlianfu.com/?p=2274
2.對輸出結果的解讀
