生信初级python 序列长度，domain数量

使用hmmer计算出domain之后，进行处理

任务

1. 计算出蛋白质氨基酸序列中含有同一domain的数量并将大于10个的domain筛选出来

2. 计算出氨基酸数量

#任务1

hmmerscan 跑出的结果是这样的以.dom结尾的文件，我们要提取jgi|Pecora1|*的蛋白质号

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) #重新写入清理完的文件

清理完之后可以得到

干净的名称，每一个重复都是一个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 计算序列长度

这样的序列我们要计算他们的长度 可以用 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)

################

结果一览

长度

domain数量