操作爆炸结果文件 Python

Question

我写了一个 Biopython 脚本，它给了我结果，我有一个这样的文件：

>NW_020169394.1_41 [10497-10619]|KE647364.1_346 [38084-37959]
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET

>NW_020169394.1_41 [10497-10619]|KE646921.1_20 [383-240]
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET

>NW_020169394.1_41 [10497-10619]|KE647277.1_227 [70875-70720]
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET

如何在这样的单个注释行上获得结果：

>NW_020169394.1_41 [10497-10619]|KE647364.1_346 [38084-37959] | KE646921.1_20 [383-240] | KE647277.1_227 [70875-70720] 
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET                                 

我尝试使用正则表达式，但它不起作用。感谢您的回答。

Answer 1

okkey 没有使用 'SeqIO.to_dict()' 不知道，也许其他人会使用 [EDITED] 发布答案

record_dict = SeqIO.to_dict(SeqIO.parse(fasta_file , "fasta"))

会提出一个

raise ValueError("Duplicate key '%s'" % key)
ValueError: Duplicate key 'NW_020169394.1_41'

我输入的“res.txt”文件是：

>NW_020169394.1_41 [10497-10619]|KE647364.1_346 [38084-37959]
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET

>NW_020169394.1_41 [10497-10619]|KE646921.1_20 [383-240]
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET

>NW_020169394.1_41 [10497-10619]|KE647277.1_227 [70875-70720]
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET

我的代码是：

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 28 15:57:59 2021

@author: Pietro

https://stackoverflow.com/questions/68160254/manipulate-blast-result-file-python


"""
from Bio import SeqIO

fasta_file = 'res.txt'


sequences = {}
sequences2 = {}

sequa = SeqIO.parse(fasta_file, "fasta")


count = 0
for seq_record in sequa   :
            sequences[count] = seq_record
            count +=1


for key,seq_record in sequences.items():
    countz = 0
    if seq_record.seq not in [valuez.seq for keiz,valuez in sequences2.items()]:
        sequences2[countz] = seq_record
        countz += 1
        
    else:
        for keiz,valuez in sequences2.items():
            if seq_record.seq == valuez.seq:
                new = valuez.description.replace('valuez.id', '')
                newnew = new.replace(valuez.name , '')
                sequences2[keiz].description = seq_record.description + (newnew)
            
            
        
print('########################')
print(sequences2)

with open('modified_res.text', 'w') as handle:
    
    
    for key,value in sequences2.items():
            
           print(value)
           SeqIO.write(value, handle, 'fasta')

我的输出“modified_res.text”文件是：

>NW_020169394.1_41 [10497-10619]|KE647277.1_227 [70875-70720] [10497-10619]|KE646921.1_20 [383-240] [10497-10619]|KE647364.1_346 [38084-37959]
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET

无法摆脱重复的[10497-10619] 并且订单与您的不同

对不起，我尽力了

此外，如果您的序列可以通过特定的 seqrecord 格式解析（即 genbank 等，这里不是专家，可能会以不同的方式使用它们的描述，请告诉我）。

Answer 2

在 FASTA 标签中，第一个空格之前的所有内容都是 ID，并且应该是唯一的。在您的示例中并非如此，因此 SeqIO.to_dict() 将不起作用。相反，将序列映射回它们的标签，然后将它们组合起来：

from Bio import SeqIO
from collections import defaultdict

seq2label = defaultdict(list)
for record in SeqIO.parse('result.fa', 'fasta'):
    seq2label[str(record.seq)].append(record.description)

for sequence, labels in seq2label.items():
    combined_label = ' | '.join(labels[:1] + [label.split('|')[1] for label in labels[1:]])
    print(f'>{combined_label}\n{sequence}\n')

输出：

>NW_020169394.1_41 [10497-10619]|KE647364.1_346 [38084-37959] | KE646921.1_20 [383-240] | KE647277.1_227 [70875-70720]
MDQLSRKLNLTYLKVGILTSQNEFVTKHLLIIKGLKIFTET