formato PDBx/mmCIF del Protein Data Bank

Hola,
aunque el formato PDB sea probablemente todavía el estándar de facto para intercambiar información de macromoléculas biológicas, se ha quedado pequeño para representar grandes complejos moleculares y por esa razón desde 2014 el formato oficial del Protein Data Bank es el PDBx/mmCIF. Hay además otro formato en liza, el PDBML, pero tiene la desventaja de ocupar mucho más espacio.

Complejo proteína-ADN 1LFU, tomado de http://www.rcsb.org/pdb/explore/explore.do?structureId=1lfu

Repasemos cómo pasan a ser los campos ATOM de un fichero PDB, los más relevantes. Veamos el siguiente ejemplo en el vetusto formato PDB, tomado de la estructura 1LFU, con columnas de ancho fijo:

ATOM      1  O5'  DA C 100      31.258  -2.296  76.212  1.00 81.62           O  
ATOM      2  C5'  DA C 100      29.867  -2.121  76.367  1.00 69.89           C  
ATOM      3  C4'  DA C 100      28.980  -3.049  77.172  1.00 67.21           C  
ATOM      4  O4'  DA C 100      29.376  -3.145  78.557  1.00 64.58           O  
ATOM      5  C3'  DA C 100      27.626  -2.376  77.196  1.00 64.41           C  
ATOM      6  O3'  DA C 100      26.569  -3.309  77.165  1.00 66.18           O  
ATOM      7  C2'  DA C 100      27.647  -1.527  78.451  1.00 63.85           C  
ATOM      8  C1'  DA C 100      28.739  -2.123  79.322  1.00 56.01           C  
ATOM      9  N9   DA C 100      29.771  -1.142  79.635  1.00 49.13           N  
ATOM     10  C8   DA C 100      30.533  -0.428  78.740  1.00 48.58           C  
ATOM     11  N7   DA C 100      31.429   0.348  79.306  1.00 43.14           N  
ATOM     12  C5   DA C 100      31.218   0.141  80.664  1.00 40.35           C  
ATOM     13  C6   DA C 100      31.837   0.679  81.794  1.00 42.42           C  
ATOM     14  N6   DA C 100      32.826   1.571  81.750  1.00 48.24           N  
ATOM     15  N1   DA C 100      31.393   0.262  82.998  1.00 42.81           N  
ATOM     16  C2   DA C 100      30.397  -0.626  83.046  1.00 45.83           C  
ATOM     17  N3   DA C 100      29.734  -1.195  82.054  1.00 41.14           N  
ATOM     18  C4   DA C 100      30.197  -0.765  80.875  1.00 41.62           C  
ATOM     19  P    DG C 101      25.116  -2.785  76.764  1.00 78.33           P  
ATOM     20  OP1  DG C 101      24.278  -3.977  76.501  1.00 74.57           O

En formato PDBx/mmCIF esto mismo se expresa así, con columnas separadas con espacios en blanco y cabeceras auto-explicativas, incluyendo una para la carga:

