Mostrando entradas con la etiqueta GD. Mostrar todas las entradas
Mostrando entradas con la etiqueta GD. Mostrar todas las entradas

28 de febrero de 2012

Leyendo el diagrama de Ramachandran

Hola,
la semana pasada les planteé a mis alumnos de la Licenciatura en Ciencias Genómicas el problema de asignar estructura secundaria a los residuos de una proteína, una vez medidos los ángulos dihedros del esqueleto peptídico psi y phi. Para esta tarea lo lógico es usar el diagrama de Ramachandran, como éste tomado de la wikipedia:


La solución naive consiste en delimitar las regiones del diagrama por medio de condiciones tales como 'si phi > X  AND psi < Y' entonces un residuo está en estado Z. Sin embargo, para tener mucha precisión sería necesario tener muchas condiciones como ésta, dada las formas irregulares del diagrama.


Una solución interesante, propuesta por F.Peñaloza, consiste en realmente leer el mapa, como lo hacemos nosotros. Lo que él hizo fue obtener un mapa de Ramachadran cuadrado con unos cuantos colores conocidos,

Diamagra de Ramachandran RGB de F.Peñaloza.


Diagrama anterior con los colores RGB empleados

para luego convertirlo a texto ASCII (con algo como text-image). Con este mapa en memoria, es sencillo obtener la estructura secundaria para unas coordenadas phi,psi dadas. 

Aprovechando el módulo GD de perl podemos leer la imagen directamente, como muestra el siguiente código:
 use strict;  
 use GD;  
   
 my $RAMAPLOT = 'ramachandran.png';  
 my $plot = GD::Image->newFromPng($RAMAPLOT);  
 my ($maxpsi,$maxphi) = checkPlot($plot);  
   
 my $inputPSI = 160;   
 my $inputPHI = -45;   
 printf("El estado de estructura secundaria para %s,%s es: %s\n",  
    $inputPSI,$inputPHI, getSimpleSS($maxpsi,$maxphi, $inputPSI,$inputPHI));  
   
 sub getSimpleSS  
 {  
   my ($maxpsi,$maxphi,$psi,$phi) = @_;  
   my %colorSS = (  
     '236,30,36','E',  
     '164,74,164','E',  
     '156,218,236','H',  
     '4,162,236','H',  
     #'180,230,28','L'   
     );  
   
   if($psi > 180 || $psi < -180 || $phi > 180 || $phi < -180)  
   {  
     die "# valid psi/phi are [-180,+180]\n";  
   }  
   $psi = (0.5 - (0.5 * $psi/180)) * $maxpsi;  
   $phi = (0.5 + (0.5 * $phi/180)) * $maxphi;  
     
   my @RGB = $plot->rgb( $plot->getPixel($phi,$psi) );  
   return $colorSS{"$RGB[0],$RGB[1],$RGB[2]"} || 'C';  
 }  
   
 sub checkPlot  
 {  
   my ($plot,$print_pixels) = @_;  
   my ($maxphi,$maxpsi) = $plot->getBounds();  
   if(!$print_pixels){ return ($maxphi,$maxpsi) }  
     
   foreach my $phi ( 0 .. $maxphi )   
   {  
     foreach my $psi ( 0 .. $maxpsi )   
    {  
       my @RGB = $plot->rgb( $plot->getPixel($phi,$psi));  
       print "$RGB[0] $RGB[1] $RGB[2]\n";  
    }  
   }  
   return ($maxphi,$maxpsi);  
 }  

Un saludo, Bruno


2 de junio de 2010

Gráfico sencillo de dominios sobre secuencia (glifo)

