使用hmmer计算出domain之后,进行处理
任务
计算出蛋白质氨基酸序列中含有同一domain的数量并将大于10个的domain筛选出来
计算出氨基酸数量
#任务1
hmmerscan 跑出的结果是这样的以.dom结尾的文件,我们要提取jgi|Pecora1|*的蛋白质号
![](https://i-blog.csdnimg.cn/blog_migrate/703e5ea30c2c97f8980a188c9d117e8a.png)
os.chdir() #更改我们的工作路径到结果存放路径
dom文件预处理成只有序列名称
for file in os.listdir(os.getcwd()): #遍历工作路径的文件
if ".dom" in file: #只有dom文件可以进入循环
f = open("%s" % (file), 'r' ) #打开文件
line_ls = []
for line in f.readlines(): #遍历文件的每一列
if "#" not in line: #不管前三列
line = line.split("-")[1].strip().strip("16").strip("\t").strip(" ")
#上面处理找到可以分割最干净的”-“,再去干净一点
line_ls.append(line)
f2 = open("%s" % (file), 'w' )
for i in line_ls:
print(i,file=f2) #重新写入清理完的文件
清理完之后可以得到
![](https://i-blog.csdnimg.cn/blog_migrate/9e4d9f17c3d808f27c0a65118f3712d1.png)
干净的名称,每一个重复都是一个hmmer命中的domain
接下来我们算每个蛋白序列中命中的domain数量
#定义函数计算dom结果数量 利用字典只能存在一个key的原理清除多余的序列号 然后每遍历到一个就db[key] += 1
def count_dom(dom):
fi = open("%s" % (dom),"r")
x = str(dom)+".txt"
f4 = open("%s" % (x), "a")
db = {}
for line in fi.readlines():
line2 = line.strip("\n")
key = line2
if line2 in db:
db[key] += 1
if line2 not in db:
db[key] = 1
#print(db)
for key,value in db.items():
print(key,end="\t",file=f4)
print(value,file=f4)
#定义函数筛选dom_result的>10的序列 简单的if操作
def dom_result_filter(txt):
f = open("%s" % (txt) , "r")
filter_ls = []
for line in f.readlines():
print(line)
count = int(line.split("\t")[1])
if count >= 10:
filter_ls.append(line)
f2 = open("%s" % (txt) , "w")
for i in filter_ls:
print(i,end="",file=f2)
任务2 计算序列长度
![](https://i-blog.csdnimg.cn/blog_migrate/64c1fb1daa09736a3b2af5edef1b1574.png)
这样的序列我们要计算他们的长度 可以用 len() 计算字符串的长度
#定义函数测序列长度
def count_len(txt):
f = open("%s" % (txt), 'r')
db = {}
for line in f.readlines():
if line.startswith('>'):
keys = line.split(" ")[0].strip()
db[keys] = []
else:
db[keys].append(line.strip())
fx = str(txt)+".len.txt"
f2 = open("%s" % (fx), 'a')
for key,value in db.items():
print(key,end="\t",file=f2)
print(len(str(value)),file=f2)
################
结果一览
![](https://i-blog.csdnimg.cn/blog_migrate/68873e7eb259a1d67bf4211b6c83171a.png)
长度
![](https://i-blog.csdnimg.cn/blog_migrate/f9925e4fc01fdc400da5891c1e362982.png)
domain数量