loop_
_atom_site.group_PDB 
_atom_site.id 
_atom_site.type_symbol 
_atom_site.label_atom_id 
_atom_site.label_alt_id 
_atom_site.label_comp_id 
_atom_site.label_asym_id 
_atom_site.label_entity_id 
_atom_site.label_seq_id 
_atom_site.pdbx_PDB_ins_code 
_atom_site.Cartn_x 
_atom_site.Cartn_y 
_atom_site.Cartn_z 
_atom_site.occupancy 
_atom_site.B_iso_or_equiv 
_atom_site.pdbx_formal_charge 
_atom_site.auth_seq_id 
_atom_site.auth_comp_id 
_atom_site.auth_asym_id 
_atom_site.auth_atom_id 
_atom_site.pdbx_PDB_model_num 
ATOM   1    O  "O5'" . DA  A 1 1  ? 31.258 -2.296 76.212  1.00 81.62  ? 100 DA  C "O5'" 1 
ATOM   2    C  "C5'" . DA  A 1 1  ? 29.867 -2.121 76.367  1.00 69.89  ? 100 DA  C "C5'" 1 
ATOM   3    C  "C4'" . DA  A 1 1  ? 28.980 -3.049 77.172  1.00 67.21  ? 100 DA  C "C4'" 1 
ATOM   4    O  "O4'" . DA  A 1 1  ? 29.376 -3.145 78.557  1.00 64.58  ? 100 DA  C "O4'" 1 
ATOM   5    C  "C3'" . DA  A 1 1  ? 27.626 -2.376 77.196  1.00 64.41  ? 100 DA  C "C3'" 1 
ATOM   6    O  "O3'" . DA  A 1 1  ? 26.569 -3.309 77.165  1.00 66.18  ? 100 DA  C "O3'" 1 
ATOM   7    C  "C2'" . DA  A 1 1  ? 27.647 -1.527 78.451  1.00 63.85  ? 100 DA  C "C2'" 1 
ATOM   8    C  "C1'" . DA  A 1 1  ? 28.739 -2.123 79.322  1.00 56.01  ? 100 DA  C "C1'" 1 
ATOM   9    N  N9    . DA  A 1 1  ? 29.771 -1.142 79.635  1.00 49.13  ? 100 DA  C N9    1 
ATOM   10   C  C8    . DA  A 1 1  ? 30.533 -0.428 78.740  1.00 48.58  ? 100 DA  C C8    1 
ATOM   11   N  N7    . DA  A 1 1  ? 31.429 0.348  79.306  1.00 43.14  ? 100 DA  C N7    1 
ATOM   12   C  C5    . DA  A 1 1  ? 31.218 0.141  80.664  1.00 40.35  ? 100 DA  C C5    1 
ATOM   13   C  C6    . DA  A 1 1  ? 31.837 0.679  81.794  1.00 42.42  ? 100 DA  C C6    1 
ATOM   14   N  N6    . DA  A 1 1  ? 32.826 1.571  81.750  1.00 48.24  ? 100 DA  C N6    1 
ATOM   15   N  N1    . DA  A 1 1  ? 31.393 0.262  82.998  1.00 42.81  ? 100 DA  C N1    1 
ATOM   16   C  C2    . DA  A 1 1  ? 30.397 -0.626 83.046  1.00 45.83  ? 100 DA  C C2    1 
ATOM   17   N  N3    . DA  A 1 1  ? 29.734 -1.195 82.054  1.00 41.14  ? 100 DA  C N3    1 
ATOM   18   C  C4    . DA  A 1 1  ? 30.197 -0.765 80.875  1.00 41.62  ? 100 DA  C C4    1 
ATOM   19   P  P     . DG  A 1 2  ? 25.116 -2.785 76.764  1.00 78.33  ? 101 DG  C P     1 
ATOM   20   O  OP1   . DG  A 1 2  ? 24.278 -3.977 76.501  1.00 74.57  ? 101 DG  C OP1   1

La diferencia más notable, además de que las columnas no son de ancho fijo, es que cada átomo tiene ahora nuevas etiquetas (label) identificativas, además de las asignadas originalmente por los autores. En el ejemplo, la desoxiadenina 100 para a llevar el número 1. Otra diferencia notable es que en estructuras de NMR, como este caso, el número de modelo se indica también para cada átomo (entity_id).

Hasta luego,
Bruno

Jugando a leer y extraer secuencias de un fichero FASTQ comprimido con GZIP (.fq.gz)

Un fichero con extensión '.fq.gz' es un fichero en formato FASTQ comprimido en formato GZIP. Los ficheros FASTQ con datos de los nuevos métodos de secuenciación (NGS) suelen ocupar decenas de Gigabytes de espacio de disco (una cantidad indecente para la mayoría de los mortales) y comprimidos con GZIP se reducen.

