ungapped Index für Biopython Alignments

Question

Mein erstes Mal mit biopython. Verzeih mir, wenn das eine grundlegende Frage ist.

Ich möchte Sequenzen eingeben, sie dann ausrichten und in der Lage sein, auf die Indexposition der ursprünglichen Sequenz (ungapped) und die ausgerichtete Sequenz (lückenhaft) zu verweisen.

Meine reale Welt Beispiel ist enolase (Uniprot P37869 und P0A6P9). Das Substrat, das Lysin bindet, ist der Index 392 in E. coli und 389 in B. subtilis. Wenn man eine Muskelausrichtung der beiden durchführt, ist der Index dieses Lysins in der Ausrichtung 394. Ich möchte in der Lage sein, leicht zwischen Lücken und ungapped Indices zu konvertieren.

Beispiel 1: Was ist der ausgerichtete Index von E. coli Rückstand # 392? (Antwort 394 in der ausgerichteten Sequenz).

Beispiel 2 Ich fand einen konservierten Rest in der Ausrichtung bei 394. Wo ist das in der ursprünglichen (ungeteilten) Sequenz? (Antwort 392 in E. coli und 389 in B. subtilis).

>>>cline = MuscleCommandline(input="enolase.txt", out="enolase.aln") 
>>>cline() 

>>> alignment = AlignIO.read("enolase.aln","fasta") 
>>> print(alignment[:,390:]) 
SingleLetterAlphabet() alignment with 2 rows and 45 columns 
AGQIKTGAPSRTDRVAKYNQLLRIEDQLAETAQYHGINSFYNLNK sp|P37869|ENO_BACSU 
AGQIKTGSMSRSDRVAKYNQLIRIEEALGEKAPYNGRKEIKGQA- sp|P0A6P9|ENO_ECOLI 
>>> print(alignment[:,394]) 
KK 

Answer 1

Fun Frage! Und soweit ich weiß nicht etwas eingebaut in BioPython. Hier ist, wie ich es lösen würde.

Lassen Sie sich mit Ihrem Beispiel 2. beginnen Wenn Sie Ihre zwei Dateien nehmen enolase.txt und enolase.aln mit jeweils ursprünglichen lückenlosen Sequenzen und abgeglichenen Sequenzen in FASTA-Format, dann können wir Schleife über die RV-Aufzeichnungen, die Anzahl der Lücken in Ihrem ausgerichteten zählen Sequenz und den Index des Rückstandes in der lückenlosen Sequenz berechnen:

from Bio import SeqIO, AlignIO 

def find_in_original(index, original_path, alignment_path): 
    alignment = AlignIO.read(alignment_path, 'fasta') 
    original = SeqIO.parse(original_path, 'fasta') 

    for original_record, alignment_record in zip(original, alignment): 
     alignment_seq = str(alignment_record.seq) 
     original_seq = str(original_record.seq) 

     gaps = alignment_seq[:index].count('-') 
     original_index = len(alignment_seq[:index]) - gaps 
     assert alignment_seq[index] == original_seq[original_index] 
     yield ("The conserved residue {} at location {} in the alignment can be" 
       " found at location {} in {}.".format(alignment_seq[index], 
       index, original_index, original_record.id.split('|')[-1])) 

Dies gibt als Ergebnis:

>>> for result in find_in_original(394, 'enolase.txt', 'enolase.aln'): 
...  print result 
The conserved residue K at location 394 in the alignment can be found at location 392 in ENO_ECOLI. 
The conserved residue K at location 394 in the alignment can be found at location 389 in ENO_BACSU. 

Für den Rückwärtsbetrieb betrachten wir alle möglichen Indizes in der Ausrichtung und sieht, welche gleich die lückenlose Sequenz, wenn wir die Lücken subtrahieren:

def find_in_alignment(index, organism, original_path, alignment_path): 
    alignment = AlignIO.read(alignment_path, 'fasta') 
    original = SeqIO.parse(original_path, 'fasta') 

    for original_record, alignment_record in zip(original, alignment): 
     if organism in original_record.id: 
      alignment_seq = str(alignment_record.seq) 
      original_seq = str(original_record.seq) 

      residue = original_seq[index] 
      for i in range(index, len(alignment_seq)): 
       if alignment_seq[i] == residue and \ 
        alignment_seq[:i].replace('-', '') == original_seq[:index]: 
        return ("The residue {} at location {} in {} is at location" 
          " {} in the alignment.".format(residue, index, 
          organism, i)) 

Dies als Ergebnis liefert:

>>> print find_in_alignment(392, 'ENO_ECOLI', 'enolase.txt', 'enolase.aln') 
The residue K at location 392 in ENO_ECOLI is at location 394 in the alignment. 
>>> print find_in_alignment(389, 'ENO_BACSU', ungapped_path, alignment_path) 
The residue K at location 389 in ENO_BACSU is at location 394 in the alignment.