fas文件序列提取:
# fas文件转换为字典dict();
def gene_to_dict(file):
import os
os.chdir("F:\\file")
seq_dict = {}
name = ''
seq=''
for line in open(file,"r"):
if line.startswith('>') and seq !='':
seq_dict[name]=seq
seq=''
if line.startswith('>'):
name=line[1:20]
else:
seq=seq+line.strip()
print(seq_dict)
gene_to_dict('00.fas') # file需要包含引号 gene_to_dict('00.fas')
·-----------------------
{'Zm00001d028804_T001': 'MLSGDIPPSQTIYLNNLNEPMAISYAKK', }
# 改进版本——提取基因蛋白序列(以玉米V4版本基因为例;[:15] 提取第一个转录本;[:20]提取全部转录本])
import os
os.chdir("F:\\file")
seq_dict = {}
name = ''
seq=''
file2 = open(r"name.fas","r") # 输入基因name
file3 = open(r"out.fas","w") # 输出文件
for line in open(r'Zea_mays.AGPv4.pep.all.fas',"r"): #玉米V4蛋白数据库文件;
if line.startswith('>') and seq !='':
seq_dict[name]=seq
seq=''
if line.startswith('>'):
name=line[0:20] ## 提取玉米V4版本基因 不同转录本[0:20]; 提取第一个基因 [0:15]
else:
seq=seq+line.strip()
#print(seq_dict,type(seq_dict))
# seq_dict为dict;
n=0
for line in file2:
line2=line.strip()
n+=1
for key,value in seq_dict.items():
if line2 in key and line2 !='':
print(key,value)
file3.write(key+'\n')
file3.write(value+'\n')
print(n)
print('序列提取完毕!')
file2.close()
file3.close()
-------------------------------
>Zm00001d033955_T001
MAAASSAAPARHSCAKLSVPVEDPKAVTAGGGTVFVKATWLPSRFSLAVTDGAGAWVADASDHEVRLRAEQWDQPVADYIALAERYLAFQQPGSTYSFHDAGKGQRRLAWTFERQGTKLEWRWKLQPSPNTQQTISEILDFLMDANIRLSEEVVRKTQSFDKLKQEAEKCLQQSERFNNEKAEFEQAAFSKFVAVLNSKKAKLRQLRDKVVELEPAVKPPKEEAGQEQEEEEENSTDRTELFEGENDKEASAKDEHSGETGSGNVPSSPGESAATSRGRGRGRGRKKARR
总结:
1. if ... in ...: 结构可以用来提取基因不同的转录本(全部还是取还是取第一个;)
2. 缺点:最后一个基因无法转换为dict;(在库文件最后一行加上一个 > ;完美解决转换问题)
网友评论