python - "NotImplementedError: SeqRecord"在使用 SeqIO 解析的 fasta 文件上使用排序时

Question

我正在尝试按文件中序列的字母顺序(而不是序列 ID)对 fasta 文件进行排序。 fasta 文件包含超过 200 个序列，我试图在一个 bit master(使用 python 代码)中找到重复项(重复项我的意思是几乎相同的蛋白质序列，但不相同的 ID)。所以我想从 fasta 文件中制作一个字典，然后对字典的值进行排序。我尝试使用的代码如下:

from Bio import SeqIO


input_file = open("PP_Seq.fasta")    
my_dict = SeqIO.to_dict(SeqIO.parse(input_file, "fasta"))
print sorted(my_dict.values())

我不断收到此消息错误:

"Traceback (most recent call last):
  File "sort.py", line 4, in <module>
    print sorted(my_dict.values())
  File "/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/Bio/SeqRecord.py", line 730, in __lt__
    raise NotImplementedError(_NO_SEQRECORD_COMPARISON)
NotImplementedError: SeqRecord comparison is deliberately not implemented. Explicitly compare the attributes of interest."

我也试图寻找如何解决这个错误，但是关于这个的信息不多，而且我读到的信息中很少有明显说存储在字典字典中的序列长度可能是一个问题？.. . 如果是的话如何在没有SeqIO的情况下对fasta文件进行排序？

Answer 1

作为mata说，你需要传递一个关键函数给sorted:

from Bio import SeqIO
import operator
input_file = open("example.fasta")    
my_dict = SeqIO.to_dict(SeqIO.parse(input_file, "fasta"))
for r in sorted(my_dict.values(), key=operator.attrgetter('seq')):
    print r.id, str(r.seq)

返回:

seq3 ABCDEFG
seq0 ABCWYXO
seq2 BCDEFGH
seq1 IJKLMNOP

现在，对于您想要完成的事情。如果您已按字母顺序对 200 个序列进行排序，您仍然需要手动扫描列表以找到接近的重复项。这很容易出错，因此最好也为此编写一些代码。

在计算机科学领域，edit distance是一种通过计算将一个字符串转换为另一个字符串所需的最少操作数来量化两个字符串(例如，单词)彼此之间的不同程度的方法。

此算法有多种实现方式。我们将从 this answer 中取出一个.

def levenshteinDistance(s1, s2):
    if len(s1) > len(s2):
        s1, s2 = s2, s1

    distances = range(len(s1) + 1)
    for i2, c2 in enumerate(s2):
        distances_ = [i2+1]
        for i1, c1 in enumerate(s1):
            if c1 == c2:
                distances_.append(distances[i1])
            else:
                distances_.append(1 + min((distances[i1], distances[i1 + 1], distances_[-1])))
        distances = distances_
    return distances[-1]

现在我们需要确定两个序列可能有多不同(有多少插入/删除/替换)的阈值。然后我们成对比较 FASTA 文件中的每两个序列:

from Bio import SeqIO
from itertools import combinations
input_file = open("example.fasta")    

treshold = 4
records = SeqIO.parse(input_file, "fasta")
for record1, record2 in combinations(records, 2):
    edit_distance = levenshteinDistance(str(record1.seq), str(record2.seq))
    if edit_distance <= treshold:
        print "{} and {} differ in {} characters".format(record1.id, record2.id, edit_distance)

这给出:

seq0 and seq3 differ in 4 characters
seq2 and seq3 differ in 2 characters