这段代码的意思是,先找出每一竖列出现次数最多的碱基,再根据出现的频率进行排序,再用特征序列的碱基位置找到原序列真实的碱基,形成一个训练集,为之后的神经网络训练做准备。

我这里使用的蛋白质序列都是事先用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()




本人是正在做毕设的低水平菜鸟。

欢迎生物信息学的各位大佬评论和指教~

Logo

CSDN联合极客时间,共同打造面向开发者的精品内容学习社区,助力成长!

更多推荐