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: }

No hay comentarios:

Publicar un comentario