使用python提取蛋白质序列的相似位点

这是本人的第一篇博客，以后还会更新一些类似学习笔记的东西

话不多少，直接上代码

#!/usr/bin/python
import numpy as np
from collections import Counter

with open('A1.fasta','r') as A1:             #打开蛋白质序列
	with open('senA1.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('.')     #用‘.’代替蛋白质名字信息

A1.close()
A1sen.close()
A1sen= open('senA1.fasta')            

ntxt=A1sen.read()                 #全部读到一个字符串里
txt=ntxt[:0]+ntxt[1:]             #去掉开头的一个点
slist=[0 for col in range(72)]     #创建含72个string的列表

for i in range (72):                #按点‘.’分割序列

	slist[i]= txt.split('.')[i]
comlist=str()                    #创建一个而用来存放矩阵竖列的变量
lsen=str()                          #创建一个用来存放特征序列的变量


for j in range(len(slist[0])):    
	for i in range(72):
		comlist=comlist+slist[i][j]

	if ( Counter(comlist).most_common(1)[0][0]!='-'):            #找出列字符串里出现频率最高的氨基酸
		if(Counter(comlist).most_common(1)[0][1]>36):            #大于一半的
			lsen=lsen+Counter(comlist).most_common(1)[0][0]      #一个一个存放起来
	comlist='-'                                                        #重置compare list

lsen=lsen[:100]                                          #取前100个
print len(lsen)                                           
print lsen
with open('A1feature.fasta','w')as feature:
	feature.write(lsen)                          #创建个fasta文件存起来

这里A1是一类蛋白质的序列，是我之前用muscle跑过的，所以长短都是一样的，是790个，有很多的gap（-），不知道muscle的同学可以上网搜一下，是一个多序列比对软件。

代码的功能是提取了这一类蛋白质前100个同一位置的最相似氨基酸序列，当然也可以不设置这100个限制

关于代码我详细的写了注释，欢迎各位研究生物信息学的大神一起讨论~~~