目的如题。
准备:Allwp.faa保存于D:\PyCharm\py_code\test
将以下代码写入MoveRedundancy.py中,并分段运行
创建文件夹D:\PyCharm\py_code\test\AllWP
备注:NCBI中不同的Organism可能含有相同的蛋白。Allwp.faa包含冗余蛋白序列,中蛋白名称的格式为:organism1(ProteinID1), organism2(ProteinID2)。欲对Allwp.faa去冗余,当ProteinID相同时,整合蛋白名称为格式organism1/organism2/organismX(ProteinID)。
#获取蛋白名称,并去除冗余的ProteinID,写入wpnameUni.txt中:
import os
import re
#获取唯一的蛋白ID,但是蛋白命名中organism只有一个,可以用于获取唯一蛋白序列
local = os.path.join(os.getcwd(),'AllWp.faa')#所有nr数据库下载的蛋白,有重复的蛋白
sequence = open(local,'r')#打开冗余蛋白序列
wpnameUni = open(os.path.join(os.getcwd(),'wpnameUni.txt'),'a+')#去掉重复的proteinID后的名称储存于此文件,但是organism只有一个
wpnamelist = []#创建空列表,用于储存proteinID
for linesR in sequence:
if '>' in linesR:
find1 = re.compile(r'[(](.*?)[)]', re.S)#定义查找proteinID的规则,proteinID在括号内
wpname = (re.findall(find1, linesR))#获取proteinID,wpname1为列表格式
find2 = re.compile(r'[>](.*?)[(]', re.S)#定义查找organism的规则,organism写在>与(之间
organism = (re.findall(find2, linesR))#获取organism,,organism1为列表格式
if wpname[0] not in wpnamelist:
wpnamelist.append(wpname[0])
wpnameUni.write('>' + organism[0] +'('+ wpname[0]+ ')' +'\n')
#拆分冗余序列,每个序列保存为一个文件。产生n个文件与AllWP中:
#拆分冗余序列,每个序列保存为一个文件,保存与文件夹AllWP,结果文件用于获取唯一蛋白序列
import os
import re
local = os.path.join(os.getcwd(),'AllWp.faa')#获取冗余序列路径,AllWp.faa包含所有冗余序列
resultPath =os.path.join(os.getcwd(),'AllWP')#拆分后的序列文件,储存与此路径
file = open(local,'r')#打开AllWp.faa
for line in file:#读取行
if '>' in line:#获取序列名称
find1 = re.compile(r'[>](.*?)[)]', re.S)
wpname = (re.findall(find1, line))#或许protein名称
Newfile = open(os.path.join(resultPath,wpname[0]),'a+')#在AllWP文件夹中创建以序列名称命名的文件
Newfile.write(line)#将序列名称写入新建文件
else:
Newfile.write(line)#将序列写入新建文件
#获取非冗余蛋白序列,但是每个序列只对应一个organism:
#获取唯一蛋白序列储存在WPsequenceUni.faa
import os
wpnameUni = open(os.path.join(os.getcwd(),'wpnameUni.txt'),'r')#去掉重复的proteinID后的名称储存于此文件,但是organism只有一个
WPsequenceUni = open(os.path.join(os.getcwd(),'WPsequenceUni.faa'),'a+')#去掉重复蛋白后,序列存储在此文件
local = os.path.join(os.getcwd(), 'AllWP')# 获取冗余蛋白序列储存路径
for name in wpnameUni: # 遍历去掉重复的proteinID后的名称,organism只有一个
name = name.strip() # 去掉换行符
for file in os.listdir(local):#遍历冗余蛋白序列
localfile = os.path.join(local,file)#获取冗余蛋白序列路径
s = open(localfile,'r')#打开冗余蛋白序列
line1 = s.readline()#读取第一行
line11 = line1.strip()#去掉第一行换行符
lines = s.readlines()#读取剩余行
if line11 == name:#判断名称是否相同
WPsequenceUni.write(line1)#写入
WPsequenceUni.write(''.join(lines))#''.join(lines)作用是将列表lines中元素连接成字符串
#获取非冗余蛋白序列名称,整合所有organism:
#重命名蛋白序列名称,整合多个organism,名称保存于UniWPMultiOrganism.txt
import os
import re
find1 = re.compile(r'[(](.*?)[)]', re.S)
find2 = re.compile(r'[>](.*?)[(]', re.S)
IDuni = []#唯一organism的proteinID将储存的列表
IDRed = []#冗余organism的proteinID将储存的列表
path = os.getcwd()
UniWPMultiOrganism = open(os.path.join(path,'UniWPMultiOrganism.txt'),'a+')#含有多个organism的proteinID将储存于此
with open(os.path.join(path,'wpnameUni.txt'),'r') as a:#打开唯一的蛋白ID,但是蛋白命名中organism只有一个
for line in a:
line = line.strip()
IDuni.append(line)#将唯一organism的proteinID全部存储在列表中
with open(os.path.join(path,'AllWp.faa'),'r') as b:#打开冗余基因名称
for lineRed in b:
if ">" in lineRed:
lineRed = lineRed.strip()
wpnameRed = (re.findall(find1, lineRed))#获取冗余wp
orgRed = (re.findall(find2, lineRed))#获取冗余org
for name in IDuni:#遍历非冗余基因名称
wpnameuni = (re.findall(find1, name))#获取非冗余wp
orguni = (re.findall(find2, name))#获取非冗余org
if wpnameRed[0] == wpnameuni[0]:#判断非冗余和冗余wp是否相同
i = IDuni.index(name)#获取相同出非冗余列表的索引位置
IDuni[i] = ('>' + orguni[0]+'/'+ orgRed[0] + '(' + wpnameuni[0] + ')')
for nameuniT in IDuni:
if '/' in nameuniT:#第一个organism重复,所有此处需删掉
findT = re.compile(r'[/](.*?)[)]', re.S)
T = (re.findall(findT, nameuniT))
UniWPMultiOrganism.write('>' + T[0] + ')'+ '\n')
else:
UniWPMultiOrganism.write(nameuniT + '\n')
#生成非冗余蛋白序列,蛋白名称整个所有含有此蛋白的Organism。生成UniProteinMultiOrganism.faa
import os
import re
#将名称中仅含有多个唯一organism非冗余蛋白序列文件——WPsequenceUni.faa,的蛋白名称置换成含有多个organism的——UniWPMultiOrganism.txt
#UniWPMultiOrganism.txt储存的proteinID成含有多个organism
#WPsequenceUni.faa储存非冗余蛋白序列,但是蛋白的organism只有一个
find1 = re.compile(r'[(](.*?)[)]', re.S)
path = os.getcwd()
UniProteinMultiOrganism = open(os.path.join(path,'UniProteinMultiOrganism.faa'),'a+')
with open(os.path.join(path,'WPsequenceUni.faa'),'r') as protein:
for lineP in protein:
if '>' in lineP:
lineP = lineP.strip()
proteinID = (re.findall(find1, lineP))#获取蛋白的wp
with open(os.path.join(path, 'UniWPMultiOrganism.txt'), 'r') as f:#含有多个organism的蛋白名称
for lineT in f:
lineT = lineT.strip()
togewp = (re.findall(find1, lineT))
if proteinID[0] == togewp[0]:#判断两个wp是否相等
UniProteinMultiOrganism.write(lineT + '\n')
else:
UniProteinMultiOrganism.write(lineP)