python 中实现 把short.fa中的序列比对到ref.fa

 

文章来源：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赋值（继续查找后续序列，一个循环能够对目标序列查找两遍）