### 10. Consensus and Profile ###
按我自己的理解翻译一下:
如果只找一条链的同源,拿这条链去跟别的链一一比对,找出字符不匹配的最小的那条链就是同源链喽
但是当寻找多条链的同源的时候,这时我们需要一条所谓的“参考基准链”(就是某一位置出现概率最大的碱基组成的链,比如:给定的8条链中第一位出现概率最大的为A(出现5次),则基准链的第一位为A,第二位出现概率最大的碱基是T,则选择T,以此类推,然后将每个位置的出现概率最大的碱基连成一条链,就是参考基准链)
那个profile 是每个碱基在每条链出现的次数组成的矩阵。(可以看出A这7条链中在第一个位置出现的次数为5,第二个位置的次数为1)。应该讲明白了
我的是python3.6
[Python] 纯文本查看 复制代码 from collections import Counter
import collections
fh = open('E:/我的自学/mypyton/profile10.txt', 'wt')
rosalind=collections.OrderedDict()#这里是py3的写法,读取序列
seqLength =0
with open('E:/我的自学/mypyton/sample10.txt') as f:#提取每一个seq为line
for line in f:
line = line.rstrip()
if line.startswith('>'):
rosalindName = line
rosalind[rosalindName] = ''
continue
rosalind[rosalindName] += line
seqLength = len(rosalind[rosalindName])#每个seq的长度为8
##
##seq是每条seq的某一位置的那个碱基连成一条seq,然后计数,排序,取最大,连成一个
A,C,G,T = [],[],[],[]#初始化
consensusSeq = ''#初始化
for i in range(seqLength):
seq = ''
for k in rosalind.keys():
seq += rosalind[k]
A.append(seq.count('A'))#计数
C.append(seq.count('C'))
G.append(seq.count('G'))
T.append(seq.count('T'))
counts = Counter(seq)#为seq计数
consensusSeq += counts.most_common()[0][0]#从多到少返回,是个list啊,只返回第一个
fh.write(consensusSeq+'\n')
fh.write('\n'.join(['A:\t'+'\t'.join(map(str,A)), 'C:\t'+'\t'.join(map(str,C)),
'G:\t'+'\t'.join(map(str,G)), 'T:\t'+'\t'.join(map(str,T))]))#
##
fh.close()
consensusSeq#原来输出该是这样 |