"NotImplementedError: SeqRecord" when using sorted on a fasta file parsed using SeqIO

Question

我正在尝试按文件中序列的字母顺序（而不是序列 ID）对 fasta 文件进行排序。 fasta 文件包含超过 200 个序列，我试图在位主控（使用 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)

returns:

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]

现在我们需要确定两个序列可能有多少不同（多少insertions/delations/substituons）的阈值。然后我们成对比较 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

Answer 2

这是另一种从 fasta 文件中消除重复项的方法，该方法基于 bioawk and fastq-tools（以及 awk 和 uniq，它们通常存在于任何类 UNIX 命令行环境中):

bioawk -c fastx '{print "@"$name"\n"$seq"\n+\n"$seq}' test.fasta \
    | fastq-sort -s \
    | bioawk -c fastx '{print $name"\t"$seq}' \
    | uniq -f 1 \
    | awk '{print ">""\n"}'

bioawk 是 awk 的修改版本，它有助于处理一些常见的生物信息学格式。

第一行将数据转换为 fastq 格式，因为 fastq-sort 可以使用这种格式。 -c fastx 选项告诉 bioawk 将数据解析为 fasta 或 fastq。这定义了可以在命令中使用的 $name 和 $seq 字段。我们使用 $seq 字段作为虚拟质量来获得有效的 fastq 格式。

第二行告诉fastq-sort（来自fastq-tools）按顺序对数据进行排序（选项-s）。

第三行使用 bioawk 提取名称和序列并将它们作为两个制表符分隔的字段，以便每条记录在一行中包含相关信息。

第四行使用 uniq 来消除重复项，在比较连续的行时忽略第一个字段（如果我正确理解我刚刚发现的 -f 选项的含义）。如果我对uniq的理解是正确的，那么融合记录的名称将是同序列系列中第一条记录的名称。

第五行使用 awk 将制表符分隔的字段重新格式化为 fasta。

我在我的一些数据上测试了这种方法，它似乎有效。

"NotImplementedError: SeqRecord" when using sorted on a fasta file parsed using SeqIO

"NotImplementedError: SeqRecord" when using sorted on a fasta file parsed using SeqIO

python

bioinformatics

fasta

biopython