miércoles, 13 de julio de 2011

Conseguir a secuencia de un PDB

Problema



Temos un archivo .pdb, e queremos extraer a secuencia de aminoácidos de unha soa letra. Atopei algunha solución por internet, pero implica importar a librería Bio.PDB.Polypeptide.PPBuilder, e resultábame un método algo escuro. Así que busquei ata atopar o diccionario que traduza código de tres letras a código de unha letra.

Librerías



from Bio.SCOP.Raf import to_one_letter_code
from Bio.PDB import PDBParser


O diccionario está algo escondido, o módulo Raf (ASTRAL RAF (Rapid Access Format) Sequence Maps) é a primeira vez que me cruzo con él. Penso que debería estar referenciado en Bio.Alphabet.IUPAC un diccionario similar de traducción tres<->unha letras.

Código



Creamos un "motor" para parsear archivos PDB (p), parseamos o archivo en s e logo extraemos a secuencia dese archivo s a unha secuencia y, que traducimos con to_one_letter_code:

p = PDBParser()

s = p.get_structure("id", "/path/to/file.pdb")
y = [x.get_resname() for x in s.get_residues()]
sequence = "".join([to_one_letter_code[z] for z in y])

¡Listo! En sequence temos un string coa secuencia do PDB. Non sei si é máis lento ou máis ineficiente que o método do PPBuilder, pero desde logo é o método que eu estaba buscando: sacar a secuencia tri-letra do PDB e traducila. Haberá outras situacións nas que teñas unha secuencia tri-letra que queiras levar a unha letra sin ter que chamar polo PPBuilder. Nese caso, importando o diccionario e facéndolle:

[to_one_letter_code[z] for z in y]

Xa o teremos.

No hay comentarios: