HMMER批量比对及结果处理

HMMER是基于隐马尔可夫模型序列比对工具，相较于blast，HMMER更为准确，尤其是在功能基因的研究方面，它可以更准确的检测到远的同源序列。
自己的分析需要对一批数据进行hmmsearch，因此写了这个python脚本，用于批量hmmer比对及后续的结果处理。

1 HMMER比对

1.1hmmbuild

使用hmmer首先自然是构建.hmm文件，如果可以从Pfam下载到待研究基因的.seed文件自然是最好的，如果不能的话，hmmer也可以接受来自blast等程序的多序列比对文件来完成构建.hmm文件。

假设我们的多序列比对文件名为nature.fas，构建.hmm文件的命令为：

hmmbuild nature.hmm nature.fas

这个nature.hmm就是hmmer基于隐马尔科夫模型构建的概要文件，我们将使用它去对数据进行下一步的搜索。

1.2 hmmsearch

hmmsearch使用的基本命令为：

hmmsearch [options] hmmfile seqdb

hmmfle就是我们上一步产生的nature.hmm文件；
seqdb是我们将要用于搜索的数据文件，最好提供全路径；
[options]是用于控制输出的选项，具体参数可以选择到hmmer官网下载使用说明查看。

通用的命令格式为：

hmmsearch -o out_put.file hmm.file seqdb

我选择设置--noali用于禁止在结果输出包含序列比对结果，并设置-E 1e-5来控制期望值，这样就可以方便后面的统计结果。

所以我们对应例子中的比对命令最后格式为：

hmmsearch -o nature.hmmout --noali -E 1e-5 nature.hmm nature

为了对多个数据文件进行批量hmmsearch，我们新建一个文件夹E:\py_work\practice\practice1\test，并将nature.hmm文件放入其中，然后在其中新建python脚本并运行。需要更改的地方有两处：一处是in_path，是我们的待比对序列数据所在目录，在windows下，需要在路径前面加“r”；还有就是将hmm文件名传入hmmfile变量中。

#执行hmmsearch的函数
#hmm文件需要放入当前文件夹 
#hmmsearch函数#############
def hmmsearch(in_path, hmmfile):
	import os
	#获取当前目录
	out_path = os.getcwd()
	#获取文件列表
	files = os.listdir(in_path)
	#遍历执行
	for file in files:
		# 获取输入输出文件的绝对路径
		in_file = os.path.join(in_path, file)
		out_file = os.path.join(out_path, file + '.hmmout')
		# hmm文件
		hmmfiles = hmmfile
		# 比对命令
		hmmsearch = "hmmsearch -o " + out_file + " --noali -E 1e-5 " + hmmfiles + " " + in_file
		os.system(hmmsearch)


##############################################
if __name__=='__main__':
	#修改位置一：数据文件夹路径
	in_path = r'E:\py_w