文章来源:https://www.jianshu.com/p/2475c3240a67
简化的短序列匹配程序 (map.py) 把short.fa中的序列比对到ref.fa, 输出短序列匹配到ref.fa文件中哪些序列的哪些位置。
f1 = r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\short.fa' f2 = r'E:\Bioinformatics\Python\practice\chentong\notebook-master\data\ref.fa' ## 读入数据 #通过生成两个字典的方式进行查找 #short字典中,基因名为去除'>'及'\n'后,剩余部分 #ref字典中,基因名为去除'>'及'\n'后,剩余部分 short = {} ref = {} for line in open(f1): if line.startswith('>'): key = line.strip('>\n') short[key] = [] else: short[key] = line.strip() ## 将short保存为字典 #----end reading f1------------------- for line in open(f2): if line.startswith('>'): key = line.strip('>\n') ref[key] = [] else: ref[key].append(line.strip()) ## 将ref保存为字典 #----end reading f2(ref)-------------- #以单个ref为参照,对所有待查找序列进行遍历 for key2, value2 in ref.items(): ## 将ref作为外层迭代 #将ref中的序列进行连接,合并为一条长序列 seqRef = ''.join(value2) ## 将ref的每一个scafflod合并为一个长的字符串 for key1, value1 in short.items(): start = seqRef.find(value1) ## 根据字符串.find()返回匹配的首字符索引 while start != -1: #表明ref中可以查找到short序列 print('{}\t{}\t{}\t{}'.format(key2, start + 1, start + len(value1), value1)) ## 输出结果 new = seqRef[start+1:].find(value1) #继续在剩余序列中查找 ## 更新匹配的起始位置 if new == -1: break start = start + new + 1 #若new不等于-1,重新对start赋值(继续查找后续序列,一个循环能够对目标序列查找两遍)