|||
题目(点我)很好理解,难点在于python擅长行计算而不是列计算,需要一些额外的循环和存储。我的无脑解法为:
from Bio import SeqIO from collections import Counter, OrderedDict import operator def cons(file): # 用biopython里面的SeqIO模块帮助读入fasta文件并存储进一个字典,key是tag,value是序列 dic = {} i = 0 length_list = [] for record in SeqIO.parse(file,'fasta'): i += 1 # 检查是否有tag一样的序列,当然这不影响我们的计算,就是一个warning if record.id in dic: print('multiple ' + record.id + ' exist!') record.id = record.id + '_' + str(i) dic[record.id] = record.seq length_list.append(len(record.seq)) # 检查是否所有序列都一样长 length_set = set(length_list) if len(length_set) > 1: print('Sequences have different lengths!') # 在Rosalind里面肯定都是干净的数据,但在实操里只能去最短的 length = sorted(length_set)[0] # 计算每个位置的共识 con = [] dic_matrix = {'A':[], 'C':[], 'G':[], 'T':[]} for i in range(length): temp_list = [] for key in dic: temp_list.append(dic[key][i]) temp_dic = Counter(temp_list) #profile for key in dic_matrix: dic_matrix[key].append(temp_dic[key]) #consensus sorted_list = sorted(temp_dic.items(), key=operator.itemgetter(1)) con.append(sorted_list[-1][0]) # 将字典按key排序,以便输出时按照ACGT的顺序,后面强人代码可见这是没有必要的 profile_dic = OrderedDict(sorted(dic_matrix.items())) print(''.join(con)) for key in profile_dic: print(key + ': ' + ' '.join([str(a) for a in profile_dic[key]]))
好了下面是卖家秀(根据Darkstar的python2版本修改而成)
from Bio import SeqIO with open('data/cons_example.fa') as f: #读入所有序列入一个list(而不是字典) strands = [str(record.seq) for record in SeqIO.parse(f,'fasta')] #最为关键的一步,保证下面的操作都是按列而不是按行 matrix = zip(*strands) #下面这行返回一个list of dictionary, 每一个item就是每一列的profile,信息量太大,代码后详解 profile_matrix = list(map(lambda x: dict((base, x.count(base)) for base in "ACGT"), matrix)) #计算每一列的共识,max()函数还可以定制key呢 consensus = [max(x,key = x.get) for x in profile_matrix] print("".join(consensus)) #按格式要求打印出profile,使用了print()函数不自动加换行符 for base in "ACGT": print(base+":", end='') for x in profile_matrix: print(str(x[base]), end='') print()
上面*的用法详见此贴的第四种,神操作啊。
至于信息量很大的那一行,首先是用map()函数对zip之后的序列进行了一个操作(reminder: zip之后iterable的每个item已经是列了哦),这是一个什么操作呢?这个操作/函数是自己写的lambda函数,具体内容是要数每个zip 之后的iterable item/列/x 里面base的数目(x.count);base是啥呢,请遍历"ACGT",也就是依次数这四个碱基啦,数好之后放到字典里dict(),key是base,value是count. 最后把map函数的运行结果用一个list存起来,就可以反复调用啦。希望解释明白了。。。
这就是买家秀和卖家秀的差距。。。用zip()函数变行为列是真正实现字符matrix的时刻,不用依赖pandas了呢。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-20 19:47
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社