Metro en man: medindo un PDB

Xa escribín antes sobre os PDBs, hoxe toca algo máis metido co formato PDB. Un PDB trae unha lista de átomos con un montón de datos sobre a estructura en xeral e cada átomo en particular. Para esta entrada só nos interesan as liñas de átomos, que son as que comezan con "ATOM", "HETATM" ou "ANISOU".

Parseando lixeiramente o archivo

Buscamos estas liñas:

ATOM     10 1HG1 ILE A   1       1.867  30.218  56.564  1.00  0.00

Que según a especificación do PDB,

31 - 38 Son as coordenadas ortogonales de X en Angstroms.
39 - 46 Son as coordenadas ortogonales de Y en Angstroms.
47 - 54 Son as coordenadas ortogonales de Z en Angstroms.

Así que ahí vamos a cazar esas liñas e gardalas nun array.

import numpy as np

atom = ["ATOM", "HETATM", "ANISOU"]
coords = []
n_atoms = 0
for line in pdb:
    if line.split()[0] in atom:
        coords.append(
            [float(s) for s in [line[31:38], line[39:46], line[47:54]]])

#Reformateamos as posicions como unha matriz numpy
positions = np.reshape(coords, (n_atoms, 3))

positions vale algo así como:

[[ 133.927  117.426  173.587]
 [ 133.405  117.299  174.43 ]
 [ 133.622  116.764  172.902]
 ..., 
 [ -27.042   56.598    9.173]
 [ -26.066   57.535    9.316]
 [ -27.61    56.209    8.   ]]

Calculando o tamaño do PDB

Agora simplemente temos que tomar os valores máximos e mínimos en cada eixe, e teremos as medidas do PDB:

box_size = []

for col in range(3):
    box_size.append(abs(positions[:,col].max() - positions[:,col].min()))

print box_size

box_size vale algo así como [220.609, 121.049, 178.267]. As medidas en Angstroms do PDB.

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.