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_fileos.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 = 0with open(file) as infile:for line in infile:if line == '\n':i += 1elif i >= 2:breakelse:outfile.write(line)outfile.close()

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

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

第二步，统计结果。

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

#########################################################################################
#获取待处理文件绝对路径函数
def Get_dir(in_path):import osfiles = 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_countreturn hits######################一个获取文件大小的函数###########################
def get_FileSize(filePath):import osfsize = 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_sizereturn  sizes##################结果输出，计算比例#########################
def JieGuo(size, hit):import pandas as pdsize = pd.Series(size)hit = pd.Series(hit)result = pd.DataFrame([size, hit], index=['size', 'hit']).Tresult['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。

HMMER批量比对及结果处理相关推荐

基因家族的鉴定-基于Windows系统上的HMMER
文章首发于简书链接(https://www.jianshu.com/p/24e4ad69f3e5),发此备份基因组的序列提取,详情请看我之前的教程:https://www.jianshu.com/p ...
php批量导出pdf文件大小,php完美导出pdf，pdf合并批量导出
使用到的工具 pdftk https://www.pdflabs.com/tools/pdftk-the-pdf-toolkit/ pdf合并工具 wkhtmltopdf ...
Redis 笔记（16）— info 指令和命令行工具（查看内存、状态、客户端连接数、监控服务器、扫描大key、采样服务器、执行批量命令等）
Info 命令返回关于 Redis 服务器的各种信息和统计数值.通过给定可选的参数 section ,可以让命令只返回某一部分的信息. 1. 显示模块 server : 一般 Redis 服务器信息, ...
Redis 笔记（15）— 管道 pipeline（客户端将批量命令打包发送用来节省网络开销）
Redis 是一种基于客户端-服务端模型以及请求/响应协议的 TCP 服务.这意味着通常情况下一个请求会遵循以下步骤: 客户端向服务端发送一个查询请求,并监听 Socket 返回,通常是以阻塞模式,等 ...
Redis 笔记（03）— string类型（设置key、获取key、设置过期时间、批量设置获取key、对key进行加减、对key值进行追加、获取value子串）
字符串 string 是 Redis 最简单的数据结构.Redis 所有的数据结构都是以唯一的 key 字符串作为名称,然后通过这个唯一 key 值来获取相应的 value 数据.不同类型的数据结构的 ...
shell 批量转换文件编码
相信大家在平时的跨平台编程中碰到过文件编码问题,比如在Windows代码字符编码方式是GB2312,然而转到Linux却只支持utf-8,虽然对代码部分没啥影响,但是很多中文注释部分,却一片乱码,很让 ...
批量梯度下降(BGD)、随机梯度下降(SGD)以及小批量梯度下降(MBGD)的理解
批量梯度下降(BGD).随机梯度下降(SGD)以及小批量梯度下降(MBGD)的理解 </h1><div class="clear"></div> ...
pip 将某包指定到某目录批量安装
pip install -r requirements/requirements.txt 包批量安装 pip 将某包指定到某目录安装: pip install --target=d:\somew ...
2021年大数据HBase（十五）：HBase的Bulk Load批量加载操作
全网最详细的大数据HBase文章系列,强烈建议收藏加关注! 新文章都已经列出历史文章目录,帮助大家回顾前面的知识重点. 目录系列历史文章 HBase的Bulk Load批量加载操作一.Bulk L ...

HMMER批量比对及结果处理

HMMER批量比对及结果处理相关推荐

最新文章

热门文章