标签：short start python fa 序列 ## ref line

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

标签：short,start,python,fa,序列,##,ref,line
来源： https://www.cnblogs.com/liujiaxin2018/p/16585766.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。