写在前面
做基因功能的人都会特别注意基因上有什么功能结构域，通常我们认为，结构域决定了这个基因的功能。随着高通量测序技术的发展，我们完全可以通过一级序列来预测该基因的结构域，pfam和smart数据库是一个不错的选择，下面我将结合我的经验，来完成对多序列文件的结构域注释。虽然网页工具也能做这部分内容(https://www.ebi.ac.uk/Tools/hmmer/search/hmmscan)，但是如果你的蛋白质大小大于5000个氨基酸，网页工具就没法完成这个事情了，毕竟太吃内存了。所以学好代码是多么的重要，可以自己随心所欲地搭环境，跑程序。
1.工具
HMMER
2.下载
#下载安装工具并添加到环境变量
wget http://eddylab.org/software/hmmer/hmmer.tar.gz
#下载pfam数据库
ftp://ftp.ebi.ac.uk/pub/databases/Pfam/releases/Pfam33.1/Pfam-A.hmm.gz
gzip -d Pfam-A.hmm.gz
# 得到 PFAM 数据库的 HMM 文件。 HMM 文件是文本文件，需要将其变成二进制格式，以加快运算速度，同时进行压缩，并建立成索引数据库。
hmmpress Pfam-A.hmm
3.使用hmmscan 进行Pfam注释
Pfam 数据库中每个编号代表一个蛋白质家族。Pfam 分 A 和 B 两个数据库，其中 A 数据库是经过手工校正的高质量数据库， B 数据库虽然质量低些，依然可以用来寻找蛋白质家族的保守位点。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
附上参数解读
$ 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 阈值。该参数和 -E 参数类似，不过是 domain 比对设定的值。
–cpu
多线程运行的CPU。默认应该是大于1的，表示支持多线程运行。但其实估计一般一个hmmscan程序利用150%个CPU。并且若进行并行化调用hmmscan，当并行数高于4的时候，会报错：Fatal exception (source file esl_threads.c, line 129)。这时，设置–cpu的值为1即可。
参考链接：
2.对输出结果的解读