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_work\practice\practice1\seqdata'
	#修改位置二：hmm文件名
	hmmfile = "nature.hmm"
	hmmsearch(in_path, hmmfile)

如果参数传入正确，那么执行上面的脚本后，便可以在我们当前所在的文件夹中生成一系列后缀为.hmmout的文件了。

完成后删除本文件夹中的python脚本及nature.hmm文件，方便下一步分析。

2 HMMSEARCH结果处理

在这里进行的主要是统计每个序列数据文件中有多少个序列比对上了我们的nature.hmm模型，然后计算了一下每Mb序列有多少个hits。

首先对结果文件进行观察，确定如何处理数据。通过观察可以发现.hmmout文件大体上可以分为两块：

第一部分是比对结果统计，一开始是表头部分，然后从第18行开始才是真正的数据结果；

第二部分序列的注释信息，这一部分我们这次的统计暂时用不到。

通过观察，我们所要做的统计就是将数据中结果统计的第一部分过滤出来，然后统计有多少行有效数据即可；还有就是计算一下每Mb数据中有多少个hit。分两步走即可：

第一步，过滤文件

通过观察，我们发现统计结果部分与后面的注释部分之间有两个空行，而文件一开始的表头部分有一个空行，也就是说当空行数大于2的时候，就到了注释部分，由此来作为数据过滤的依据即可。

我们新建一个文件夹E:\py_work\practice\practice1\test2，在其中新建python脚本，为in_path变量传入上一步产生的.hmmout文件所在文件夹路径，然后执行即可：

#处理hmmer比对文件结果，将结果部分提出来，注释部分
#另起一个新文件夹，在新文件夹中运行脚本，只需修改in_path即可

import os
#输入文件所在文件夹的路径
in_path = r'E:\py_work\practice\practice1\test'
#获取当前文件件路径
out_path = os.getcwd()
#获取待处理文件全路径
files = os.listdir(in_path)
infiles = []
for file in files:
	out = os.path.join(in_path, file)
	infiles.append(out)

#结果处理
for file in infiles:
	#获取待待输出文件绝对路径
	file1 = file.split('\\')[-1]
	outfile1 = os.path.join(out_path, file1 + '.chuli')
#处理hmmer结果文件，当空行数量大于2时，停止输入
	outfile = open(outfile1 ,'w')
	i = 0
	with open(file) as infile:
		for line in infile:
			if line == '\n':
				i += 1
			elif i >= 2:
				break
			else:
				outfile.write(line)
	outfile.close()

执行后，我们会在当前文件下获得一系列后缀为.hmmout.chuli的文件，这里面值保存了第一部分数据。

完成后删除本文件夹中的python脚本，方便下一步分析。

第二步，统计结果。

在这一步，我们需要获取序列文件的大小（以Mb为单位），并统计每个序列文件中的有效比对数hit，然后计算每Mb序列文件中有多少个hit。

#########################################################################################
#获取待处理文件绝对路径函数
def Get_dir(in_path):
	import os
	files = os.listdir(in_path)
	input_files = []
	for file in files:
		out = os.path.join(in_path, file)
		input_files.append(out)
	return input_files

#############定义计数函数########################
def CountNum(file):
	import pandas as pd
	# 设置参数engine='python'以正确识别正则表达式sep='\s{2,}'
	read_file = pd.read_table(file, skiprows=range(0, 16), header=None, sep='\s{2,}',
					   engine='python',
					   names=['f-E-value', 'f-score', 'f-bias', 'b-E-value', 'b-score', 'b-bias', 'exp', 'N',
							  'Sequence', 'Description'])
	conut_num = read_file.loc[:, 'f-E-value'].count()
	return conut_num

#定义名称处理函数
def ChuLiName(file):
	file = file.replace('.hmmout.chuli', '')
	return file

#计数函数
def GetCounts():
	hits = {}
	result_infiles = Get_dir(result_in_path)
	for file in result_infiles:
		# linux下要将'\\'改为'/'
		file_name1 = file.split('\\')[-1]
		file_name = ChuLiName(file_name1)
		file_count = CountNum(file)
		hits[file_name] = file_count
	return hits

######################一个获取文件大小的函数###########################
def get_FileSize(filePath):
	import os
	fsize = os.path.getsize(filePath)
	#以Mb为单位
	fsize = fsize / float(1024 * 1024)
	#保留小数点后两位
	return round(fsize, 2)

#获取文件大小并存进字典
def GetSize():
	sizes = {}
	data_infiles = Get_dir(data_in_path)
	for file in data_infiles:
		# linux下要将'\\'改为'/'
		file_name = file.split('\\')[-1]
		file_size = get_FileSize(file)
		sizes[file_name] = file_size
	return  sizes

##################结果输出，计算比例#########################
def JieGuo(size, hit):
	import pandas as pd
	size = pd.Series(size)
	hit = pd.Series(hit)
	result = pd.DataFrame([size, hit], index=['size', 'hit']).T
	result['per'] = result['hit'] / result['size']
	return result


#########################################################
if __name__=='__main__':
	# 输入序列数据文件所在文件夹的路径
	data_in_path = r'E:\py_work\practice\practice1\seqdata'
	# 输入过滤后文件所在文件夹路径
	result_in_path = r'E:\py_work\practice\practice1\test2'
	hits = GetCounts()
	sizes = GetSize()
	result = JieGuo(sizes, hits)
	result.to_csv("result.csv")

只最后传入序列数据文件夹路径及过滤后的文件夹路径，执行脚本后就可以在当前文件夹下获得一个result.csv文件，其内容为：

size的单位为Mb；
hit为比对到的有效序列个数，即期望值E<= 1e-5；
per为每Mb中有多少hit。