El problema que nos ocupa hoy es el de representar gráficamente dominios o regiones sobre una secuencia de referencia. Quizás el ejemplo más típico sea el de mapear los resultados de una búsqueda BLAST sobre la secuencia de partida, que no siempre alinean contra la secuencia completa.
Por ejemplo, dada una secuencia de 935 aminoácidos, BLAST me devuelve un alineamiento que solapa desde la posición 23 a la 512. Cómo podemos representar esto gŕaficamente de manera sencilla? Es posible usar BioPerl para este fin, como se explica en el taller de (bio)perl, pero si queremos un diagrama como este,


que representa la longitud total de la secuencia y marca el segmento alineado a escala, es mejor recurrir al módulo GD, como se muestra en este programa:

1:  #!/usr/bin/perl -w
2: # prog2.pl: crea barra PNG con uno o más dominios sobre secuencia de referencia
3:
4: use strict;
5: use GD;
6:
7: # variables globales que definen el tamaño de la figura rectangular
8: my $ANCHOFIGURA = 100;
9: my $ALTOFIGURA = 8;
10: my @COLORDOMINIO= (150,150,150); # en 3 canales RGB
11:
12: # ejemplo: secuencia de 935 aminoácidos, con un dominio reconocido del 23 al 512
13: if(pinta_dominios_secuencia(935,'23:512','panel.png'))
14: {
15: print "# archivo de salida: panel.png\n";
16: }
17:
18: sub pinta_dominios_secuencia
19: {
20: # argumentos:
21: # 1) tamaño de secuencia de referencia
22: # 2) cadena con dominios separados por compas, como: '23:78,123:189'
23: # 3) nombre de archivo PNG de salida
24: my ($length,$domain_definitions,$outPNGfilename) = @_;
25:
26: my ($n_of_domains,$dom,@domains,$startd,$endd) = (0);
27:
28: # comprueba los dominios
29: foreach $dom (split(/,/,$domain_definitions))
30: {
31: if($dom =~ /(\d+):(\d+)/)
32: {
33: ($startd,$endd) = ($1,$2);
34: if($startd<1 || $endd > $length)
35: {
36: print "# pinta_dominios_secuencia: dominio incorrecto ($dom)\n";
37: next;
38: }
39:
40: push(@domains,[$startd,$endd]);
41: $n_of_domains++;
42: }
43: else{ print "# pinta_dominios_secuencia: dominio incorrecto ($dom)\n"; }
44: }
45:
46: if(@domains)
47: {
48: my $imagen = GD::Image->new($ANCHOFIGURA,$ALTOFIGURA);
49:
50: # define colores básicos en 3 canales RGB
51: my $blanco= $imagen->colorAllocate(255,255,255);
52: my $negro = $imagen->colorAllocate(0,0,0);
53: my $color = $imagen->colorAllocate(@COLORDOMINIO);
54:
55: # crea panel rectangular con dimensiones globales $ANCHOFIGURA y $ALTOFIGURA
56: $imagen->filledRectangle(0, 0, $ANCHOFIGURA, $ALTOFIGURA, $blanco); # fondo
57: $imagen->rectangle(0, 0, $ANCHOFIGURA-1, $ALTOFIGURA-1, $negro); # borde
58:
59: # añade dominios
60: foreach $dom (@domains)
61: {
62: # calcula dimensiones del dominio a escala, eje X
63: $startd = int( $ANCHOFIGURA * ( $dom->[0] / $length ) );
64: if($startd < 1){ $startd = 1 } # evita pisar el borde
65: $endd = int( $ANCHOFIGURA * ( $dom->[1] / $length ) );
66: if($endd > $ANCHOFIGURA-2){ $endd = $ANCHOFIGURA-2 } # evita pisar el borde
67: # con 1 y -2 evitamos el borde en eje Y
68: $imagen->filledRectangle($startd, 1, $endd, $ALTOFIGURA-2, $color);
69: }
70:
71: open(PNG,">$outPNGfilename")
|| die "# pinta_dominios_secuencia: no puedo crear $outPNGfilename\n";
72: binmode PNG;
73: print PNG $imagen->png();
74: close(PNG);
75: }
76:
77: return $n_of_domains;
78: }