这段代码的意思是,先找出每一竖列出现次数最多的碱基,再根据出现的频率进行排序,再用特征序列的碱基位置找到原序列真实的碱基,形成一个训练集,为之后的神经网络训练做准备。
我这里使用的蛋白质序列都是事先用muscle跑过的(不知道muscle的同学可以搜一下‘多序列比对软件muscle ’)
#!/usr/bin/python
#coding=utf-8
import string
import numpy as np
from collections import Counter
nb_seq=0
with open('A6.fasta','r') as A1:
with open('sequence.fasta','w') as A1sen:
for line in A1.readlines():
line=line.strip('\n')
if'>' not in line:
A1sen.write(line)
if'>'in line :
A1sen.write('.')
nb_seq+=1
A1.close()
A1sen.close()
A1sen= open('sequence.fasta')
ntxt=A1sen.read()
txt=ntxt[:0]+ntxt[1:]
slist=[0 for col in range(nb_seq)]
for i in range (nb_seq):
slist[i]= txt.split('.')[i] #原序列
comlist=str()
lsen=str()
#print slist
wei_list=[[0 for col in range(3)] for row in range(1000)]#一个二维列表,存放碱基,碱基出现的频率以及位置
k=0
for j in range(len(slist[0])): #这一个大循环给二维数组赋了值
for i in range(nb_seq):
comlist=comlist+slist[i][j]
if ( Counter(comlist).most_common(1)[0][0]!='-'):
if(Counter(comlist).most_common(1)[0][1]>(nb_seq/2)):
lsen=lsen+Counter(comlist).most_common(1)[0][0]
wei_list[k][0]=Counter(comlist).most_common(1)[0][0]
wei_list[k][1]=Counter(comlist).most_common(1)[0][1]
wei_list[k][2]=j
k+=1
#print comlist
comlist=''#重置
wei_list.sort(reverse=True,key=lambda x:x[1])#逆序排列
seq_real=str()
with open('A6train.fasta','w')as A1train:
for i in range(nb_seq):
for j in range(100):
seq_real=seq_real+slist[i][wei_list[j][2]]
A1train.write(seq_real)#一句一句写到最终输出的文件里去
print len(seq_real)
A1train.write('\n')
seq_real=''#重置
A1train.close()
本人是正在做毕设的低水平菜鸟。
欢迎生物信息学的各位大佬评论和指教~