全基因组基因家族的分析系列之HMMER3.1使用
最近在探索全基因组基因家族的分析方法，而找到某一家族需要通过保守结构域来预测。
从而找到物种的某一基因家族，从而进行之后的分析。
这里就需要用到HMMER3.1软件，鉴定物种某一基因家族。下面开始HMMER3.1基本使用教程。
1. 软件的下载与安装
系统环境：Ubuntu16.04 LTS
HMMER3.1官方下载地址：http://hmmer.org/download.html。
HMMER3.1使用手册：http://eddylab.org/software/hmmer3/3.1b2/Userguide.pdf。
# 直接下载匹配OS的二进制包，根本就不需要进行安装，只要稍微设置一下PATH变量就可以使用了，非常方便。
# 在home目录下创建biosoft目录，一般的生物分析软件都在这里
mkdir biosoft
cd biosoft
# 下载二进制包解压，速度好像有些慢，只能无奈等待
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目录添加到环境变量中，能够全局环境使用
cd ~
# 我这里使用的shell是zsh，有些不太一样，一般需要修改.bashrc文件，我这里是.zshrc，打开文件在最后面添加一句
vim ./zshrc
##(添加的语句) 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.
主要的蛋白质domain分析和注释管道流就是这四个二进制文件——hmmbuild/hmmsearch/hmmscan/hmmalign
2.软件使用
本文使用的数据例子都是安装包中提供的，在
~/biosoft/hmmer-3.1b2-linux-intel-x86_64/tutorial目录中。
2.1 hmmbuild构建HMM文件
该程序需要提供sto格式的序列比对文件，格式如图所示，它的主要区别就是# STOCKHOLM 1.0开头和//结尾。
一般我们要去获得已知基因组的某个家族成员，常用的方法就是从已知的某个物种的基因家族成员的保守结构序列去比对。而HMMER3.1需要sto格式的序列比对结果，因此我们要将比对之后的序列转换成sto格式。在线工具:sequence conversion可以将很多格式的比对结果转换成sto格式。需要注意的是：比对后的序列长度一定要一致，否则没有办法进行转换。
有了比对文件之后，可以开始构建HMM文件，hmmbuild +要输出的文件名（.hmm）+sto文件（相对路径），用法如下：
hmmbuild globins4.hmm tutorial/globins4.sto
然后会得到如下信息：
解释一下其中的idx这一行，nseq表示一共四条序列，alen表示比对氨基酸一共171个，mlen表示最大比对上149个氨基酸，eff_nseq表示比对效率0.96存在22个gap，re/pos表示每个位置的相对熵0.589（这个就表示看不太懂）。
HMM文件已经生成了，我们来看一下，其中的内容，但是好像看不懂，使用手册里说这是ASCII码文件，之后还会提到。
2.2 hmmsearch搜索序列数据库
首先，得有一个本地的数据库让你来搜索，这个蛋白质数据文件比较大，是fasta格式，uniprot sprot.fasta，我已经是提前下载好了,解压之后大概300M左右的大小，但是你也可以不用下载，用一个比较小的现有子文件globin45.fa，大概有45条蛋白序列。hmmsearch可以识别的格式包括fasta，EMBL/UniProt文本格式和GENBANK格式。
用法：hmmsearch globins4.hmm uniprot sprot.fasta > globins4.out，hmm文件是上一步生成的，蛋白质库是下载的，重定向输出globins4.out文件。文件内容：
比对结果的说明：最后两列是对序列的说明和描述信息：
第一列的E-vaule是最终要的一个参数，越小越有可能是同源的序列
第二列的score也可以用来评估可能性，而且不依赖于比对数据库的大小，只依赖HMM文件和目标序列
第三列bias是score的偏差，比如score是222.7，那么原始值就是加上3.2，等于225.9，但是这个数值一般不重要，可以忽略不看。
后面的3列也是一样的参数，只是对于best 1 domain而言，前三列是对于full sequence而言。最后的#dom内容是有多少个domain，exp是均值，可以是小数，而N是真正的整数个数。
这是最基本的使用，后续如何从比对到的序列中再进一步进行相关序列的挖掘还需进一步探索。请各位指正，难免有错误。