"NotImplementedError: SeqRecord" bei der Sortierung nach einer Fasta-Datei, die mit SeqIO

Question

analysiert wurde Ich versuche, eine Fasta-Datei nach der alphabetischen Reihenfolge der Sequenzen in der Datei zu sortieren (nicht die ID der Sequenzen). Die Fasta-Datei enthält mehr als 200 Sequenzen und ich versuche, Duplikate zu finden (mit Duplikaten meine ich fast dieselbe Proteinsequenz, aber nicht dieselbe ID) innerhalb eines Bit-Masters (mit einem Python-Code). Also wollte ich ein Wörterbuch aus der Fasta-Datei machen und dann die Wörterbuchwerte sortieren. Der Code, den ich zu verwenden Ich versuche ist die folgende:

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()) 

Ich halte diese Fehlermeldung bekommen:

"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." 

Ich habe auch versucht zu suchen, wie man diesen Fehler fin aber es ares't viel Informationen darüber, und einige der Informationen, die ich gelesen habe, wo scheinbar gesagt wird, dass die Länge der im Wörterbuch gespeicherten Sequenzen ein Problem sein könnte? ... Wenn ja, wie wird die Fasta-Datei ohne SeqIO sortiert?

Answer 1

Wie mata sagte, müssen Sie eine Schlüsselfunktion zu sorted weitergeben müssen:

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) 

kehrt:

seq3 ABCDEFG 
seq0 ABCWYXO 
seq2 BCDEFGH 
seq1 IJKLMNOP 

Nun, für das, was Sie erreichen wollen. Wenn Sie die 200 Sequenzen alphabetisch sortiert haben, müssen Sie die Liste noch manuell scannen, um die Beinahe-Duplikate zu finden. Das ist fehleranfällig, also schreiben Sie auch etwas Code dafür. In der Computerwissenschaft ist edit distance eine Art zu quantifizieren, wie ungleich zwei Strings (z. B. Wörter) zueinander sind, indem die minimale Anzahl von Operationen gezählt wird, die erforderlich sind, um einen String in den anderen zu transformieren.

Es gibt verschiedene Implementierungen dieses Algorithmus. Wir nehmen die eine von 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] 

Jetzt müssen wir auf einem treshold, wie unähnlich (wie viele Einträge/delations/substituons) zwei Sequenzen sein kann entscheiden. Dann paarweise wir alle zwei Sequenzen in unserer FASTA Datei vergleichen:

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) 

Das gibt:

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

Answer 2

Hier ist eine andere Art und Weise Duplikate aus einer fasta Datei zu beseitigen basierend auf bioawk und fastq-tools (und auch awk und uniq die in der Regel in jeder UNIX-ähnlichen Befehlszeilenumgebung):

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 ">"$1"\n"$2}' 

bioawk ist eine modifizierte Version von awk, die die Manipulation einiger gängiger Bioinformatikformate erleichtert.

Die erste Zeile konvertiert die Daten in Fastq-Format, weil das fastq-sort mit arbeiten kann. Die Option -c fastx teilt bioawk mit, die Daten als fasta oder fastq zu analysieren. Dies definiert ein $name und ein $seq Feld, das in dem Befehl verwendet werden kann. Wir verwenden das Feld $seq als Dummy-Qualität, um ein gültiges fastq-Format zu erhalten.

In der zweiten Zeile wird fastq-sort (von fastq-tools) angezeigt, um die Daten nach Sequenz zu sortieren (Option -s).

Die dritte Zeile verwendet bioawk, um den Namen und die Sequenz zu extrahieren und sie als zwei durch Tabulator getrennte Felder zu setzen, damit die relevanten Informationen in einer Zeile pro Datensatz enthalten sind.

Die vierte Zeile verwendet uniq Duplikate zu eliminieren das erste Feld ignoriert, wenn aufeinanderfolgende Zeilen zu vergleichen (wenn ich richtig die Bedeutung von -f Option verstehe ich gerade entdeckt). Wenn mein Verständnis von uniq korrekt ist, wird der Name der fusionierten Datensätze der des ersten Datensatzes in einer Sequenz der gleichen Sequenz sein.

Die fünfte Zeile verwendet awk, um die durch Tabulatoren getrennten Felder in fasta neu zu formatieren.

Ich testete diesen Ansatz auf einige meiner Daten und es schien zu funktionieren.