新媒体管家
这是《所有人问生信媛》的第一期,主要是解答读者的问题,不定期更新。
如果你想向我们请教问题,欢迎发送邮件到[email protected] ,问题能否被答复以及何时答复只与回答者的心情有关。
比如说我会优先回答那些加入小密圈,以及有过打赏记录粉丝的问题。
如果你要提问,本期的问题就是模板。
徐老师:
您好!
我是上次在生信技能树上您指导我python入门的那个家伙。经过上次您的指导,感觉收获很大,感觉算是python摸到了入门。不过我现在又遇到了一个问题,想了很久没有很好的方法。我有一个blastsx的结果文件,其中从左往右数第一列是序列的ID号,然后有另外一个文件是序列文件里面有序列ID和相应的序列。请问如何根据blastx文件里的ID号找出序列文件里的对应的序列并计算出相应序列的长度,然后把对应的长度数值放到blastx文件的第二列,让每个ID后面都显示其对应序列的长度。我自己的想法是将文件读入列表中,计算出序列长度。但是不知道怎么样将列表中的内容按行写入文件
这是BLASTX文件
这是对应序列
这个世界都是你的已知条件
为每一步都找到最佳的工具,而不是纠结于一定要用某一个工具做所有事情
需求分析:其实很简单,就是给定一组ID,返回这些ID所在序列的长度是多少。
第一步,计算出序列文件,ID和对应的序列长度,命名为seq_len.txt 工具是bioawk
第二步,导入seq_len.tx 和 BLAST结果文件(blast_result.txt), 工具是Python
第三步,将seq_len.txt 构建字典。
第四步,用blast_result的ID查询字典
第一步: 在shell下,用bioawk获取每个ID及其序列长度
bioawk -c fastx '{len=length($2); print $1 "\t" len}' getorf_xiaomi_candidate_lncRNA > seq_len.txt
head -n 3 sqe_len.txt
PB.6.3|scaffold_1:110833-114871(-)|c/49810/f1p99 2763
PB.25.2|scaffold_1:269913-271715(+)|c/79312/f1p99 1718
PB.441.3|scaffold_1:7174735-7175436(+)|c/55975/f1p99 552
第二步: 导入两个序列
file1 = open('C:/Users/DELL/Desktop/请教徐老师关于python的问题/seq_len.txt','r')
file2 = open("C:/Users/DELL/Desktop/请教徐老师关于python的问题/blastx_getorf_xiaomi_candidate_lncRNA",'r')
seq_len = [line for line in file1]
blastx_res = [line for line in file2]
file1.close()
file2.close()
第三步,将seq_len.txt 构建字典
请熟练掌握字典推导式, 列表推导式
seq_dict = {line.split("\t")[0]: line.split("\t")[1] for line in seq_len }
第四步, 根据字典获取对应ID的长度,并对ID进行替换
# 建立一个列表,存放最后的结果
results = []
for res in blastx_res:
id_res = res.split("\t")[0]
length = seq_dict[id_res]
new_res = res.replace(id_res, id_res+'\t'+length.strip())
results.append(new_res)
最后,输出结果
with open("C:/Users/DELL/Desktop/请教徐老师关于python的问题/blast.txt",'w') as f:
for line in results:
f.write(line)
建议: 不要指望写那么多循环,既然是很多的步骤,就一步一步拆开来。