【HMMSCAN】使用pfam数据库对多序列文件进行结构域注释

写在前面

做基因功能的人都会特别注意基因上有什么功能结构域，通常我们认为，结构域决定了这个基因的功能。随着高通量测序技术的发展，我们完全可以通过一级序列来预测该基因的结构域，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即可。

参考链接：

1.http://www.chenlianfu.com/?p=2274

2.对输出结果的解读

https://www.omicsclass.com/article/499