最近在探索全基因組基因家族的分析方法，而找到某一家族需要通過(guò)保守結(jié)構(gòu)域來(lái)預(yù)測(cè)，從而找到物種的某一基因家族，從而進(jìn)行之后的分析。這里就需要用到HMMER3.1軟件，鑒定物種某一基因家族。下面開始HMMER3.1基本使用教程。
1.軟件的下載與安裝
系統(tǒng)環(huán)境：Ubuntu16.04 LTS
HMMER3.1官方下載地址：http://hmmer.org/download.html。
HMMER3.1使用手冊(cè)：http://eddylab.org/software/hmmer3/3.1b2/Userguide.pdf。
# 直接下載匹配OS的二進(jìn)制包，根本就不需要進(jìn)行安裝，只要稍微設(shè)置一下PATH變量就可以使用了，非常方便。
# 在home目錄下創(chuàng)建biosoft目錄，一般的生物分析軟件都在這里
mkdir biosoft
cd biosoft
# 下載二進(jìn)制包解壓，速度好像有些慢，只能無(wú)奈等待
wget http://eddylab.org/software/hmmer3/3.1b2/hmmer-3.1b2-linux-intel-x86_64.tar.gz
tar -zxvf hmmer-3.1b2-linux-intel-x86_64.tar.gz
# 將binaries目錄添加到環(huán)境變量中，能夠全局環(huán)境使用
cd ~
# 我這里使用的shell是zsh，有些不太一樣，一般需要修改.bashrc文件，我這里是.zshrc，打開文件在最后面添加一句
vim ./zshrc
##(添加的語(yǔ)句) PATH=$PATH:~/biosoft/hmmer-3.1b2-linux-intel-x86_64/binaries
source .zshrc
hmmbuild/hmmsearch/hmmscan/hmmalign are the core functionalities for protein domain analysis and annotation pipelines.
主要的蛋白質(zhì)domain分析和注釋管道流就是這四個(gè)二進(jìn)制文件——hmmbuild/hmmsearch/hmmscan/hmmalign
2.軟件使用
本文使用的數(shù)據(jù)例子都是安裝包中提供的，在
~/biosoft/hmmer-3.1b2-linux-intel-x86_64/tutorial目錄中。
2.1 hmmbuild構(gòu)建HMM文件
該程序需要提供sto格式的序列比對(duì)文件，格式如圖所示，它的主要區(qū)別就是# STOCKHOLM 1.0開頭和//結(jié)尾。
一般我們要去獲得已知基因組的某個(gè)家族成員，常用的方法就是從已知的某個(gè)物種的基因家族成員的保守結(jié)構(gòu)序列去比對(duì)。而HMMER3.1需要sto格式的序列比對(duì)結(jié)果，因此我們要將比對(duì)之后的序列轉(zhuǎn)換成sto格式。在線工具:sequence conversion可以將很多格式的比對(duì)結(jié)果轉(zhuǎn)換成sto格式。需要注意的是：比對(duì)后的序列長(zhǎng)度一定要一致，否則沒(méi)有辦法進(jìn)行轉(zhuǎn)換。
有了比對(duì)文件之后，可以開始構(gòu)建HMM文件，hmmbuild +要輸出的文件名（.hmm）+sto文件（相對(duì)路徑），用法如下：
hmmbuild globins4.hmm tutorial/globins4.sto，然后會(huì)得到如下信息：
HMM文件已經(jīng)生成了，我們來(lái)看一下，其中的內(nèi)容，但是好像看不懂，使用手冊(cè)里說(shuō)這是ASCII碼文件，之后還會(huì)提到。解釋一下其中的idx這一行，nseq表示一共四條序列，alen表示比對(duì)氨基酸一共171個(gè)，mlen表示最大比對(duì)上149個(gè)氨基酸，eff_nseq表示比對(duì)效率0.96存在22個(gè)gap，re/pos表示每個(gè)位置的相對(duì)熵0.589（這個(gè)就表示看不太懂）。
2.2 hmmsearch搜索序列數(shù)據(jù)庫(kù)
首先，得有一個(gè)本地的數(shù)據(jù)庫(kù)讓你來(lái)搜索，這個(gè)蛋白質(zhì)數(shù)據(jù)文件比較大，是fasta格式，uniprot sprot.fasta，我已經(jīng)是提前下載好了,解壓之后大概300M左右的大小，但是你也可以不用下載，用一個(gè)比較小的現(xiàn)有子文件globin45.fa，大概有45條蛋白序列。hmmsearch可以識(shí)別的格式包括fasta，EMBL/UniProt文本格式和GENBANK格式。
用法：hmmsearch globins4.hmm uniprot sprot.fasta > globins4.out，hmm文件是上一步生成的，蛋白質(zhì)庫(kù)是下載的，重定向輸出globins4.out文件。文件內(nèi)容：
比對(duì)結(jié)果的說(shuō)明：最后兩列是對(duì)序列的說(shuō)明和描述信息，第一列的E-vaule是最終要的一個(gè)參數(shù)，越小越有可能是同源的序列，第二列的score也可以用來(lái)評(píng)估可能性，而且不依賴于比對(duì)數(shù)據(jù)庫(kù)的大小，只依賴HMM文件和目標(biāo)序列，第三列bias是score的偏差，比如score是222.7，那么原始值就是加上3.2，等于225.9，但是這個(gè)數(shù)值一般不重要，可以忽略不看。后面的3列也是一樣的參數(shù)，只是對(duì)于best 1 domain而言，前三列是對(duì)于full sequence而言。最后的#dom內(nèi)容是有多少個(gè)domain，exp是均值，可以是小數(shù)，而N是真正的整數(shù)個(gè)數(shù)。
這是最基本的使用，后續(xù)如何從比對(duì)到的序列中再進(jìn)一步進(jìn)行相關(guān)序列的挖掘還需進(jìn)一步探索。請(qǐng)各位指正，難免有錯(cuò)誤。