Para empezar intentaremos saber cuanto espacio ocupa realmente un archivo 'fq.gz' comprimido, para ello usaremos el mismo comando 'gzip':

> gzip --list reads.fq.gz
         compressed        uncompressed  ratio uncompressed_name
        18827926034          1431825024 -1215.0% reads.fq

Parece que GZIP en vez de comprimir expande, pero no es verdad, simplemente que la opción '--list' de 'gzip' no funciona correctamente para archivos mayores de 2GB. Así que hay que recurrir a un método más lento y esperar unos minutos:

> zcat reads.fq.gz | wc --bytes
61561367168

Si queremos echar un vistazo al contenido del archivo podemos usar el comando 'less' o 'zless':

> less reads.fq.gz
> zless reads.fq.gz

Y para saber el número de secuencias que contiene simplemente hay que contar el número de líneas y dividir por 4 (le costará unos minutos):

> zcat reads.fq.gz | echo $((`wc -l`/4))
256678360

Podemos buscar una determinada secuencia en el archivo y contar cuántas veces aparece, por ej. ATGATGATG:

> zcat reads.fq.gz | awk '/ATGATGATG/ {nlines = nlines + 1} END {print nlines}'

398065

A veces nos interesará tomar un trozo del fichero para hacer pruebas, por ejemplo las primeras 1000 secuencias (o 4000 líneas):

> zcat reads.fq.gz | head -4000 > test_reads.fq
> zcat reads.fq.gz | head -4000 | gzip > test_reads.fq.gz

O extraer un rango de líneas (1000001-1000004):

> zcat reads.fq.gz | sed -n '1000001,1000004p;1000005q' > lines.fq 

También nos puede interesar dividir el fichero en varios más pequeños de por ejemplo 1 miĺlón de secuencias (o 4 millones de líneas):

> zcat reads.fq.gz | split -d -l 4000000 - reads.split
> gzip reads.split*
> rename 's/\.split(\d+)/.$1.fq/' reads.split*

Y posteriormente reunificarlos:

> zcat reads.*.fq.gz | gzip > reads.fq.gz


Finalmente os recomiendo leer un post similar en inglés:
http://darrenjw.wordpress.com/2010/11/28/introduction-to-the-processing-of-short-read-next-generation-sequencing-data/

Convertir un archivo PostScript en PDF con Perl

Normalmente convertimos los archivos PostScript en formato PDF y viceversa abriéndolos con nuestro visor de documentos (por ejemplo Okular para los que somos fans de KDE) y usando el menú de Impresión donde tenemos la opción de imprimir nuestro archivo en formato PostScript o PDF. Los usuarios un poco más expertos en linux prefieren usar el comando ps2pdf (también existe pdf2ps):

Si queremos automatizar esta tarea en nuestros scripts de Perl, existe el módulo PostScript::Convert que nos permitirá hacerlo de una forma sencilla:

 use PostScript::Convert;  
 psconvert($infile, filename => $outfile, format => 'pdf');  

Sin embargo, todos los métodos anteriores generarán PDFs recortados cuando el tamaño del PostScript a convertir no es estándar, por ello os propongo una pequeña subrutina que solucionará este problema:

 # Converts a PostScript file into a PDF  
 sub convert_ps_to_pdf{  
      my ($infile,$outfile) = @_;  
      my ($height,$weight);  
      open(IDENTIFY, "identify $infile |")|| die "# $0 : cannot run 'identify $infile'\n";  
      my $ps_properties = join('',<IDENTIFY>);  
      if ($ps_properties =~ /$infile PS (\d+)x(\d+)/){  
           $weight = $1;  
           $height = $2;  
      } else {  
           die "# convert_ps_to_pdf failed to identify PostScript dimensions.\n";  
      }  
      close IDENTIFY;  
      `ps2pdf -dDEVICEWIDTHPOINTS=$weight -dDEVICEHEIGHTPOINTS=$height $infile $outfile`;  
      return $outfile;  
 }