Bruno Contreras-Moreira
Centro de Ciencias Genómicas, UNAM
México
(versión actualizada en www.eead.csic.es/compbio/material/bioinfoPerl)
Este curso es una introducción a las aplicaciones del lenguaje Perl y su entorno en la bioinformática, el campo donde se aplican herramientas computacionales a problemas biológicos. Es un curso para alumnos que tengan ya algunas nociones de programación, de algoritmos y de manejo de sistemas UNIX*, por lo que apenas se discuten conceptos, más bien se muestran ejemplos y aplicaciones prácticas.
Este curso es parte de la materia "Computación II" de la licenciatura de Ciencias Genómicas de la Universidad Nacional Autónoma de México y cubre 5 semanas. Está escrito principalmente en español, aunque igual encontráis fragmentos en inglés repartidos entre las páginas. De todas maneras, como muchos otros lenguajes de programación, Perl está inspirado en la lengua inglesa, así que un poco de inglés viene muy bien para dominarlo.
Aquí encontraréis este curso en un sólo documento HTML, fácil de guardar o imprimir.
Podéis mandarme comentarios o sugerencias através de la
dirección de correo contrera_arroba_ccg.unam.mx.
Perl en español
Perl en inglés
Bioperl
Interfaz DBI
Bases biológicas y experimentales de la bioinformática Interfaz DBI
Búscadores específicos para código fuente
Perl es un acrónimo en inglés para Practical Extraction and Report Language, pero se ha convertido en un lenguaje con múltiples aplicaciones. Podemos consultar perl.org para hacernos una idea sobre las características del lenguaje Perl:
Perl Facts
Supported Operating Systems
Perl Features
Claro está que Perl también tiene defectos y enemigos. A mi me parece que sus virtudes son:
Por otro lado, si no se adoptan ciertas normas de estilo, se pueden escribir en Perl programas difíciles de depurar, con resultados impredecibles o imposibles de comprender. Pero estos problemas son normalmente responsabilidad del programador.
En cuanto a velocidad de ejecución, Perl tiene las desventajas de los lenguajes interpretados, y es por tanto más lento que por ejemplo C. Pero en realidad la velocidad de ejecución no es tan importante en la mayoría de las aplicaciones en bioinformática.
Espero que con este curso aprendáis a usar Perl sin estos problemas.
Para instalar Perl en tu sistema operativo, os sugiero que consultéis las instrucciones dadas en perl.com para cada sistema operativo. Debéis saber que al instalar Perl estáis instalando un programa intérprete, cuya misión es convertir el código Perl de un programador en un código ejecutable por vuestro procesador. El código fuente del intérprete está escrito en C y es público, aquí podéis ver su función principal.
Perl es un lenguaje interpretado, que deber ser compilado cada vez que se desea ejecutar un programa. Esto lo diferencia de lenguajes como C, que se compilan una sola vez y luego se pueden ejecutar N veces. Por esta razón a lenguajes como Perl en inglés se les llama scripting languages y a los programas escritos en Perl se les llama scripts.
Durante la conversión o compilación, el intérprete va revisando la sintaxis del código que le ha pasado el usuario y va informando de los errores que vaya encontrando. Esos mensajes de error son muy útiles para guiaros a lo largo del proceso de programación y prueba de vuestro código, normalmente es imprescindible leerlos.
En este curso programaremos Perl desde un terminal de un sistema UNIX* (incluyendo Linux, MacOS X, BSD y derivados), que sería equivalente a programar bajo Windows con el intérprete MS-DOS.
Además necesitaremos un editor de texto, para ir guardando en archivos el código que vayamos generando. Por ejemplo, bajo Linux yo uso Nedit o vi y bajo Windows WordPad.
Primero comprobemos qué versión de Perl tenemos. Para ellos escribimos en el terminal:
$ perl -v
Que a mi me devuelve:
This is perl, v5.8.5 built for i386-linux-thread-multi (with 1 registered patch, see perl -V for more detail) Copyright 1987-2004, Larry Wall Perl may be copied only under the terms of either the Artistic License or the GNU General Public License, which may be found in the Perl 5 source kit. Complete documentation for Perl, including FAQ lists, should be found on this system using `man perl' or `perldoc perl'. If you have access to the Internet, point your browser at http://www.perl.com/, the Perl Home Page.
Además podéis recordar con qué opciones se puede invocar el intérprete si escribís:
$ perl -h
Las opciones más importantes podrían ser:
Usage: perl [switches] [--] [programfile] [arguments] ... -d[:debugger] run program under debugger -e program one line of program (several -e's allowed, omit programfile) -v print version, subversion (includes VERY IMPORTANT perl info) -w enable many useful warnings (RECOMMENDED) -W enable all warnings ...
El proceso de creación de un programa Perl comienza normalmente abriendo un archivo de texto en un editor y escribiendo algún código. Por ejemplo, si creamos un archivo llamado programa1.pl y dentro escribimos:
print "Hola, este es mi primer programa, por supuesto con errores
Ahora podemos guardar el archivo e invocar al intérprete:
$ perl programa1.pl
La respuesta inmediata es:
Can't find string terminator '"' anywhere before EOF at test.pl line 1.
Que nos recuerda que se nos olvidó cerrar la cadena de caracteres impresa con unas comillas. Si corregimos el código, terminando la cadena con una terminador de línea y cerrando la frase con un punto y coma, tenemos:
print "Hola, este es mi primer programa sin errores\n";
Y al ejecutarlo con el intérprete, activando los avisos, obtenemos:
$ perl -w programa1.pl Hola, este es mi primer programa sin errores
En sistemas UNIX* podemos además averiguar la ruta absoluta (path) del intérprete con el comando which:
$ which perl
Que normalmente devuelve:
/usr/bin/perl
Si ahora modificamos programa1.pl:
#!/usr/bin/perl -w # ruta absoluta de perl, esto es un comentario print "Hola, este es mi primer programa sin errores\n";
Y damos al archivo privilegios de ejecución:
$ chmod +x programa1.pl
Podemos ahora ejecutar este programa directamente:
$ ./programa1.pl
Antes de profundizar en el lenguaje Perl creo importante que discutamos aspectos de lo que podríamos llamar un estilo de programación . Básicamente esto significa que al escribir un programa debemos tener en cuenta si nos podrá ser util otro día y si otras personas podrán utilizarlo y modificarlo. Si la respuesta a alguna de estas preguntas es sí, entonces se hace necesario comentar el código para que quede explicado a lo largo del código qué va haciendo cada segmento de programa.
Como veremos más adelante, otra manera de facilitar la comprensión de un programa es aislar partes del código en procedimientos y funciones. De hecho ésta es toda una estrategia de programación. Podemos escribir en lenguaje natural, por ejemplo como comentarios, las tareas que queremos haga nuestro programa y poco a poco ir sustituyendo esas descripciones de tareas por funciones y procedimientos que en efecto las ejecuten.
Finalmente, dado que el intérprete Perl puede ser demasiado flexible con la sintaxis de los programas que escribimos, os recomiendo que uséis el módulo strict , que, por ejemplo, no os dejara usar variables no declaradas o sin inicializar. De esta manera no obtendremos resultados inexplicables o inesperados en nuestros programas.
Podemos volver a modificar programa1.pl con estas nociones de estilo:
#!/usr/bin/perl -w # Programa programa1.pl, escrito por Bruno Contreras en Julio de 2005 # Este programa simplemente imprime un mensaje por la salida estándar use strict; # usa módulo strict, intérprete muy estricto con la sintaxis print "Hola, este es mi primer programa sin errores\n";
Como con cualquier lenguaje, debemos conocer al menos parte de la sintaxis de Perl para escribir programas que el intérprete pueda entender. Incluso aunque ya tengas experiencia en Perl, se te pueden olvidar a veces detalles de sintaxis. Como en la guía de referencia de Perl de Rex Swain, podemos dividir la sintaxis en partes. En este curso lo haremos de esta manera:
Vamos a ver qué son estas cosas en el programa sintaxis.pl:
#!/usr/bin/perl -w
# Programa sintaxis.pl, escrito por Bruno Contreras en Julio de 2005
# elementos de la sintaxis de Perl
# llamada al módulo strict
use strict;
# declaramos la variable $fecha
my $fecha;
# llamada al sistema operativo para ejecutar el programa date y asignar el resultado a $fecha con el operador =
$fecha = `date`;
# estructura de control if-else para ejecutar partes del programa de forma condicional
# la condición es que dentro del contenido de $fecha encontremos la subcadena
# "mar" seguida de uno ó más espacios \s+ y un dígito \d
# esto es una expresión regular, fácilmente identificable por / /
if($fecha =~ /mar\s+\d/)
{
# llamada a la función print con una cadena de caracteres como parámetro
# normalmente las funciones reciben parámetros separados por comas entre paréntesis
print("Estamos en Marzo, el mes de mi cumpleaños\n");
}
else
{
# pero perl no siempre requiere los paréntesis en las llamadas a funciones
print "Todavía queda para mi cumpleaños\n";
}
Perl tiene básicamente dos tipos de variables para contener datos: las que contienen un sólo valor y las que pueden contener más de uno.
$escalares
(scalars), variables que contienen un sólo valor
o dato. Pueden ser cadenas
de caracteres, números enteros, números reales e incluso direcciones de memoria de otras variables.
@arreglos
(arrays), matrices N-dimensionales (donde 
) que contienen más
de un escalar, que no tienen por qué ser del mismo tipo. Permiten accceder a cada valor
por medio de subíndices o coordenadas en N. El primer valor de un arreglo tiene siempre subíndice 0.
%tablas_asociativas
(hashes), que contienen parejas llave-valor. Estas tablas permiten
el acceso a los valores contenidos por medio de las llaves, que a su vez son escalares, por lo que
sólo pueden contener un valor para cada llave usada. Si se introducen dos parejas con la misma llave,
en la práctica sólo la última en asignarse es almacenada.
Un caso algo especial son las referencias, las direcciones en la memoria del computador que ocupan las variables, dadas en hexadecimal por comodidad.
A pesar de esta flexibilidad, el intérprete sí lleva un control de contexto sobre las variables y protestará si intentas aplicar, por ejemplo, una función de arreglos sobre un escalar.
En ocasiones el intérprete no puede decidir qué tipo de datos está manejando y es necesario que el
usuario diga explícitamente en qué contexto deben evaluarse, usando
${},@{},%{}
en torno a la variable en cuestión.
Por cierto, las variables que declaréis en Perl deben comenzar por una letra o un subrayado, como en
$_variable.
my $escalar = 1234; # declaración de $escalar y asignación de un entero
$escalar = 1_234; # el mismo entero, marcando el millar por comodidad
$escalar = 123.4; # real
$escalar = 0xff; # hexadecimal
$escalar = 'abc'; # cadena de caracteres
$escalar = "abc\n"; # cadena de caracteres que interpreta caracteres especiales,
# como el de nueva línea \n o tabuladores \t
$escalar = \$escalar; # referencia, dirección de $escalar
$escalar = `un_comando_sistema_operativo`; # salida de un comando ejecutado en el sistema
my @arreglo = (); # declaración de un arreglo vacío
my (@arreglo,$escalar) = ((),0); # declaración de un arreglo vacío y un escalar con valor 0
@arreglo = (1,2,3); # arreglo de enteros
@arreglo = (1,2,'tres'); # arreglo de enteros y cadenas
@arreglo = (1..25); # arreglo con los 25 primeros naturales
$arreglo[12]; # elemento número 13 de @arreglo
$arreglo[12] = 'escalar'; # asigna un escalar al elemento número 13 de @arreglo
@arreglo =[(1,2),(3,4)]; # arreglo con dos dimensiones 2x2
$arreglo[0][1]; # elemento con coordenadas 0,1 de la matriz bidimensional @arreglo
@arreglo = @copia; # copia de arreglos
$escalar = \@arreglo; # referencia, dirección de @arreglo
$escalar->[12]; # elemento número 13 de @arreglo
$$escalar[12]; # lo mismo, elemento número 13 de @arreglo
@{$escalar}; # es lo mismo que @arreglo
$arreglo[$#arreglo]; # úlimo elemento de @arreglo
@arreglo[3,4,5]; # subarreglo con 3 elementos
scalar(@arreglo); # en contexto escalar, el número de elementos del arreglo
my %tabla = (); # declaración de un tabla vacía
%tabla = ('A'=>'Alanina', 'C'=>'Cisteína'); # asigna dos parejas llave=>valor a tabla
%tabla = ('A','Alanina','C','Cisteína'); # lo mismo, asigna dos parejas llave=>valor a tabla
$tabla{'A'}; # valor de la pareja cuya llave es 'A', en este caso 'Alanina'
$tabla{'L'} = 'Leucina'; # asigna una nueva pareja a la tabla, con llave 'L' y valor 'Leucina'
$escalar = \%tabla; # referencia, dirección de %tabla
$escalar->{'L'}; # valor de la pareja cuya llave es 'L'
%{$escalar}; # lo mismo que %tabla
En Perl hay un conjunto de variables globales, accesibles al programas en ejecución, que debemos conocer para aprovechar mejor las capacidades del lenguaje. El listado completo podéis verlo en la guía de referencia de Perl, pero con que conozcáis éstas podréis ir programando:
@ARGV # guarda los argumentos con que se ejecuta un programa perl
Si invocamos en el terminal
$ perl programa.pl -f archivo.txt,
mientras dure la ejecución del programa
@ARGV
contendrá:
$ARGV[0] contiene '-f' $ARGV[1] contiene 'archivo.txt'
@_ # hace el papel de @ARGV al invocar una subrutina, función o procedimiento $_ # al leer un archivo desde un programa Perl, cada línea leída se guarda en $_
Si dentro de un programa invocamos a la función
print("Hola\n"),
la función recogerá el parámetro
"Hola\n"
como el primer elemento de
@_,
es decir,
$_[0].
Estudia el siguiente código en Perl y encuentra qué contendrá
$resultado
al terminar.
#!/usr/bin/perl -w
# Ejercicio escrito por Bruno Contreras en Julio de 2005
use strict;
my ($dato1,@datos2,%datos3,$resultado);
$datos2[0] = 'A';
$datos2[1] = 'G';
$datos2[2] = 'T';
$datos2[3] = 'C';
$dato1 = scalar(@datos2);
$datos3{ $datos2[0] } = 'ATP';
$datos3{ $datos2[1] } = 'GTP';
$datos3{ $datos2[2] } = 'TTP';
$datos3{ $datos2[3] } = 'CTP';
$resultado = \%datos3;
$resultado = $resultado->{ $datos2[$dato1-2] };
print("resultado = $resultado\n");
El repertorio completo de operadores de Perl es amplio y podéis consultarlo en la guía de referencia de Perl. Yo os menciono los que más he usado. Se sobreentiende que son operadores binarios, que relacionan dos variables, a menos que se indique lo contrario:
= # asignación
** # potencia
+ - * / # suma, resta, multiplicación y división
% # resto de la división dividendo%divisor
|| && # OR y AND, disyunción y conjunción lógica
. # concatenación de dos cadenas
-> # ya lo hemos usado para las referenciar elementos de arreglos y tablas
\ # toma la dirección de una variable (unario)
! # negación lógica (unario)
++ -- # autoincremento y autodecremento de escalares numéricos (unario)
== != # igualdad y desigualdad entre números
eq ne # igualdad y desigualdad entre cadenas (equal, not equal)
< > # mayor o menor numérico
<= >= # mayor o igual, menor o igual numérico
=~ !~ # operadores de búsqueda de patrones: comprueba que la variable izquierda contiene (no contiene)
# cierto patrón a la derecha del operador
Perl permite combinar el operador asignación con los operadores
+ , - , * , / , . ,
de forma que por ejemplo
$suma = $suma + $sumando es lo mismo que $suma += $sumando.
Como los operadores, aquí sólo mencionaré las estructuras de control de Perl que más uso, muy parecidas a las de lenguajes como C. Podéis consultar guía de referencia de Perl para ver toda la variedad.
while(condicion) # ejecuta un bloque de código mientras se cumpla cierta condición
{
# bloque de código
last if(condicion2); # termina el bucle si se da condicion2
if(condicion3) # si se da condicion3
{
# bloque de código
}
elsif(condicion4) # alternativa si se da condicion4, en caso contrario ejecuta else
{
# bloque de código alternativo
}
else
{
# otro bloque de código
}
next if(condicion5); # salta a la siguiente iteración si se da condicion5 sin ejecutar bloque4
# bloque de código 4
}
for($i=0;$i<20;$i++) # ejecuta el bloque de código 20 veces
{
# bloque de código
die if(condicion6) # termina el programa, no sólo el bucle, si se da condicion6
}
foreach $elemento (@arreglo) # bucle sobre cada uno de los elementos de un @arreglo
{
# procesa cada $elemento
}
Ahora que hemos visto variables, operadores y estructuras de control, un ejemplo más complicado.
¿Qué contiene
$tabla_amino_1_a_3 al final?
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
my @amino = ('A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V','X');
my %amino1a3 = ('A','ALA','C','CYS','D','ASP','E','GLU','F','PHE','G','GLY','H','HIS',
'I','ILE','K','LYS','L','LEU','M','MET','N','ASN' ,'P','PRO','Q','GLN',
'R','ARG','S','SER','T','THR','V','VAL','W','TRP','Y','TYR');
my $tabla_amino_1_a_3 = ""; # declara cadena vacía
my ($aa); # variable para recorrer un bucle
foreach $aa (@amino)
{
if($aa eq 'Z' || $aa eq 'B') # B y Z son formas ionizadas de asp y glu
{
next;
}
elsif($aa eq 'X')
{
$tabla_amino_1_a_3 .= "$aa equivale a cualquier aminoácido\n";
}
else
{
$tabla_amino_1_a_3 .= "$aa equivale a $amino1a3{ $aa }\n";
}
}
#print("Tabla de nomenclatura de aminoácidos\n$tabla_amino_1_a_3");
Las expresiones regulares se usan para analizar el contenido de cadenas de caracteres por medio de patrones. Son muy útiles y por ello deberéis aprender a definir vuestros propios patrones.
Como en otras secciones, os remito a la guía de referencia de Perl para ver todo el repertorio de símbolos y que se pueden usar para construir patrones.
Un patrón está delimitado por dos barras inclinadas,
/patrón/ y los caracteres
que insertéis en él pueden tener un significado literal u otro especial cuando vaya precedido
del modificador
\. Por ejemplo,
/n/
es un patrón que coincidirá con cualquier aparición de un carácter
n en una cadena. Sin embargo, /\n/
sólo coincidirá con los caracteres de nueva línea.
Los símbolos más habituales en expresiones regulares podrían ser:
/./ # cualquier carácter excepto \n, comodín (wildcard)
/\s/ # es un espacio en blanco (space)
/\t/ # es un tabulador
/\w/ # es un carácter alfanumérico (word), incluyendo '_'
/\W/ # no es un carácter alfanumérico (word)
/\d/ # es un dígito
/\D/ # no es un dígito
/\A/ # es el principio de una cadena
/\Z/ # es el final de una cadena
/^/ # es el principio de una línea de archivo
/$/ # es el final de una línea de archivo
/\// # es el carácter /
/[...]/ # es una clase de carácteres que hay que buscar
Los cuantificadores son:
+ una ó más veces
? cero ó una vez
* cero ó más veces
{n,m} mínimo y máximo de veces
Los paréntesis se usan:
/(\d\s\w)/ # para agrupar varios símbolos en un sólo patrón
/(A|B|C)/ # para definir un patrón formado por varias alternativas: coincide con 'A', con 'B' o con 'C'
Expresiones más complejas, que memorizan los patrones encontrados:
/(11+)/ # guarda la primera aparición de al menos dos unos seguidos
/(11+)\s+\1/ # guarda la primera aparición de al menos dos unos seguidos y mira a ver si
# hay otra igual separada por al menos un espacio en blanco
# \1 corresponde al primer (patrón), \2 al segundo y así sucesivamente
# sus apariciones se almacenan en las variables especiales $1, $2, ...
# $1, $2, ... son asignadas desde cada patrón agrupado con paréntesis
Veamos unos ejemplos:
if($cadena =~ /\s+/) { # ejecuta este código si hay uno ó más espacios en blanco en $cadena }
if($cadena =~ /\s{2,}/) { # ejecuta este código si hay al menos dos espacios en blanco en $cadena }
if($cadena =~ /(\d|\s)/) { # ejecuta este código si encuentra un dígito o un espacio en blanco en $cadena }
if($cadena =~ /[A-Za-z]/) { # ejecuta este código si encuentra cualquier letra, mayúscula o minúscula, en $cadena }
Otra posible forma de uso:
my $cadena = "34 KDa";
$cadena =~ /^(\d\s+KDa)/;
print "$1\n"; # en caso de encontrar el (patrón) en $cadena, $1 contendrá su primera aparición
El repertorio de aplicaciones de las expresiones regulares se ha merecido que hayan aparecido libros
dedicados sólo a ellas, por lo que no pretendo aquí mostraros todo el potencial que tienen. Valga de
ejemplo que se han usado incluso para comprobar si un número es
primo
mediante (1 x $numero) !~ /^(11+)\1+$/ .
Pero sí vamos a ver un par de funciones de Perl, las primeras después de
print y scalar,
que están asociadas a las expresiones regulares y son muy útiles. Ambas usan el operador
=~.
La primera es la sustitución (s), cuya sintaxis simplificada es:
$cadena =~ s/patrón/sustitución/[g|i]
donde g
es un modificador opcional para que sustituya todas las apariciones del patrón, no sólo la primera,
y i
es un modificador opcional para que letras iguales coincidan aunque estén mayúsculas o minúsculas.
Veamos un ejemplo:
my $cadena = "Letras y espacios \n"; $cadena =~ s/\s//; # patrón es /\s/ y sustitución es // , ahora $cadena contiene "Letrasy espacios \n"; $cadena =~ s/\s//g; # ahora $cadena contiene "Letrasyespacios\n";
La otra función es la traducción (tr) y su sintaxis es muy similar:
$cadena =~ tr/lista_caracteres/lista_caracteres_sustitución/
que reemplaza cada aparición de un carácter de la lista de caracteres de la izquierda por el carácter que ocupa la misma posición en la lista de la derecha. Veamos otro ejemplo:
$DNAcomplementario =~ tr/ATGC/TACG/;
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
my $RNA = "CAUACUAAGAUCGCGAUAUUAUUAGCGAUAUACGACU";
my $stop1 = 'UAA';
my $stop2 = 'UGA';
my $stop3 = 'UAG';
my ($total1,$total2,$total3,$total) = (0,0,0,0);
while( $RNA =~ /$stop1/g )
{
$total1++;
}
while( $RNA =~ /$stop2/g )
{
$total2++;
}
while( $RNA =~ /$stop3/g )
{
$total3++;
}
while( $RNA =~ /($stop1|$stop2|$stop3)/g)
{
$total++;
}
print "$total1|$total2|$total3|$total\n";
#!/usr/bin/perl -w
# Ejercicio escrito por Bruno Contreras en Julio de 2005
use strict;
my @enz_sitios = ('GGATCC','CC(A|T)GG'); # también sirve 'CC[AT]GG' o 'CC[A|T]GG'
my %enz_nombres = ('GGATCC','BamHI','CC(A|T)GG','EcoRI');
my $DNA = "ATCGCGACCAGGTTAGCGCGCGATTATTATATATCGATATCGCGGCGATATATCTCGAGAGTATATGGATCCTGGATAGCTAGATCATAGA";
# 1) obtengo la secuencia de la hebra de DNA complementaria, no olvidando invertirla con la función reverse
# 2) recorro el vector de enzimas de restricción mediante un bucle foreach, buscando sitios
# para cada enzima, mostrando y contando los sitios encontrados
# PISTA: para mostrar los sitios encontrados usar la variable especial $1
# 2.1) primero en la hebra + del DNA
# 2.2) ahora en la hebra complementaria
Al igual que pasa con otros lenguajes de programación, aprender a programar en Perl implica
aprender a usar buena parte de la biblioteca de funciones y procedimientos que el intérprete
incorpora por defecto, como la función
print, que ya hemos usado.
Vamos a hacer un repaso de algunas funciones esenciales en Perl, y como siempre os remito a la guía de referencia de Perl para que recordéis el repertorio completo cuando estéis programando.
scalar
Devuelve el tamaño en número de elementos de un arreglo.
my @arreglo = (1,2,3,"cuatro"); print scalar(@arreglo); # imprime 4
push
Ya en la sección 2.3.2 vimos cómo asignar elementos a un arreglo usando subíndices.
La función
push
permite añadir elementos o incluso arreglos al final de un arreglo dado.
push(@arreglo,"elemento"); push(@arreglo,$elemento); my @arreglo = (1,2,3); my @otro_arreglo = (4,5,6); push(@arreglo,@otro_arreglo); # ahora @arreglo contiene 1,2,3,4,5,6
unshift
Añade elementos o arreglos al principio de un arreglo dado.
my @arreglo = (1,2,3); my @otro_arreglo = (4,5,6); unshift(@arreglo,@otro_arreglo); # ahora @arreglo contiene 4,5,6,1,2,3
pop
Extrae y devuelve el último elemento de un arreglo (un escalar), que reduce en uno su tamaño.
my @arreglo = (1,2,3); my $ultimo = pop(@arreglo); # $ultimo = 3, @arreglo = (1,2)
shift
Extrae y devuelve el primer elemento de un arreglo (un escalar), que reduce en uno su tamaño.
my @arreglo = (1,2,3); my $primer = shift(@arreglo); # $primer = 1, @arreglo = (2,3)
reverse
Toma un arreglo y devuelve otro arreglo con los mismos elementos invertidos.
my @arreglo = (1,2,3); my @invertido = reverse(@arreglo); # @invertido = (3,2,1)
grep
Recibe dos parámetros, una expresión regular y un arreglo, y compara el patrón con cada
uno de los elementos del arreglo. Devuelve un arreglo con los elementos que coincidieron
con el patrón. Es muy útil para comprobar si ciertos valores están contenidos en un arreglo.
my @arreglo = ("uno","dos","tres");
my @coinc = grep(/s/,@arreglo); # @coinc = ("dos","tres")
sort
Ordena los elementos de un arreglo de menor a mayor alfabéticamente, devolviendo otro arreglo
ya ordenado. Se puede cambiar el criterio de ordenación.
my @arreglo = (3,7,1.1,11.8);
my @ordenado = sort(@arreglo); # @ordenado contiene 1.1, 11.8, 3, 7
@ordenado = sort {$a<=>$b} (@ordenado); # ahora contiene 1.1, 3, 7, 11.8, criterio numérico ascendente
@ordenado = sort {$b<=>$a} (@ordenado); # ahora contiene 11.8, 7, 3, 1.1, criterio numérico descendente
Las cadenas de caracteres en Perl, como sabéis, las manejamos como variables escalares y tienen también unas cuantas funciones específicas asociadas.
substr
Devuelve una subcadena de la cadena pasada como parámetro. Puede recibir dos o tres parámetros.
El primero es una cadena, el segundo el punto (offset, empezando en el 0) donde empieza
la nueva subcadena y el tercero, opcional, la longitud de la subcadena a partir del punto
anterior.
my $cadena = "cadena"; my $subcadena = substr($cadena,0,4); # $subcadena contiene "cade"; $subcadena = substr($cadena,1,4); # $subcadena contiene "aden"; $subcadena = substr($cadena,2); # $subcadena contiene "dena";
index
Devuelve la posición de una subcadena dentro de una cadena, el valor devuelto en caso contrario
es -1. Puede recibir un tercer parámetro opcional que marca en qué posición debe
empezar la búsqueda.
my $cadena = "cadena";
my ($subc1,$subc2) = ('en','no');
print index($cadena,$subc1); # imprime 3
print index($cadena,$subc1,4); # imprime -1
print index($cadena,$subc2); # imprime -1
length
Devuelve la longitud de una cadena de caracteres.
print length("cadena"); # imprime 6
chomp
Elimina el último terminador de línea (símbolo de nueva línea)
\n de una cadena
y devuelve 1 si de hecho lo eliminó. Se puede aplicar también a un arreglo de cadenas.
my $cadena = "cadena\n";
my @cadenas=("cadena1\n","cadena2\n");
print chomp($cadena); # imprime 1
print chomp(@cadenas); # imprime 2
lc, uc
Devuelven versiones en mayúsculas (upper case, uc) o minúsculas (lower case, lc)
de la cadena pasada como parámetro.
my $cadena = "May y min"; print lc($cadena); # imprime 'may y min' print uc($cadena); # imprime 'MAY Y MIN'
split
Función que recibe un patrón regular y una cadena y trata de extraer las diferentes subcadenas
que haya separadas por tal patrón, devolviéndolas en un arreglo. Se puede pasar un tercer
parámetro que indica el número máximo de subcadenas deseadas.
my $cadena = "cadena de caracteres|otra cadena\n"; my @subcadenas1 = split(/\s+/,$cadena); my @subcadenas2 = split(/\|/,$cadena); # separa por el carácter '|' print scalar(@subcadenas1); # contiene 4 elementos print scalar(@subcadenas2); # contiene 2 elementos
keys
Devuelve un arreglo con las llaves de la tabla asociativa pasada como parámetro. Es muy útil
para realizar bucles iterativos sobre tablas.
my $codon;
my %codones = ('GCC'=>'Ala','GCG'=>'Ala','ATG'=>'Met');
my @llaves = keys(%codones);
foreach $codon (@llaves)
{
print "$codon $codones{ $codon } \n";
}
foreach $codon (sort (@llaves)) # itera sobre las llaves ordenadas alfabéticamente
{
print "$codon $codones{ $codon } \n";
}
values
Devuelve un arreglo con los valores de la tabla asociativa pasada como parámetro.
my $codon;
my %codones = ('GCC'=>'Ala','ATG'=>'Met');
my @valores = values(%codones); # contiene 'Ala' y 'Met'
Aprovecho aquí para mostraros como se pueden ordenar las tablas asociativas por valor. Como es lógico, la solución pasa por ordenar las llaves por sus valores asociados:
my @llaves_ordenadas = sort {$hash{$a} <=> $hash{$b}} ( keys(%hash) );
Las funciones aritméticas más usuales son quizás éstas.
abs
Devuelve el valor absoluto del escalar numérico pasado como parámetro.
my $valorabs = abs($escalar);
exp
exp($exponente) devuelve la exponencial de
$exponente.
int
Devuelve la parte entera del número real pasado como parámetro.
my $valorint = int(8.2); # $valorint contiene 8
log
Devuelve el logaritmo natural del o escalar numérico.
sqrt
Devuelve la ráiz cuadrada del parámetro escalar numérico.
my $raiz = sqrt(25); # $raiz contiene 5
rand
Esta función es un generador de números aleatorios que por defecto devuelve un valor
real comprendido entre 0 y 1, aunque puede delimitarse el límite superior pasándolo como
parámetro. Funciona con un valor semilla que puedes definir llamando a la función
srand,
por ejemplo para obtener resultados replicables en tus programas aleatorios. Con la misma
semilla, dos llamadas a
rand
devolverán el mismo valor aleatorio.
my $semilla = 12345; srand($semilla); my $azar = rand(14); # 3.15459917914819 , entre 0 y 14 my $azar = rand(); # 0.919183068533556 , entre 0 y 1
Llegados a este punto creo que podemos ya empezar con algunos ejercicios que tocan de cerca problemas de la bioinformática.
Muchos de vosotros ya sabréis de la importancia de los alineamientos de secuencias en bioinformática, puesto que tienen aplicaciones en cosas tan diferentes como la predicción de genes en genomas o la predicción de estructura secundaria en proteínas.
Los algoritmos de alineamiento de secuencias, como BLAST (Altschul et al., 1997), usan matrices de sustitución para dirigir el proceso de alineamiento. Para proteínas, una de las matrices más conocidas es BLOSUM62 (Henikoff & Henikoff, 1993).
|
En este ejemplo vamos a construir una matriz de sustitución de aminoácidos, como BLOSUM62, pero basada en el código genético. La hipótesis que vamos a manejar es que los aminoácidos cuyos codones se parecen más son más fácilmente intercambiables, olvidándonos por un momento de sus propiedades bioquímicas.
#!/usr/bin/perl -w
# Ejercicio escrito por Bruno Contreras en Julio de 2005
# Programa que calcula una matriz de sustitución de aminoácidos basándose únicamente en el parecido de
# sus codones de tRNA
use strict;
my %tripletes = (
'F' => ['TTT','TTC'],
'L' => ['TTA','TTG','CTT','CTC','CTA','CTG'],
'S' => ['TCT','TCC','TCA','TCG','AGT','AGC'],
'Y' => ['TAT','TAC'],
'C' => ['TGT','TGC'],
'W' => ['TGG'],
'P' => ['CCT','CCC','CCA','CCG'],
'H' => ['CAT','CAC'],
'Q' => ['CAA','CAG'],
'R' => ['CGT','CGC','CGA','CGG','AGA','AGG'],
'I' => ['ATT','ATC','ATA'],
'M' => ['ATG'],
'T' => ['ACT','ACC','ACA','ACG'],
'N' => ['AAT','AAC'],
'K' => ['AAA','AAG'],
'V' => ['GTT','GTC','GTA','GTG'],
'A' => ['GCT','GCC','GCA','GCG'],
'D' => ['GAT','GAC'],
'E' => ['GAA','GAG'],
'G' => ['GGT','GGC','GGA','GGG'],
'X' => ['XXX'],
);
my @aa = ('A','R','N','D','C','Q','E','G','H','I','L','K','M','F','P','S','T','W','Y','V','X');
my $MAXPUNTUACION = 0; # máximo valor de la matriz de sustitución
# es una matriz de costes, con números negativos
my ($aa1,$aa2,$cod1,$cod2,@codones1,@codones2);
my ($dist,$distmedia);
print "# Matriz de sustitución de aminoácidos basada en el código genético\n";
print " A R N D C Q E G H I L K M F P S T W Y V X\n";
foreach $aa1 (@aa)
{
@codones1 = @{ $tripletes{$aa1} }; # toma todos los codondes correspondientes a $aa1
print "$aa1 ";
foreach $aa2 (@aa)
{
if($aa1 eq $aa2){ print " $MAXPUNTUACION "; }
else
{
@codones2 = @{ $tripletes{$aa2} }; # toma todos los codondes correspondientes a $aa2
$distmedia = 0;
foreach $cod1 (@codones1)
{
foreach $cod2 (@codones2)
{
# cuenta el número de cambios que hay que hacer para pasar de un codón a otro (máximo 3)
$dist = 0;
if(substr($cod1,0,1) ne substr($cod2,0,1)){ $dist++; }
if(substr($cod1,1,1) ne substr($cod2,1,1)){ $dist++; }
if(substr($cod1,2,1) ne substr($cod2,2,1)){ $dist++; }
# recuerda la distancia entre estos dos codones para calcular la media entre aa1 y aa2
$distmedia += $dist;
}
}
# calcula la distancia media entre aa1 y aa2
$distmedia /= (scalar(@codones1)*scalar(@codones2));
$distmedia = int($distmedia);
# imprime el coste medio para pasar de aa1 a aa2
print "-$distmedia ";
}
}
print "\n";
}
¿Se te ocurre cómo modificar el algoritmo?
Ahora os toca a vosotros. La tarea es que uséis la matriz de sustitución generada en el ejemplo anterior para evaluar este alineamiento en formato FASTA:
>unnamed ------KNATRESTSTLKAWLSEHRKNPYPTKGEKIMLAIITKMTLTQVSTWFANARRRLKKENKMTWTPRNRTDEEGNVYSS------ >PDB1LFU M.musculus DNA-binding protein MARRKRRNFNKQATEILNEYFYSHLSNPYPSEEAKEELAKKSGITVSQVSNWFGNKRIRYKKN-------IGKFQEEANIYAAKTAVTA
Os sugiero que modifiquéis el código del ejemplo anterior 2.7.5 para convertir la matriz impresa en una tabla asociativa donde la clave sea un par de aminoácidos y el valor su coste. Sería conveniente imprimir el alineamiento con el coste de cada pareja de aminoácidos alineados, para localizar las zonas más confiables del alineamiento.
Una función de Perl que os puede ser de utilidad es
rand,
el generador de números
aleatorios. Los números aleatorios se han usado en biología para modelar eventos de mutación
y de recombinación genética . Pero se han llevado más
allá en computación y son la base de varios tipos de algoritmos como los algoritmos genéticos
y de Monte Carlo.
En este ejemplo que os propongo vamos a implementar un algoritmo genético muy sencillo que va a recombinar una serie de secuencias de ADN para obtener una secuencia final recombinante con un contenido en GC que se aproxima a uno definido.
#!/usr/bin/perl -w # Ejercicio escrito por Bruno Contreras en Julio de 2005 # Programa que recombina una serie de moléculas de ADN que codifican para genes homólogos, # y obtiene un gen recombinante que optimiza su contenido en %GC # REQUISITOS: # los genes deben estar alineados y si tienen diferente longitud deberán incorporar inserciones # mira http://www.cbs.dtu.dk/services/RevTrans/ use strict; @genes = ( # en este ejemplo hay 3 genes homólogos del regulador CRP #CRP_HAEIN 'ATGTCAAATGAATTAACCGAAATTGATGAAGTTGTAACCTCCTCTCAAGAAGAAGCAACT CAACGAGATCCCGTTTTAGATTGGTTTCTTACTCACTGCCATTTGCATAAATATCCTGCA AAATCAACTTTAATTCATGCTGGGGAAGATGCGACCACGCTGTATTATGTAATTAAAGGT TCTGTAATGGTATCTTCAAAAGATGATGAAGGCAAAGAGATGATCCTCACTTACTTAGGT GCAGGACAATTTTTTGGCGAAGCGGGATTATTTGATGAAGGTTCAAAACGATCAGCTTGG GTAAAAACAAAAACAACATGTGAAATTGCTGAAATTTCCTATAAGAAATATCGCCAGTTG ATTCAGGCAAACCCTGAAATCTTAATGTTTCTCACTGCACAATTGGCAAGACGTTTGCAA AATACATCACGTCAAGTCACGAATTTGGCATTTTTAGACGTCGCAGGTCGCATCGCTCAA ACTTTAATGAACTTAGCTAAACAGCCTGAAGCAATGACGCATCCTGATGGTATGCAAATC AAAATTACACGCCAAGAAATAGGGCAAATGGTGGGTTGTTCACGAGAAACTGTGGGGCGC ATTATTAAGATGTTGGAGGATCAGAATCTTATCCACGCTCATGGAAAAACAATCGTTGTA TATGGCGCAAGA', #CRP_KLEAE 'ATGGTGCTTGGCAAACCGCAAACA------------------------------------ ------GACCCTACCCTTGAATGGTTCTTGTCTCATTGCCACATTCATAAGTACCCATCA AAGAGCACGCTGATCCACCAGGGTGAAAAAGCAGAAACGCTGTACTACATCGTTAAAGGC TCCGTGGCTGTACTCATCAAGGATGAAGAAGGTAAAGAGATGATCCTCTCCTACCTCAAC CAGGGCGATTTCATCGGTGAATTAGGCCTGTTTGAAGAGGGTCAGGAGCGTAGCGCCTGG GTACGGGCGAAAACCGCATGTGAAGTGGCCGAAATCTCTTATAAGAAATTCCGTCAGCTG ATCCAGGTGAACCCGGACATTCTGATGCGTCTCTCTTCGCAAATGGCTCGCCGTCTGCAG GTCACGTCTGAGAAAGTGGGCAACCTCGCCTTCCTCGACGTGACGGGCCGTATCGCCCAG ACGCTGCTGAACCTGGCGAAGCAACCGGATGCCATGACCCACCCGGACGGTATGCAAATT AAAATTACACGCCAGGAAATTGGTCAGATCGTCGGCTGCTCTCGTGAAACCGTTGGTCGT ATTTTGAAAATGCTGGAAGATCAGAACCTGATTTCCGCGCACGGTAAAACCATCGTCGTC TACGGCACCCGT', #CRP_PASMU 'ATGCAAACTACACCATCGATA--------------------------------------- ------GATCCTACGTTAGAGTGGTTTCTATCTCACTGCCATATTCACAAATATCCTTCC AAAAGTACATTGATTCATGCAGGTGAAAAAGCCGAGACGCTGTATTATCTGATTAAAGGT TCCGTTGCTGTCTTAGTCAAAGATGAAGATGGCAAAGAAATGATTTTGACTTATTTAAGT CAAGGTGATTTTTTTGGTGAGGCTGGTCTTTTCGAAGAGGGACAATTACGTTCAGCTTGG ATTAAGGCAAAAAGCCCTTGTGAAATTGCTGAAATTTCTTATAAAAAATTTCGTCAATTA ATTCAGGTCAACCCAGATATTTTAATGCATTTATCAGCTCAGCTAGCACGTCGCTTACAA AATACGTCGAGACAAGTAAGTAATCTAGCCTTTTTAGATGTGACGGGGCGTATCGCACAA ACCTTATTAAATCTGGCTAAAATGCCAGAAGCCATGACCCATCCAGATGGGATGCAAATT AAAATCACGCGCCAAGAAATTGGGCAAATGGTGGGATGTTCGCGTGAAACGGTAGGACGA ATCCTCAAAATGTTAGAAGATCAGCATTTAATCTCTGCACATGGTAAAACCATTGTCGTT TACGGTACAAGA' ); # 0) establecer el %GC deseado # 1) abrir un bucle while que recombine aleatoriamente estos genes, # conservando el marco de lectura y controle el GC de las secuencias # recombinantes que van saliendo # 2) calcula el GC con algo como # my $totalGC = $gen =~ tr/gc/gc/; # my $GC = $totalGC/length($gen); # 3) sal del bucle cuando GC lleve n generaciones sin cambiar # 4) elige al azar dos secuencias para recombinar # 5) para recombinar usa rand con la longitud del alineamiento # restringiendo los puntos de sobrecruzamiento a los que sean múltiplos # de tres
Si queremos manejar el fichero
fichero.txt
haríamos:
open(MIFICH,"fichero.txt"); # abrimos un fichero existente para lectura
# y le asignamos el descriptor MIFICH
open(MIFICH,">fichero.txt"); # creamos un nuevo fichero y le asignamos
# el descriptor MIFICH
open(MIFICH,">>fichero.txt"); # abrimos un fichero existente para añadir
# contenido al final, conservando el contenido
# original, y le asignamos el descriptor MIFICH
close(MIFICH); # así cerramos el archivo en cualquier caso
Una versión más segura hubiera sido:
open(MIFICH,"fichero.txt") || die "lo siento, no puedo encontrar fichero.txt\n";
Para leer un archivo normalmente se usa un bucle
while:
open(MIFICH,"fichero.txt")|| die "lo siento, no puedo encontrar fichero.txt\n";
while(<MIFICH>) # significa: para cada línea ejecuta el bloque de código
{ # entre los corchetes
print $_; # la variable especial $_ contiene la línea actual
# ahora podríamos usar expresiones regulares para procesar esta línea de texto
}
close(MIFICH);
Para escribir en un fichero que acabamos de abrir:
open(MIFICH,">fichero.txt")|| die "lo siento, no puedo crear fichero.txt\n";
print MIFICH "Una línea dentro del archivo\n";
close(MIFICH);
Cuando usamos la función
print
sin especificar un descriptor de archivo, como hemos hecho en varios
ejemplos hasta ahora, en realidad estamos escribiendo al archivo de salida por defecto, llamado STDOUT.
De la misma manera podemos acceder al teclado, el archivo de entrada por defecto, por medio de <STDIN>.
Por otro lado, existen en Perl operadores especiales para operar sobre ficheros (filetest operators) que pueden ser muy útiles (guía de referencia de Perl):
-r -w -x # comprueba si hay permiso para leer, escribir o ejecutar un archivo
-e -z # comprueba si el archivo existe / tiene tamaño nulo
-s # compruebe si el archivo tiene tamaño no nulo y devuelve el tamaño
-f -d # comprueba si el archivo es un fichero / es un directorio
Ejemplo de uso:
if(-r 'mifichero.txt')
{
# abro y leo el archivo
}
Hay toda una batería de funciones específicas que podéis consultar en la guía de referencia de Perl, pero creo que estas son las esenciales.
Perl también tiene funciones para manejar directorios, como se muestra en el siguiente ejemplo:
$ ls -la ruta/directorio
-rw-r--r-- 1 contrera users 52849 jul 20 19:15 archivo.txt
-rw-r--r-- 1 contrera users 715 jul 19 11:08 datos
-rw-r--r-- 1 contrera users 129 jul 7 16:49 .latex
drwxr-xr-x 2 contrera users 1024 jul 7 13:38 png
opendir(DIR,"ruta/directorio") || die "lo siento, no puedo abrir ruta/directorio\n";
my @entradas = grep {!/^\./} readdir(DIR); # añade a @entradas archivos y subdirectorios que no empiecen por '.'
# @entradas contiene: archivo.txt, datos, png
closedir(DIR);
Para crear un directorio se usa la función
mkdir 'nombre_directorio_nuevo .
El Protein Data Bank (Berman et al., 2000) es un repositorio de descripciones experimentales de las estructuras moleculares de proteínas y ácidos nucleicos resueltos hasta el momento. Cada descripción es un archivo de texto que contiene las coordenadas atómicas de la molécula en cuestión en un formato que se llama PDB.
Aquí podéis ver un ejemplo sencillo, la estructura de una proteína reguladora de ratón en complejo con una molécula de ADN, en el archivo 1lfu.pdb. Podéis desplegar esta molécula con programas como Rasmol y, en cualquier caso, debéis abrir este archivo con vuestro editor de texto para verle las tripas.
En este ejemplo es muestro el programa
extraeProtPDB.pl
que lee un archivo PDB, como 1lfu.pdb
y separa el esqueleto peptídico (backbone) de las proteínas
presentes y las deposita en un nuevo archivo.
![\begin{figure}
\begin{center}
\includegraphics[width=0.8\textwidth]{1lfu}
\end{center}
\end{figure}](img2.png)
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
# Programa extraeProtPDB.pl que lee un archivo PDB pasado como argumento,
# separa el esqueleto peptídico de las proteínas presentes y lo deposita
# en formato PDB en un archivo llamado esquelto.pdb
use strict;
## variables importantes ###########################################################
my ($archivoPDB,$nombre_atomo,@esqueleto);
## revisa los argumentos pasados al invocar el programa ############################
if(scalar(@ARGV) < 1)
{
die "uso: extraeProtPDB.pl <archivoPDB>\n";
}
else
{
$archivoPDB = $ARGV[0]; # lee el nombre del argumento y asígnalo a
}
## abre $archivoPDB en modo lectura y lo recorre línea a línea ###############
open(PDB,$archivoPDB) || die "# no puedo abrir $archivoPDB\n";
while(<PDB>)
{
if(/^ATOM/) # líneas de coordenadas atómicas, si abrís 1lfu.pdb sabréis a qué me refiero
{
$nombre_atomo = substr($_,13,4); # toma la parte de la línea que contiene el nombre del átomo
# en el formato PDB
# toma sólo los átomos correspondientes al esqueleto peptídico: N,CA,C,O
if($nombre_atomo =~ /N / || $nombre_atomo =~ /CA/ || $nombre_atomo =~ /C /|| $nombre_atomo =~ /O /)
{
push(@esqueleto,$_);
}
}
}
close(PDB);
# crea un nuevo archivo para depositar ahí las coordenadas del esqueleto peptídico ###
open(ESQPDB,">esqueleto.pdb") || die "lo siento, no puedo crear esqueleto.pdb\n";
print ESQPDB "HEADER esqueleto peptídico de $archivoPDB\n";
print ESQPDB @esqueleto;
close(ESQPDB);
Ahora, como ejercicio, os pido que modifiquéis este programa para que cree una cadena de caracteres con la secuencia de la proteína 1lfu.pdb en código de una letra. Por ejemplo, Ala deberá ser 'A' y Trp 'W'. El programa deberá además crear un nuevo archivo e imprimir dentro la secuencia en formato FASTA, como el que vísteis en 2.7.5.
Con el archivo FASTA generado en el ejercicio anterior 1lfu.fas hacemos dos búsquedas BLAST (Altschul et al., 1997) con los siguientes comandos (necesitas tener instalado BLAST localmente y bajarte la biblioteca de secuencias nr):
blastall -p blastp -i 1lfu.fas -d nr -o 1lfu.blast -v 20 -b 20 # -o es el archivo de resultados blastpgp -i 1lfu.fas -j 2 -Q 1lfu.pssm -d nr -o 1lfu.blast -v 20 -b 20 # -o y -Q son archivos de resultados
Los archivos de resultados generados son 1lfu.blast y 1lfu.blastpgp, donde se encuentra la información relativa a las secuencias potencialmente homólogas encontradas en la biblioteca nr (non redundant), que puedes obtener en el NCBI.
El último archivo de resultados generado es 1lfu.pssm, que contiene información también muy valiosa, entre la que destacamos una matriz de sustitución de aminoácidos específica para la familia de proteínas 1lfu y el contenido en información evolutiva de cada aminoácido de la secuencia de entrada 1lfu.fas.
En el siguiente ejemplo os muestro un trozo de código que extrae de la salida de BLAST (indistintamente 1lfu.blast o 1lfu.blastpgp) la información relevante de cada uno de los alineamientos reportados.
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
# Programa extraeBLAST.pl que lee un archivo generado por BLAST pasado como argumento
# y extrae la información relevante de cada alineamiento reportado
use strict;
## variables ######################################################################
my ($archivoBLAST,@hits);
# otras variables manejadas, donde Q significa 'query' y S ' subject',
# las dos secuencias de un alineamiento en BLAST
my ($aux0,$aux1,$aux2,$aux3); # variables auxiliares
my ($Q,$S); # variables booleanas que indican si se está leyendo Q o S
my ($firstQ,$lastQ,$seqQ); # primer resíduo, último resíduo y secuencia de 'query'
my ($firstS,$lastS,$seqS);
my ($longQ,$name,$id,$ev); # longitud de 'query', nombre se secuencia, % de identidad y e-value
## revisa los argumentos pasados al invocar el programa ############################
if(scalar(@ARGV) < 1)
{
die "uso: extraeBLAST.pl <archivoBLAST>\n";
}
else
{
$archivoBLAST = $ARGV[0]; # lee el nombre del argumento y asígnalo a
}
## abre $archivoBLAST en modo lectura y recórrelo línea a línea ####################
$Q = $S = 0;
open(BLAST,$archivoBLAST) || die "# no puedo abrir $archivoBLAST\n";
while(<BLAST>)
{
if(/Query:/)
{
($aux0, $aux1, $aux2, $aux3) = split; # split separa $_ en 4 escalares separados por espacios en blanco
if(!$Q){ $firstQ = $aux1; }
$lastQ = $aux3;
$seqQ .= $aux2;
$Q++; # estamos leyendo 'query'
}
elsif(/Sbjct:/)
{
($aux0, $aux1, $aux2, $aux3) = split;
if(!$S){ $firstS = $aux1; }
$lastS = $aux3;
$seqS .= $aux2;
$S++; # estamos leyendo 'subject'
}
elsif((/^>/)||(/^ Database:/)) ## procesa el alineamiento previo
{
if($Q)
{
# crea un arreglo con los datos de este alineamiento y guárdalo en @hits
# @hits es por tanto un arreglo de arreglos
push(@hits, [($firstQ,$lastQ,$firstS,$lastS,$name,$id,$ev,$seqQ,$seqS)] );
}
# reiniciamos valores para el siguiente alineamiento
$Q = $S = 0;
$seqQ = $seqS = $name = $id = $ev = "";
# nombre del siguiente alineamiento
$name = $_;
chomp($name);
}
elsif(/^ Score =/)
{
if($Q) ## secuencias que tienen varios alineamientos
{
push(@hits, [($firstQ,$lastQ,$firstS,$lastS,$name,$id,$ev,$seqQ,$seqS)] );
$Q = $S = 0;
$seqQ = $seqS = $id = "";
}
$ev = (split)[7];
}
elsif(/^ Identities/)
{
# toma nota del % de identidad del alineamiento
$id = (split)[3];
$id =~ s/(\(|\)|\%|,)//g;
}
elsif(/^Sequences/) # reinicia todo para la siguiente ronda de resultados (sólo blastpgp)
{
undef @hits;
$S = $Q = 0;
$seqQ = $seqS = $name = $id = $ev = "";
}
elsif(/letters\)/)
{
# toma nota de la longitud de la secuencia 'query'
$longQ = (split)[0];
$longQ =~ s/\(//;
}
}
close(BLAST);
## ahora como ejemplo muestra el primer alineamiento guardado #################################
# cada alineamiento se guardó así: [($firstQ,$lastQ,$firstS,$lastS,$name,$id,$ev,$seqQ,$seqS)]
# el primer alineamiento será $hits[0], que a su vez es un arreglo, por eso necesitamos otra
# coordenada para llegar a cada elemento
if($hits[0])
{
print "Posible homólogo:\n$hits[0][4]\n";
print "Fragmento de $hits[0][0] a $hits[0][1], e-value $hits[0][6] y \%ID $hits[0][5]\n";
print "Alineamiento\n$hits[0][7]\n$hits[0][8]\n";
}
Vuestra tarea es ahora escribir un programita que lea 1lfu.pssm y genere una lista de los aminoácidos cuyo contenido en información evolutiva sea elevado. Os sugiero que uséis la función split para leer el archivo línea a línea.
Cuando queremos que un programa consulte al usuario normalmente necesitaremos leer su respuesta en el archivo estándar de entrada, STDIN . Lo haremos por ejemplo así:
print "Contenido en GC deseado: ";
my $GCusuario = <STDIN>; # cuando el usuario teclee algo y presione <ENTER>
# el programa prosigue
Ya hemos visto antes que la salida estándar de datos en Perl es un descriptor de archivo que
se llama STDOUT y es a donde escribe
print por defecto.
Por último, hay una salida estándar para los errores que se produzcan en tiempo de ejecución
que se llama STDERR, a donde podemos hacer que escriban nuestros programas:
print STDERR "Error en mi programa, división por cero\n";
Normalmente STDERR se muestra en el terminal, al igual que STDOUT.
Ya hemos usado en diferentes ejemplos la función
print.
Esta función puede usarse para
dar formato a los datos que queremos imprimir en un programa. Una manera práctica muchas veces
utilizada en bioinformática es imprimir datos en un archivo por columnas, de forma que todas
las líneas tengan el mismo número de columnas. Como separadores de columnas pueden utilizarse
carácteres como
\t (tabulador), espacios en blanco o comas:
print "$nombreGen \t $GC \t $longitud \t $especie\n"; print "$nombreGen $GC $longitud $especie\n"; print "$nombreGen,$GC,$longitud,$especie\n";
Este tipo de formatos de salida son ideales para que otros programas en Perl, por medio de
la función
split,
puedan leerlos. De hecho son también adecuados para que aplicaciones como Excel u OpenOffice puedan
abrir estos archivos.
Sin embargo, a veces es realmente importante, como en el formato PDB, que cada elemento de datos
ocupe un grupo de columnas específico de una línea de archivo. En este caso es recomendable
el uso de la función
printf .
printf
es muy similar (no sé si idéntica) a la misma función del lenguaje C y su uso es relativamente
sencillo. Su síntaxis es:
printf("cadena formateada %s %1.2f",$escalar1,$escalar2,...);
Dentro de la cadena de caracteres formateada se puede hacer referencia a escalares que queramos
formatear. Cada una de estas referencias comienza con el símbolo % y deberá haber tantas como
escalares separados por comas en la parte izquierda de la llamada a
printf.
Además, cada referencia debe especificar el formato de salida deseado:
%c # el carácter correspondiente a un código ASCII
%s # una cadena de caracteres (string)
%d # un número entero
%f # un número real con coma fija
%g # un real con decimales
%e # un real en notación científica exponencial, como 1.204e-23
%x # un hexadecimal
Modificadores:
- # justifica a la izquierda
+ # obliga a imprimir el signo
. # fuerza un mínimo número de columnas para un campo
0 # usa ceros para justificar a la derecha
Ejemplos:
printf("%1.2f",$real); # un real de que ocupe al menos 1 columna
# y tenga un máximo de 2 decimales
printf("%.10s",$cadena); # cadena de no más de 10 columnas
Finalmente, Perl incluye la posibilidad de definir formatos de otra manera, usando las funciones
format y write,
que entenderemos mejor con un ejemplo:
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
my ($nombre,$especie,$peso);
## defino un formato y le pongo nombre
# < justifica a la izquierda, > a la derecha , #. definen los campos de números
# | se usaría para centrar contenido
format MIFORMATO =
Proteína=@<<<<<<<< Organismo=@>>>>>>>>>>>> Peso molecular=@###.#KDa
$nombre, $especie, $peso
.
## abro el archivo de salida con un descriptor cuyo nombre coincida
## con el formato de salida
open(MIFORMATO,">salida.txt") || die "no puedo crear salida.txt\n";
## aquí asigno valores a las variables del formato, normalmente será un bucle
$nombre = "1lfu";
$especie = "M.musculus";
$peso = "124.456";
write (MIFORMATO);
## imprime en el archivo de salida
## 'Proteína=1lfu Organismo= M.musculus Peso molecular= 124.5KDa'
## cierro el archivo de salida
close(MIFORMATO);
Yo reconozco que nunca uso la combinación
format/write,
así que os recomiendo, si os interesa usarlos, leer
este
documento.
Para pedir al sistema operativo que ejecute otros programas desde un programa en Perl podemos hacerlo de varias maneras. Al hacerlo estamos en realidad llamando al sistema operativo y esperando su respuesta:
## insertando el comando y sus parámetros entre dos `acentos graves`
my $resultado = `ls -l *.fas`; # muestra todos los archivos FASTA del directorio actual
$resultado = `df -k`; # muestra el espacio libre en los sistemas de ficheros de tu sistema
$resultado = `cat archivo.pdb | grep ATOM`; # imprime las líneas 'ATOM' de un PDB
## usando la función system()
$resultado = system("cat archivo.pdb | grep ATOM");
$resultado = system("pwd"); # devuelve el directorio actual
$resultado system("zcat archivo.fas.gz"); # descomprime un archivo y muestra el contenido
Perl incorpora todo un conjunto de funciones que nos permiten interaccionar con el sistema operativo, que podéis ver íntegro en la guía de referencia de Perl. Yo destacaría:
chdir $otroDir; # cambia de directorio durante la ejecución de un programa kill $IDproceso; # termina otros procesos que hayas creado anteriormente sleep $segundos; # para temporalmente la ejecución del programa que la invoca
Hay otra manera muy útil de ejecutar otros programas desde Perl y directamente procesar su salida, usando las tuberías '|' (pipes) de UNIX*. De esta manera, la salida generada por el proceso ejecutado es vista desde Perl como un archivo. Veamos ejemplos:
open(BLAST,"blastpgp -i secuencia.fas |") || die "no puedo ejecutar blastpgp -i secuencia.fas\n";
while(<BLAST>)
{
# aquí lees la salida de BLAST directamente, sin guardarla en un archivo
}
close(BLAST);
open(ZCAT,"zcat gran_archivo.gz |") || die "no puedo descomprimir gran_archivo.gz\n";
while(<ZCAT>)
{
# aquí voy leyendo el contenido de gran_archivo, pero el archivo sigue comprimido en el sistema
}
close(ZCAT);
Os propongo ahora un ejercicio basado en dos programas muy importantes de la bioinformática, BLAST y CLUSTALW(Thompson et al., 1994), que deberéis tener instalados en vuestro sistema. Las tareas que os pido son:
formatdb de BLAST.
Para terminar este capítulo os muestro un programa que usa el programa externo
curl (hay otras posibilidades, como
wget)
para hacer múltiples peticiones a un servidor web de forma automática. A veces hay que recurrir
a este tipo de soluciones, por ejemplo si la versión web es la única accesible para un programa. Pero recordad que
antes de hacer este tipo de cosas es una buena idea escribir a la persona responsable de ese servidor web y
pedirle permiso. Es posible que os recomienden qué hora del día es mejor para estos trabajos grandes, para no saturar
sus máquinas.
#!/usr/bin/perl -w
# written by Bruno Contreras, July 2005
# program get_thermo_gibbs_energies.pl, that submits a series of dna sequences in FASTA format
# to the webserver http://wings.buffalo.edu/gsa/dna/dk/WEBTHERMODYN and retrieves their Gibbs
# energy profile
# This shows how scripts can be written to serialize web form submissions
# For different webservers you need to read the source code of the the submission
# page (rigth mouse button on any browser) to check out the variables involved
# These are the variables for this example, that I got from my browser
#<form enctype="multipart/form-data" method="POST" action="webthermodyn.cgi">
#<input type=text name=temperature value=37>
#<input type=text name=saltconc value=10>
#<input type=text name=molecule value="">
#<select name="shape" value="Linear">
#<INPUT TYPE="radio" NAME="inputmode" VALUE="input">
#<textarea name="sequence" value="">
#<INPUT TYPE="radio" NAME="allpart" VALUE="ALL">
#<input type=text name=step value=50>
#<input type=text value=100 name=windowsize>
#<input type=text value=1 name=marknum>
#<input type=submit value="SUBMIT QUERY">
use strict;
######################################################################
my $progname = "get_thermo_freeenergies.pl";
my $thermodynCGI = "http://wings.buffalo.edu/gsa/dna/dk/WEBTHERMODYN/webthermodyn.cgi";
my $waittime = 5; # wait 5 seconds between submissions
my $curl_options = "-F temperature=37 -F sign=positive -F saltconc=10 -F shape=Linear -F inputmode=input -F allpart=ALL -F step=25 -F windowsize=25 -F marknum=5 -F timelimit=90 ";
## parse DNA sequences file in fasta ################################
my (@DNA,@NAME);
my ($dnaseq,$dnaseqname,$n_of_seqs) = ("","",0);
open(FASTADNA,$ARGV[0]) || die "#$progname : need a valid DNA fasta file to run properly, exit... \n";
while(<FASTADNA>)
{
next if(/^$/);
if(/\>/)
{
if($dnaseq)
{
push(@DNA,$dnaseq);
push(@NAME,$dnaseqname);
$dnaseq = "";
$dnaseqname = "";
}
$dnaseqname = substr((split)[0],1);
$n_of_seqs++;
}
else { $dnaseq .= (split)[0];}
}
close(FASTADNA);
print "#$progname : sequences read $n_of_seqs (from $ARGV[0])\n";
## get free energy profiles for each sequence ########################
my ($seq,$s);
for($seq=0;$seq<scalar(@DNA);$seq++)
{
my $request = $curl_options . "-F sequence=$DNA[$seq] -F molecule=$NAME[$seq]";
open(CURL,"curl $request $thermodynCGI |") || die "#$progname : cannot run curl $request $thermodynCGI ...exit\n";
while(<CURL>)
{
if(/DATADIR/)
{
my $results = (split(/\"\>ASCII/,$_))[0];
$results = (split(/HREF=\"/,$results))[2];
open(RES,"curl $results |" ) || die "#$progname : cannot read $results ...exit\n";
while(<RES>)
{
print;
}
close(RES);
}
}
close(CURL);
sleep $waittime; # wait this time before the next submission is delivered
}
Se podría decir que con lo visto hasta ahora podrías ya lanzaros a atacar problemas bioinformáticos con Perl. Sin embargo, todavía no sabéis cómo escribir programas de forma estructurada, con bloques de código separados en forma de subrutinas, ni cómo aprovechar módulos escritos por otra gente y disponibles en Internet. También es útil crear tus propias estructuras de datos, más allá de los arreglos y las tablas, para representar de forma adecuada datos complejos. Por todas estas razones vamos a seguir profundizando un poco más en Perl, para que saquéis el máximo partido al lenguaje y a los millones de líneas de código Perl que hay distribuidas por Internet.
Una subrutina puede definirse como un programa dentro de un programa, un subprograma.
Es un trozo de código que tiene vida propia, con sus propio algoritmo y sus propias
variables, que puede estar declarado dentro de otro programa Perl, dentro de un módulo o
dentro del intérprete Perl, como
print,
por ejemplo.
Para ejecutar una subrutina hay que invocarla por su nombre, y si fuera necesario sus
parámetros. Por ejemplo ya hemos invocado
printf("cadena %s\n",$cadena);
antes, donde
printf
es el nombre de la subrutina y lo que va entre paréntesis sus parámetros.
Si invocamos una función y el intérprete no encuentra su declaración, se interrumpirá la
ejecución en ese momento con un mensaje de error como este:
Undefined subroutine &main::subrutina_sin_de called at programa.pl line 25.
Declarar una subrutina en Perl es muy sencillo, es un bloque de código precedido de la palabra reservada sub y un nombre, que, como las variables, debe empezar por un carácter alfanumérico o un '_'.
Puedes declarar subrutinas en cualquier posición de un programa en Perl, pero por cuestión de estilo es costumbre declararlas todas juntas por ejemplo al final del programa principal.
![\begin{figure}
\begin{center}
\includegraphics[width=0.8\textwidth]{subrutinas}
\end{center}
\end{figure}](img3.png)
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
mi_subrutina(); # invocación de la subrutina: el intérprete busca
# ahora la declaración de mi_subrutina y la encuentra
# más adelante en el mismo archivo; en ese momento
# interrumpe la ejecución del programa principal y
# ejecuta el código de la subrutina. El programa
# principal se reanuda al terminar mi_subrutina
######## declaraciones de subrutinas ######################################
sub mi_subrutina
{
# código propio de la subrutina
print "mi primera subrutina\n";
}
Si la ejecución de la subrutina requiere el paso de datos externos a la subrutina, normalmente esos datos se pasan como parámetros en el momento de la invocación:
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
my ($arg1,$arg2) = (1,2);
my @arg3 = (3,4,5);
mi_subrutina($arg1,$arg2,@arg3); # paso de tres parámetros
######## declaraciones de subrutinas ######################################
sub mi_subrutina
{
my($sarg1,$sarg2,@sarg3) = @_; # así recibe la subrutina los parámetros pasados,
# copiándolos del arreglo global @_ a variables locales
print "subrutina con parámetros\n";
}
El paso de parámetros a una subrutina
depende de la variable especial
@_,
que ya mencionamos en 2.3.4, que es un ejemplo
de variable global , a la que se puede acceder desde cualquier
lugar de un programa.
A diferencia de
@_,
las variables que se declaran dentro de una subrutina, como
my($sarg1,$sarg2,@sarg3) ,
son por defecto variables locales,
que dejan de existir cuando se termina la ejecución de la subrutina y que no pueden ser accedidas
si no es desde la misma subrutina.
Aunque en teoría una subrutina puede acceder a variables globales, declaradas fuera de ella, es conveniente limitar el accesos a aquellas variables que se pasen como parámetros, para facilitar la comprensión del código que escribáis y la corrección de posibles errores.
A parte de
@_ ,
las subrutinas pueden acceder a otras variables globales. Para probar
esto os pido que ejecutéis este programita, observéis la salida y luego lo modifiquéis para
demostrar que la subrutina puede hacer exactamente los mismos cálculos recibiendo como
parámetros los datos necesarios. ¿No os parece que se hace más fácil comprender el
funcionamiento de
printBLAST
en la versión modificada, ya que no hay que buscar en el
código de donde vienen las variables utilizadas?
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
my ($evalue,$identities,$queryseq,$subjectseq);
# 100 líneas de código
$evalue = 4e-34;
$identities = 87;
$queryseq = "MARRKRRNFNKQATEILNEYFYSHLSNPYPSXXXXXXXXXXSGITVSQVSNWFGNKRIRYKKNIGKFQEEANIYAAKTAVTA";
$subjectseq = "MARRKRRNFNKQATEILNEYFYSHLSNPYPSEEAKEELAKKSGITVSQVSNWFGNKRIRYKKNIGKFQEEANIYAAKTAVTA";
# programa principal
my @blastdata = ($evalue,$identities,$queryseq,$subjectseq);
# 100 líneas de código
printBLAST();
# 500 líneas de código
######## declaraciones de subrutinas ######################################
sub printBLAST
{
print "evalue $blastdata[0]\n";
print "alignment:\n$blastdata[2]\n$blastdata[3]\n";
}
Ya que sabemos cómo declarar e invocar subrutinas, ahora veremos cómo podemos escribir subrutinas que devuelvan valores. Una vez más, creo que un ejemplo deja todo más claro:
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
## variables del programa principal
my $queryseq = "MARRKRRNFNKQATEILNEYFYSHLSNPYPSXXXXXXXXXXSGITVSQVSNWFGNKRIRYKKNIGKFQEEANIYAAKTAVTA";
my $subjectseq = "MARRKRRNFNKQATEILNEYFYSHLSNPYPSEEAKEELAKKSGITVSQVSNWFGNKRIRYKKNIGKFQEEANIYAAKTAVTA";
## llamada a printf que a su vez llama a dos subrutinas que devuelven valores
printf("Alineamiento %2.2f :\n%s\n%s\n%s\n",
calculaIdentidad($queryseq,$subjectseq), # llamada a una subrutina
decoraAlineamiento($queryseq,$subjectseq)); # llamada a otra subrutina
## no siempre imprimiremos los valores devueltos por una subrutina, a veces simplemente los
## guardaremos:
my ($dev1,$dev2,$dev3) = decoraAlineamiento($queryseq,$subjectseq);
######## declaraciones de subrutinas ######################################
sub calculaIdentidad
{
# parámetros
my ($sec1,$sec2) = @_;
# variables locales
my ($id,$pos,$long) = (0,0,length($sec1));
for($pos=0;$pos<$long;$pos++)
{
if(substr($sec1,$pos,1) eq substr($sec2,$pos,1))
{
$id++;
}
}
return sprintf("%2.2f",$id/$long); # devuelve un número real
}
###########################################################################
sub decoraAlineamiento
{
# parámetros
my ($sec1,$sec2) = @_;
# variables locales
my ($decoracion,$pos,$long) = ("",0,length($sec1));
for($pos=0;$pos<$long;$pos++)
{
if(substr($sec1,$pos,1) eq substr($sec2,$pos,1))
{
$decoracion .= '|';
}
else
{
$decoracion .= ' ';
}
}
return ($sec1,$decoracion,$sec2); # devuelve tres cadenas de caracteres
}
Es común llamar funciones a las subrutinas que devuelven valores
y procedimientos a las que no. Si una subrutina contiene
la función
return,
entonces será una función.
A diferencia de otros lenguajes, en Perl una función puede devolver más de un valor, como habéis
comprobado con
decoraAlineamiento.
A diferencia de otros lenguajes, en Perl los procedimientos pueden invocarse de manera diferente
a las funciones, usando
&procedimiento;
en vez de
procedimiento();.
En los ejemplos vistos hasta ahora, cuando invocamos una subrutina y le pasamos parámetros
en realidad le estamos pasando una lista de datos, que se llama
@_ .
Cuando pasamos arreglos o tablas como parámetros en realidad los estamos vertiendo
elemento por elemento en
@_ ,
lo cual parece un poco estúpido a veces, con el
trabajo que nos dio previamente llenarlos. De hecho, si queremos pasar
por ejemplo dos arreglos a la vez, será imposible que la subrutina que los reciba
sepa donde termina uno y empieza el otro, porque lo que está recibiendo es
@_ .
Por esta razón, cuando pasamos varios parámetros y uno de ellos es un arreglo, debemos
poner al arreglo en última posición, para poder identificarlos después en la subrutina.
Lo mismo se aplica a los valores que devuelve una subrutina con
return.
## forma correcta de pasar más de un parámetro y recibir más de un resultado
my ($result1,@result2) = subrutina($arg1,$arg2,@arg3);
sub subrutina
{
my($sarg1,$sarg2,@sarg3) = @_;
my ($sresult1,@sresult2);
return ($result1,@result2);
}
## en cambio, esto no funcionaría
my (@result2,$result1) = subrutina($arg1,@arg3,$arg3); # $result1 queda sin asignar
sub subrutina
{
my($sarg1,@sarg3,$sarg3) = @_; # $sarg3 queda sin asignar
my ($sresult1,@sresult2);
return (@result2,$result1);
}
La forma más económica de pasar parámetros a subrutinas es pasarlos como referencias (ver sección
2.3). De esta manera no pasamos los datos en sí, si no la dirección que ocupan en
memoria. Así no es necesario copiar todos los datos a
@_
y la subrutina que reciba las
referencias puede acceder a los datos en su organización original, ya sean escalares, arreglos
o tablas asociativas.
Otra aplicación de las referencias parámetro es que permiten modificar de forma permanente, no local, los datos pasados como argumento. Finalmente, el paso por referencias permite pasar como parámetro otras subrutinas, lo cual puede ser útil en ocasiones.
Veamos un ejemplo:
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
## variables del programa principal
my ($escalar , @arreglo1 , @arreglo2);
## asignación de valores
$escalar = 0;
@arreglo1 = (1,2,3);
@arreglo2 = (4,5,6);
## llamada a subrutina y conversión de la referencia recibida a un arreglo
my @resultado = @{ lee_arreglos( $escalar , \@arreglo1 , \@arreglo2 ) };
print "resultado: $resultado[0] , $resultado[1] , $resultado[2]\n";
############### subrutinas ##################
sub lee_arreglos
{
# parámetros pasados
my ($sarg1,$rarreglo1,$rarreglo2) = @_;
# variables locales
my @res = ( $sarg1, $rarreglo1->[0], $rarreglo2->[0] );
return \@res; # devuelve una referencia a un arreglo como resultado
}
Como en otros lenguajes, en Perl se pueden programar funciones recursivas, que se llaman a si mismas hasta que el problema que intentan resolver tiene un solución obvia. La recursión es un mecanismo de descomposición de problemas que puede hacer más sencillo llegar a una solución, pero a cambio, normalmente, su coste computacional es relativamente grande, ya que implica muchas llamadas a funciones seguidas.
De ejemplo os muestro la función factorial en forma recursiva:
sub factorial
{
# toma parámetro
my $fact = $_[0];
if($fact < 0) # error, la función factorial no está definida para negativos
{
return 0;
}
elsif($fact == 0) # caso trivial
{
return 1;
}
else # descompón el problema un paso más
{
return $fact * factorial($fact -1);
}
}
Para aplicar lo que hemos aprendido de subrutinas os pediría que modifiquéis la solución del problema 2.10.1 ahora separando todas las tareas indivisibles en subrutinas. Podéis además comparar la velocidad de ejecución de ambos programas, el viejo y el modificado, usando el programa time en el terminal UNIX*.
Otra manera de valorar el rendimiento de un programa que use subrutinas es ejecutarlo de esta
manera: perl -d:DProf programa.pl [argumentos que necesite] . Tras su ejecución se habrá
generado un archivo llamado tmon.out que pasaremos como argumento a otro programa:
dprofpp tmon.out , generando una salida como esta, donde se ve qué coste empírico
tiene cada subrutina durante la ejecución del programa:
Exporter::export has -1 unstacked calls in outer Exporter::Heavy::heavy_export has 1 unstacked calls in outer Total Elapsed Time = 34.29820 Seconds User+System Time = 16.72820 Seconds Exclusive Times %Time ExclSec CumulS #Calls sec/call Csec/c Name 64.9 10.85 10.858 124763 0.0000 0.0000 dnalib::get_sq_dist 23.0 3.852 7.964 3 1.2841 2.6548 dnalib::get_template_contacts 14.4 2.410 7.221 588 0.0041 0.0123 dnalib::check_base_contact 6.31 1.055 2.990 3 0.3517 0.9968 dnalib::calc_DNA_complementarity 3.03 0.507 11.056 3 0.1690 3.6853 dnalib::predict_sites_complex 0.72 0.121 0.121 28932 0.0000 0.0000 dnalib::tobase4 0.50 0.084 0.084 30720 0.0000 0.0000 dnalib::num2compbase 0.42 0.070 0.070 1 0.0700 0.0700 dnalib::extract_BLAST_TF_alignments 0.41 0.069 0.069 581 0.0001 0.0001 dnalib::calc_Pvalue 0.41 0.069 0.069 598 0.0001 0.0001 dnalib::calc_family_preferences 0.24 0.040 0.050 5 0.0080 0.0100 main::BEGIN 0.23 0.039 8.113 1 0.0391 8.1126 dnalib::build_CM_complexes 0.20 0.034 0.034 30720 0.0000 0.0000 dnalib::num2base 0.18 0.030 0.030 3 0.0100 0.0100 dnalib::get_protein_coordinates 0.17 0.029 0.029 581 0.0001 0.0001 dnalib::calc_def_cost
En esta sección veremos cómo crear y manipular estructuras de datos más complejas que las que hemos visto hasta ahora, basadas sobre todo en el uso de arreglos, tablas asociativas y referencias. Como guía alternativa de referencia para este contenido os recomiendo Data structures cookbook.
En Perl es muy fácil anidar arreglos dentro de arreglos, por lo que podemos crear matrices N-dimensionales, como ya mencionamos en la sección 2.3. Hay varias maneras de crearlas y consultarlas, que os mostraré con ejemplos:
my @matriz2_3 = ( [1,2,"tres"], ## declaración con asignación explícita
[4,5,"seis"]);
my @matriz;
for($i=0;$i<5;$i++) ## asignación elemento a elemento desde bucles
{
for($j=0;$j<10;$j++)
{
$matriz[$i][$j] = $i * $j;
}
}
for $i ( 1 .. 10 ) ## asignación de filas enteras
{
# ambas expresiones son válidas
push( @matriz , [ @dimension ];
$matriz[$i] = [ @dimension ];
# $matriz[$i] contiene una referencia a un arreglo
}
while(<ARCHIVO>) ## leyendo un archivo línea por línea
{
push( @matriz , split(/\s+/,$_);
}
Las tablas de tablas son también sencillas de manejar, de forma muy parecida a los arreglos de arreglos.
my %blosum62 = ( ## declaración con asignación explícita
'A' =>
{
'A' => -4,
'R' => -1,
...
'X' => -1
}
...
'X' =>
{
'A' => --,
'R' => -1,
...
'X' => -1
}
);
my %tabla;
for($i=0;$i<5;$i++) ## asignación elemento a elemento desde bucles
{
for($j=0;$j<10;$j++)
{
$tabla{$i}{$j} = $i * $j;
}
}
for $i ( 1 .. 10 ) ## asignación de subtablas enteras
{
$tabla[$i] = { %subtabla };
# tabla[$i] contiene una referencia a una tabla
}
La posibilidad de hacer tablas de arreglos y arreglos de tablas es muy útil en ocasiones,
puesto que nos permite combinar las funciones de Perl sobre arreglos, como
push,
con la facilidad de de uso de las tablas asociativas. La única complicación es que
al crear una estructura de datos como esta debemos tener en cuenta sus anidamientos
para manipularla, utilizando funciones como
keys o
values
para las subtablas y
sort o
foreach para los subarreglos:
my @PDBcoordenadas = (
# una tabla para cada resduo de una proteína PDB
{
'N' => 'ATOM 888 N MET P 0 13.559 11.810 -20.485',
'CA' => 'ATOM 889 CA MET P 0 12.699 10.612 -20.293',
'C' => 'ATOM 890 C MET P 0 12.893 10.012 -18.905',
'O' => 'ATOM 891 O MET P 0 11.938 9.870 -18.140'
},
...
);
# podríamos recorrer @PDBcoordenadas recorrería así:
foreach $residuo (@PDBcoordenadas)
{
# cada valor de $residuo es una referencia a una tabla asociativa
print $residue->{'CA'};
}
my %codones = (
# un arreglo para los codones de cada aminoácido
'F' => ['TTT','TTC'],
'L' => ['TTA','TTG','CTT','CTC','CTA','CTG'],
'S' => ['TCT','TCC','TCA','TCG','AGT','AGC'],
'Y' => ['TAT','TAC'],
...
# $codones{'L'} sería una referencia a un arreglo
);
# que recorreríamos:
foreach my $amino (keys(%codones)) # recorre llave a llave
{
print "$amino ";
foreach my $codon ( @{ $codones{$amino} } ) # recorre los codones asociados a una llave
{
print "$codon,";
}
print "\n";
}
Los registros son estructuras de datos que engloban componentes de diferente tipo y
equivalen a los que en C serían los objetos
struct.
Por supuesto puedes crear arreglos y tablas de registros:
my $proteina = {
NOMBRE => $nombre,
ANOTACION => $anot,
PM => $PM,
PI => $PI,
NOMBRES_HOMOLOGOS => [ @homologos ],
ALINEAMIENTOS_HOMOLOGOS => { %alineamientos }
};
# así consultamos algún valor
print $proteina->{'NOMBRE'};
# así asignamos
$proteina->{'NOMBRE'} = 'CRP_ECOLI';
# así la guardamos en un arreglo
push(@proteinas,$proteina);
Las colas (queue, FIFO, First In First Out) y las pilas (stack, LIFO, Last In First Out) son dos ejemplos clásicos de estructuras de datos en los libros de texto. Por supuesto se pueden implementar en Perl por medio de registros que se apuntan unos a otros, pero lo más sencillo es usar para arreglos y algunas de las funciones vistas en la sección 2.7.1.
Si usamos la combinación
shift y
unshift, o
push y
pop,
tenemos una pila.
Si usamos
shift y
push
tenemos una cola.
Podemos utilizar los registros que acabamos de ver para representar árboles y, de modo más general, grafos. Cada nodo será un registro que apuntará a sus vecinos según la topología del grafo correspondiente.
Os muestro un ejemplo muy sencillo para un nodo de un árbol binario, donde cada nodo tiene a lo sumo dos hijos, el derecho y el izquierdo:
![\begin{figure}
\begin{center}
\includegraphics[width=0.8\textwidth]{colaypila}
\end{center}
\end{figure}](img4.png)
![\begin{figure}
\begin{center}
\includegraphics[width=0.8\textwidth]{arbolBin}
\end{center}
\end{figure}](img5.png)
my $nodoBin = {
VALOR => $valor,
IZQ => $rizq, # referencia al nodo hijo izquierdo
DER => $rder # referencia al nodo hijo derecho
};
Si un nodo puede conectarse a un número indeterminado de nodos, entonces es mejor incluir las referencias a esos nodos en un arreglo.
Os muestro aquí un sencillo programa para manipular
árboles binarios de búsqueda,
tomado de aquí.
¿Podéis entender el código?¿Qué tiene de especial la subrutina
search?
#!/usr/bin/perl -w
# bintree - binary tree demo program
use strict;
my($root, $n);
# first generate 20 random inserts
while ($n++ < 20) { insert($root, int(rand(1000)))}
# now dump out the tree all three ways
print "Pre order: "; pre_order($root); print "\n";
print "In order: "; in_order($root); print "\n";
print "Post order: "; post_order($root); print "\n";
# prompt until EOF
for (print "Search? "; <>; print "Search? ") {
chomp;
my $found = search($root, $_);
if ($found) { print "Found $_ at $found, $found->{VALUE}\n" }
else { print "No $_ in tree\n" }
}
exit;
#########################################
# insert given value into proper point of
# provided tree. If no tree provided,
# use implicit pass by reference aspect of @_
# to fill one in for our caller.
sub insert {
my($tree, $value) = @_;
unless ($tree) {
$tree = {}; # allocate new node
$tree->{VALUE} = $value;
$tree->{LEFT} = undef;
$tree->{RIGHT} = undef;
$_[0] = $tree; # $_[0] is reference param!
return;
}
if ($tree->{VALUE} > $value) { insert($tree->{LEFT}, $value) }
elsif ($tree->{VALUE} < $value) { insert($tree->{RIGHT}, $value) }
else { warn "dup insert of $value\n" }
# XXX: no dups
}
# recurse on left child,
# then show current value,
# then recurse on right child.
sub in_order {
my($tree) = @_;
return unless $tree;
in_order($tree->{LEFT});
print $tree->{VALUE}, " ";
in_order($tree->{RIGHT});
}
# show current value,
# then recurse on left child,
# then recurse on right child.
sub pre_order {
my($tree) = @_;
return unless $tree;
print $tree->{VALUE}, " ";
pre_order($tree->{LEFT});
pre_order($tree->{RIGHT});
}
# recurse on left child,
# then recurse on right child,
# then show current value.
sub post_order {
my($tree) = @_;
return unless $tree;
post_order($tree->{LEFT});
post_order($tree->{RIGHT});
print $tree->{VALUE}, " ";
}
# find out whether provided value is in the tree.
# if so, return the node at which the value was found.
# cut down search time by only looking in the correct
# branch, based on current value.
sub search {
my($tree, $value) = @_;
return unless $tree;
if ($tree->{VALUE} == $value) {
return $tree;
}
search($tree->{ ($value < $tree->{VALUE}) ? "LEFT" : "RIGHT"}, $value)
}
Como ejercicio de esta sección os sugiero que escribáis un programa que tome un grupo
de proteínas o de ácidos nucleicos de la misma familia, en formato FASTA, y calcule
mediante el programa
bl2seq,
del paquete BLAST, una matriz de similaridad de
todas las secuencias contra todas. Finalmente, basándoos en la matriz, el programa
deberá crear un árbol que represente las distancias de todas las secuencias a una
elegida por el usuario.
Las estructuras de datos que hemos visto en la sección anterior permiten manejar cómodamente relaciones abstractas entre datos, como por ejemplo un grafo. Los objetos llevan un poco más ella esa abstracción, pues son agregados de variables (atributos) con identidad propia que además incluyen las subrutinas (métodos) que permiten manipularlos.
Un objeto es una instancia de la clase a la que pertenece, como por ejemplo un Mercedes con placas 1234HJ pertenecería a la clase de los automóviles de lujo, o la proteína 1lfu pertenece a la clase PDB.
Los atributos de un objeto normalmente sólo pueden ser manipulados por medio de los métodos de su clase.
Una clase clase puede heredar componentes de otras clases, por lo que son adecuadas para representar conceptos y objetos de forma jerárquica. Por ejemplo, la clase PDB será seguramente una clase derivada de la clase proteína y ésta a su vez de la clase macromolécula, de donde heredarán el atributo 'secuencia'.
Os recomiendo que miréis en Google si no os resultan familiares estos conceptos, en el libro Advanced Perl Programming, pero recordad que la programación orientada a objetos es un estilo de programación, que tiene sus ventajas e incovenientes, que conviene conocer al menos para poder comprender programas escritos por otras personas. Como ya sabéis, Perl os dejará programar usando programación estructura, orientada a objetos y combinaciones de ambas sin problema.
Una manera de implementar los atributos de un objeto en Perl sería muy similar a los registros de la sección 3.2.4 y se podrían generar así:
![\begin{figure}
\begin{center}
\includegraphics[width=0.8\textwidth]{clase}
\end{center}
\end{figure}](img6.png)
my $pdb1 = constructor_PDB('1lfu',3,'1lfu.pdb');
sub constructor_PDB # un método que genera nuevos objetos se llama constructor
{
# valores de los atributos como parámetros
my($nombre,$n_de_cadenas,$archivoPDB) = @_;
my $rPDB = {
'nombre' => $nombre,
'ncadenas' => $n_de_cadenas,
'archivo' => $archivoPDB,
'coordenadas' => ()
};
return $rPDB; # devuelve una referencia al objeto
}
Los métodos de un objeto permiten manipular sus atributos desde el exterior, sin necesariamente conocer las tripas del objeto:
leer_archivoPDB($pdb1,0);
sub leer_archivoPDB
{
# parámetros
my($rPDB,$comprimido) = @_;
# variables locales
my $res_leidos = 0;
if($comprimido == 1)
{
open(PDB,"zcat $rPDB->{'archivo'} |") || warn "no puedo leer zcat $rPDB->{'archivo'}\n";
# rellena $rPDB->{'coordenadas'}
# $res_leidos++;
close(PDB);
}
else
{
open(PDB,"$rPDB->{'archivo'}") || warn "no puedo leer $rPDB->{'archivo'}\n";
# rellena $rPDB->{'coordenadas'}
# $res_leidos++;
close(PDB);
}
return $res_leidos;
}
Acabamos de implementar un objeto de la clase PDB usando registros y subrutinas, ahora lo
haremos usando la sintaxis propia de la programación orientada a objetos. Para ello debemos
crear un paquete
.pm
que contendrá el código de la clase. Como ejemplo os muestro la
clase PDB, contenida en el archivo PDB.pm, donde veréis que usamos la función
bless
que permite distinguir al intérprete entre referencias y objetos:
package PDB; # se incluye en otro programa con 'use PDB.pm;'
use strict;
sub PDB::constructor # en inglés este método suele llamarse 'new'
{
# autoreferencia a la clase y atributos como parámetros
my($clase,$nombre,$n_de_cadenas,$archivoPDB) = @_;
# implementación de atributos como una tabla asociativa
my $objeto = {
'nombre' => $nombre,
'ncadenas' => $n_de_cadenas,
'archivo' => $archivoPDB,
'coordenadas' => [],
};
# podría haber sido como un arreglo, más eficiente, pero menos obvio el contenido
#$objeto->[0] = $nombre;
#$objeto->[1] = $n_de_cadenas;
#$objeto->[2] = $archivoPDB;
#$objeto->[3] = [];
bless($objeto,$clase); # la función bless da carácter de objeto a la referencia $objeto
return $objeto;
}
sub PDB::pon_nombre_archivo
{
my($objeto,$archivoPDB) = @_;
$objeto->{'archivo'} = $archivoPDB;
}
sub PDB::leer_archivoPDB
{
# autoreferencia al objeto y parámetros
my($objeto, $comprimido) = @_;
# variables locales
my @coordenadas = ();
if(!$objeto->{'archivo'})
{
print "'archivo' no está definida\n";
return 0;
}
if($comprimido == 1)
{
open(PDB,"zcat $objeto->{'archivo'} |") || warn "no puedo leer zcat $objeto->{'archivo'}\n";
}
else
{
open(PDB,"$objeto->{'archivo'}") || warn "no puedo leer $objeto->{'archivo'}\n";
}
while(<PDB>)
{
if(/^ATOM/||/^TER/)
{
push(@coordenadas,$_);
}
}
close(PDB);
$objeto->{'coordenadas'} = [ @coordenadas ];
return scalar(@coordenadas);
}
sub PDB::imprime_coordenadas
{
my($objeto) = @_;
print @{ $objeto->{'coordenadas'} };
}
1;
A continuación os muestro la clase derivada PDBdna, contenida en un archivo llamado
PDBdna.pm, que hereda de la clase PDB y sólo modifica el método
leer_archivoPDB:
package PDBdna; # se incluye en otro programa con 'use PDBdna.pm;'
use strict;
# clase derivada de PDB, especializada en representar moléculas de ADN
use PDB;
use vars '@ISA';
@ISA = 'PDB';
# redefino aquí el método leer_archivoPDB, para tener en cuanta sólo los
# átomos de ADN
sub PDBdna::leer_archivoPDB
{
# autoreferencia al objeto y parámetros
my($objeto, $comprimido) = @_;
# variables locales
my @coordenadas = ();
if(!$objeto->{'archivo'})
{
print "'archivo' no está definida\n";
return 0;
}
if($comprimido == 1)
{
open(PDB,"zcat $objeto->{'archivo'} |") || warn "no puedo leer zcat $objeto->{'archivo'}\n";
}
else
{
open(PDB,"$objeto->{'archivo'}") || warn "no puedo leer $objeto->{'archivo'}\n";
}
while(<PDB>)
{
# lee sólo átomos de ADN
if((/^ATOM/ && substr($_,18,1) eq " " && substr($_,19,1) =~ /A|G|C|T/) || /^TER/)
{
push(@coordenadas,$_);
}
}
close(PDB);
$objeto->{'coordenadas'} = [ @coordenadas ];
return scalar(@coordenadas);
}
1;
Finalmente, os muestro un ejemplo que deberéis probar donde se utilizan estas clases:
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Julio de 2005
use strict;
use PDB; # incluye PDB.pm
use PDBdna; # incluye PDBdna.pm, en el mismo directorio
# os sugiero que modifiquéis el código para probar también la clase PDBdna
my $pdb1 = PDB->constructor('1lfu',3,'bioinfoPerl/problemas/1lfu.pdb');
$pdb1->leer_archivoPDB(0);
$pdb1->imprime_coordenadas();
En la sección 3.3.3 creamos un paquete de Perl en un archivo
.pm .
En concreto estábamos creando una clase, pero los paquetes de Perl sirven en general para agrupar código,
en forma de variables globales, estructuras de datos, subrutinas y programas. De esta forma,
tal cómo decíamos en la sección 3.1.5, favorecemos la reutilización y distribución
de nuestro código y facilitamos la corrección de errores.
En UNIX*, si tecleáis
$ man perlmodlib
podéis ver el listado y documentación sobre
los paquetes estándar de Perl, entre los que está, ya lo recordaréis,
strict.
Es muy recomendable que echéis un vistazo aquí antes de poneros a escribir nuevos módulos, no vaya
a ser que ya existan. La colección es impresionante, amplia las capacidad que ya conocéis de
Perl de forma considerable.
En el archivo en red CPAN (Comprehensive Perl Archive Network) podréis encontrar, con fecha de Julio de 2005, más de 8000 paquetes que casi 5000 autores han escrito y documentado para un sinfín de aplicaciones distintas. Este repositorio es sin duda uno de los mayores atractivos de Perl.
Una vez que hayáis encontrado un paquete que os interese, lo descargáis y tendréis un archivo
modulo.tar.gz
que deberéis descomprimir con
tar xvfz modulo.tar.gz.
En el subdirectorio creado tenéis un archivo
README
donde se explican los pasos para instalar
y probar el módulo. Sólo para el último de ellos necesitaréis privilegios de administrador,
pero con el
make y
make test
ya podréis usarlo.
Si tenéis privilegios de administrador, la manera más fácil es teclear:
perl -MCPAN -e "install 'Modulo'"
Podemos crear un módulo con subrutinas o estructuras de datos que nos hayan sido especialmente
útiles y que sigamos usando de forma repetida. En teoría basta con que creéis un archivo
paquete.pm
que contenga las siguientes líneas:
package paquete; # todo lo que venga aquí es el contenido del módulo 1;
Ahora os muestro un ejemplo de paquete genérico que contiene dos módulos:
package paquete;
################### módulo mod1 ########################
$mod1::var1 = 'mod1::variable1';
sub mod1::subrutina1
{
# código de subrutina
}
# ...
################### módulo mod2 ########################
$mod2::var1 = 'mod2::variable1';
sub mod2::subrutina1
{
# código de subrutina
}
########################################################
1; # convención para terminar un paquete
Para utilizar este módulo desde otro programa podríamos:
#!/usr/bin/perl -w # Ejemplo escrito por Bruno Contreras en Julio de 2005 use strict; use paquete; print $mod1::var1; # acceso a una variable del módulo mod1 de paquete mod2::subrutina1(); # acceso a una subrutina del módulo mod2 de paquete
Al instalar el intérprete Perl viene incluido el programa
h2xs
que nos ayuda a crear paquetes siguiendo una plantilla genérica. Para crear el paquete ejemplo haríamos:
$ h2xs -Xn ejemplo
.
Que a mi me devuelve:
Defaulting to backwards compatibility with perl 5.8.5 If you intend this module to be compatible with earlier perl versions, please specify a minimum perl version with the -b option. Writing ejemplo/lib/ejemplo.pm Writing ejemplo/Makefile.PL Writing ejemplo/README Writing ejemplo/t/ejemplo.t Writing ejemplo/Changes Writing ejemplo/MANIFEST
El programa crea un subdirectorio con el nombre del paquete y coloca ahí una plantilla
ya iniciada de
ejemplo.pm
donde deberemos ir añadiendo el código que deseemos.
Además, el programa genera documentación que acompaña al paquete e instrucciones
para compilarlo y finalmente instalarlo en tu sistema.
Si queréis contribuir con algún módulo a CPAN deberán ir en el formato generado por
h2xs.
Como ya vimos, un paquete se incluye en un programa por medio de
use paquete; .
Cuando el intérprete Perl llega a esa línea, normalmente en el encabezado de un programa,
busca un archivo que se llame
paquete.pm
en una serie de directorios, que incluyen el directorio actual, almacenados en el arreglo especial
@INC.
Podemos comprobar el contenido de
@INC
simplemente imprimiéndolo con
print @INC; .
Si queremos utilizar un módulo o paquete que no se encuentra en
@INC,
deberemos añadir estas
líneas antes de llamar a
use:
use lib "/home/usuario/modulos/" use strict; use mimodulo;
Para que vayáis practicando con estas cosas os voy a pedir que escribáis el paquete
contactos.pm.
Se trata de un paquete relacionado con el cálculo de mapas de
contactos de proteínas. Un mapa de contactos para una proteína con N aminoácidos
es una matriz binaria NxN donde en cada celda (coordenadas i,j) habrá un '+' si los
aminoácidos i y j están a menos de 7Å uno del otro y un ' ' en caso contrario:
001 ++++ + 002 +++++ + 003 +++++ + ...
Vuestro paquete deberá contener al menos 4 subrutinas que:
Si queréis usar el programa
h2xs
os recomiendo esta
guía.
Una vez que tengáis el módulo listo y documentado, escribid un programa que lo invoque y utilice que sirva de ejemplo.
En esta sección veremos ejemplos de algunos paquetes y módulos que me parecen muy útiles para distintas aplicaciones en bioinformática. Algunos vienen incorporados con el intérprete, otros son de CPAN y otros los obtuve directamente de las páginas web de los autores.
Este módulo facilita la creación de una interfaz para tu programa en la línea de comando, definiendo las diferentes opciones dispibles y manejando su tratamiento. Lo mostraré mejor con un ejemplo:
#!/usr/bin/perl -w
# program: example.pl
# author: Bruno Contreras Moreira
use strict;
use Getopt::Std;
my %opts;
getopts('hsi:l:', \%opts);
if(($opts{'h'})||(scalar(keys(%opts))==0))
{
print "usage: progname.pl [options]\n";
print "-h \t this message\n";
print "-i \t input transcription file\n";
print "-s \t skip preprocessing (optional)\n";
print "-l \t oligo length (default 1, optional)\n\n";
exit;
}
if(!defined($opts{'i'})) { die "# need an input file\n"; }
if(!defined($opts{'l'}) || ($opts{'l'} < 1) )
{
$opts{'l'} = 1;
}
if(!defined($opts{'s'}))
{
$opts{'s'} = 0;
}
printf("# example.pl -i %s -s %s -l %d\n",$opts{'i'},$opts{'s'},$opts{'l'});
Esta colección de paquetes de CPAN nos permite desarrollar clientes web en Perl, os recomiendo que investiguéis la documentación en CPAN para ver todas sus posibilidades.
Aquí os mostraré el uso del paquete LWP::Simple para descargarnos una página web, lo cual puede sernos de utilidad alguna vez:
#!/usr/bin/perl -w
# cómo descargarse una página web de forma no interactiva
use strict;
use LWP::Simple;
my $contenido = get("http://www.perl.com");
El XML es un lenguaje descriptivo que sirve para crear documentos estructurados y representar con ellos datos abstractos. Se está usando cada vez más como formato para distribuir documentos e información en bioinformática, ya que se pueden manipular desde muchas plataformas y lenguajes de programación. Esto sería un ejemplo de documento XML:
<?xml version='1.0' standalone='yes' ?> <documento> <encabezado> <titulo>Módulos sobre XML en Perl</titulo> <autor>CPAN</autor> </encabezado> <texto> XML::API Perl extension for creating XML documents XML-API-0.06 - 20 Jul 2005 - Mark Lawrence XML::ASCX12 ASCX12 EDI to XML Module XML-ASCX12-0.03 - 28 Sep 2004 - Brian Kaney XML::ASX An ASX file - methods for everything from <ASX> to </ASX> XML-ASX-0.01 - 26 Jun 2002 - Allen Day XML::ApacheFOP Access Apache FOP from Perl to create PDF files using XSL-FO. XML-ApacheFOP-0.01 - 01 Jun 2005 - Ken Prows ... </texto> </documento>
En Perl, sólo en CPAN hay sobre 1000 paquetes relacionados con la manipulación de documentos XML. Podéis ver ejemplos en xmlbio.org y en este tutorial.
Módulo para buscar sitios de restricción como cliente del servidor RestrictionMapper, muy sencillo de usar.
#!/usr/bin/perl -w
# program that gets the E.coli restriction maps of a series of E.coli sequences in fasta format
# based on test.pl from http://www.restrictionmapper.org/remote.htm
use lib "/home/compu2/Perl/remoterestmap/";
use strict;
use RemoteRestMap;
my %settings = (
isoschizomers => "yes",
enzymelist => ['EciI','M.Eco67Dam','EcoHAI','EcoNI','Eco29kI','EcoprrI']
);
# instantiation of map object
my $map_obj = new RemoteRestMap("http://www.restrictionmapper.org/cgi-local/sitefind3.pl", "atgc");
$map_obj->sequence("ATCGATCGATACGCGATACTAGCGCGATCATCGCGCGCTATATACTCAGCTCGCTCTAGATATCTA");
my $map = $map_obj->get_map(\%settings);
my @sites;
while (my $cuts = $map->get_next_enz())
{
push( @sites , split(/\,/,$cuts->{CUTLIST}->[0]) );
}
Este paquete permite hacer diagramas y gráficas para representar datos numéricos de manera
sencilla. Puedes obtenerlo en CPAN y deberás instalar el paquete
GD
antes de instalar
GD::Graph,
aunque este es un ejemplo de paquete que puede dar problemas al instalar,
porque requiere a su vez otros componentes. Pero dada su utilidad creo que vale la pena
el trago.
La documentación para explotar este paquete la encontráis en GD::Graph, en CPAN.
Aquí tenéis el código para hacer un diagrama X/Y de dos variables numéricas, que se muestra en el archivo ejemplo.png:
#!/usr/bin/perl -w # Ejemplo escrito por Bruno Contreras en Julio de 2005 use strict; use GD::Graph::lines; my @datos = ( [-3,-2,-1,0,1,2,3], # variable X [.2,.4,1.9,2.4,1.2,.3,.3] # variable Y ); # crea un nueva instancia my $grafica = new GD::Graph::lines(); # configura el diagrama $grafica->set( x_label => 'sigma', y_label => 'frecuencia', title => 'Ejemplo de gráfica X/Y' ) or die $grafica->error; my $diagrama = $grafica->plot( \@datos ); open(GRAF_PNG,">ejemplo.png") || die "no puedo crear ejemplo.png\n"; binmode GRAF_PNG; print GRAF_PNG $diagrama->png(); close(GRAF_PNG);
En esta sección veremos como se puede utilizar el poderoso ambiente de programación estadística R, que tantas aplicaciones tiene en el tratamiento estadístico de datos, incluyendo datos biológicos y *ómicos. Para ello debes descargar e instalar el paquete RSPerl, con dos sencillos pasos:
./configure --enable-R-shlib make
tar zxf RSPerl_0.9-0.tar.gz seguido de
R CMD INSTALL --configure-args='--with-in-perl' RSPerl
use lib "/path/to/R-2.5.0/library/RSPerl/perl/i386-linux";
En el siguiente programa muestro un ejemplo:
#!/usr/bin/perl -w
use strict;
use lib "/path/to/R-2.5.0/library/RSPerl/perl/i386-linux/";
use R;
use RReferences;
R::initR("--silent");
R::library("RSPerl");
# input data
my @x = (44.4, 35.9, 41.9, 100.3, 60.7, 84.1, 105.7, 45.2, 60.1);
my @y = ( 2.2, 1.7, 2.0, 5.0, 3.6, 4.0, 5.2, 2.8, 3.8);
# correlation
my $test = "spearman";
my $corr = R::callWithNames("cor.test",{'x'=>\@x,'y'=>\@y,'method'=>"$test"});#pearson,spearman
my $corr_coef = $corr->getEl("estimate"); # hacer edit(corr) en R para ver sus atributos
my $pvalue = $corr->getEl("p.value");
print "# $test : coef=$corr_coef p=$pvalue (h0: coef==0)\n";
# plot scatterplot & regression line
my $plotfile = "corr.png";
my $dataframe = R::callWithNames("data.frame", {'x', \@x, 'y', \@y});
my $linearfit = R::call("lm","y~x",$dataframe);
R::call("png","$plotfile");
R::callWithNames("plot", {'x',\@x,'y',\@y,'ylab'=>'y','xlab'=>'x','main'=>'title'});
R::call("abline",$linearfit);
R::call("dev.off");
Aquí os muestro algunos programas escritos en Perl que os pueden ser útiles. Están escritos en Perl y podéis por lo tanto curiosear en el código.
El programa mview
está escrito
en Perl y sirve para representar en forma de alineamiento múltiple la salida de BLAST.
Por ejemplo, esta salida de
blastpgp
1lfu_nr.blast
podemos convertirla a 1lfu_nr.html.
ps_scan.pl es un programa para hacer búsquedas de motivos PROSITE en secuencias de proteínas.
Hasta ahora hemos visto la potencia de Perl para manipular cadenas de caracteres, que de hecho es para lo que se diseñó el lenguaje. Podemos preguntarnos además qué capacidad tiene el lenguaje para las operaciones aritméticas con números reales, algo que normalmente no se ve en el código, más bien queda escondido en los detalles de la compilación. Para entender la cuestión veamos un ejemplo programado en C y en Perl;
/* Ejemplo de Julio Freyre */
#include <stdio.h>
void main()
{
for(float x = 0; x < 1; x += .1)
printf("%f es menor que 1\n", x);
}
Ahora la versión en Perl:
# este programa muestra el error que acumula perl al sumar reales
for(my $x = 0; $x < 1; $x += .1)
{
print("$x es menor que 1\n");
}
# el error es de 1e-16, este programa ya funciona como se espera en arquitecturas i586 y SPARC
my $error = 1e-16;
for(my $x = 0; $x+$error < 1; $x += .1)
{
print("$x es menor que 1\n");
}
Si probáis estos programitas veréis la limitación de Perl al manejar números reales. ¿Qué se puede hacer para resolver esto? La solución más sencilla es utilizar un módulo de CPAN que se llama Math::BigFloat, que permite usar una precisión arbitaria, elegida por el usuario, para representar reales y operar con ellos, con un coste en el tiempo de ejecución.
Por cierto, hay otro módulo llamado Math::BigInt que permite extender la capacidad intrínseca del lenguaje y la arquitectura física subyacente para representar enteros.
En este capítulo veremos con un poco más de detalle dos paquetes muy importantes, Bioperl y DBI, hechos para facilitarnos la vida en multitud de aplicaciones bioinformáticas y no reinventar la rueda con muchas tareas habituales.
Bioperl es un proyecto comunitario de código abierto que lleva ya más de 10 años produciendo aplicaciones escritas en Perl para la bioinformática, la genómica y las ciencias biológicas en general. El recurso bioinformático más llamativo hecho con Bioperl sea probablemente la base de datos Ensembl.
Según reza la el tutorial de Bioperl:
Bioperl does provide reusable perl modules that facilitate writing perl scripts for sequence manipulation, accessing of databases using a range of data formats and execution and parsing of the results of various molecular biology programs including Blast, clustalw, TCoffee, genscan, ESTscan and HMMER. Consequently, bioperl enables developing scripts that can analyze large quantities of sequence data in ways that are typically difficult or impossible with web based systems.
El proyecto distribuye programas y paquetes Perl (incluidos en CPAN) y tiene su propio tutorial y su documentación. El listado completo de los 500 módulos de Bioperl puedes encontrarlo aquí y ahí verás además documentación y ejemplos de cada módulo.
En este capítulo vamos a recorrer algunas de las utilidades de Bioperl, sobre todo a base de ejemplos y ejercicios.
Podéis descargaros la versión más reciente de
Bioperl. Para escribir estas páginas
yo descargé bioperl-1.5.tar.gz (core) y lo descomprimí con
tar xvfz bioperl-1.5.tar.gz
en un directorio. Una vez hecho esto puedes correr el programa
bptutorial.pl
que te muestra algunas de las cosas que puedes hacer con Bioperl, con estos argumentos:
1 => sequence_manipulations 2 => seqstats_and_seqwords 3 => restriction_and_sigcleave 4 => other_seq_utilities 5 => run_perl 6 => searchio_parsing 7 => bplite_parsing 8 => hmmer_parsing 9 => simplealign 10 => gene_prediction_parsing 11 => access_remote_db 12 => index_local_db 13 => fetch_local_db (NOTE: needs to be run with demo 12) 14 => sequence_annotation 15 => largeseqs 16 => liveseqs 17 => run_struct 18 => demo_variations 19 => demo_xml 20 => run_tree 21 => run_map 22 => run_remoteblast 23 => run_standaloneblast 24 => run_clustalw_tcoffee 25 => run_psw_bl2seq
Una buena manera de iniciarse en las capacidades de Bioperl es revisar el código de bptutorial.pl.
Para instalar el paquete seguid las instrucciones que ya vimos en la sección 3.4.2.
Una de las cosas más útiles de Bioperl es que te permite acceder a bases de datos en la red y a herramientas web desde un programa Perl, permitiendo así hacer análisis bioinformáticos complejos con muy poca infraestructura, tan sólo una conexión a Internet. La única precaución que debes tener es no sobrecargar de forma abusiva los servicios web que utilices, sólo podría perjudicarte.
En este ejemplo os muestro cómo obtener desde la red secuencias de GenBank(Benson et al., 2002) .
Se pueden consultar por esta vía además SwissProt (Boeckmann et al., 2003) , GenPept, EMBL y RefSeq.
#!/usr/bin/perl -w
# Ejemplo donde obtenemos por la red una secuencia de GenBank y la guardamos en formato FASTA
use strict;
use Bio::Perl;
my $gi = "NP_417816";
# ahora una de GenBank
my $secuenciaGenBank = get_sequence( 'genbank', '$gi' );
# escribe la secuencia en formato FASTA a un archivo
write_sequence(">$gi\.fas", 'fasta', $secuenciaGenBank);
Llegados a este punto, después de tanto rollo, creo que lo mejor es proponer un ejercicio semejante a los problemas que tendréis que resolver alguna vez en vuestros proyectos, donde normalmente la única ayuda disponible será Google o estos tutoriales: este, o este otro.
En este ejercicio os voy a pedir el diseño e implementación de un sistema sencillo que intente anotar de forma automática secuencias de proteínas y ácidos nucleicos según son producidas por un proyecto de secuenciación. Anotar es poner en forma de texto una caracterización de su función, de su historia evolutiva y de sus detalles. El sistema deberá:
Por supuesto, mejor si hacéis todo esto usando herramientas de Bioperl. Quizás este tutorial os sea de utilidad.
DBI es un conjunto de módulos Perl que puedes descargar de CPAN si no los tienes ya en tu sistema. DBI significa Database Independent Interface , que en español es una interfaz para bases de datos independiente de la plataforma, similar a la interfaz JDBC del lenguaje java, por ejemplo.
Hay todo un mundo de diferentes conceptos e implementaciones de bases de datos. DBI nos permite interaccionar desde Perl con muchas de ellas, ya sean propietarias o de código abierto como MySQL. Uno de las razones del éxito de DBI es que es sencillo escribir los controladores de un tipo de base de datos dado para que DBI pueda interaccionar con ella.
Las aplicaciones de DBI incluyen sobre todo las bases de datos relacionales, construidas a base de tablas y manejadas normalmente por medio de un lenguaje estándar llamado SQL, del inglés Structured Query Language.
DBI es un conjunto de módulos escritos en forma de clases.
En la sección 1 se citan una serie de documentos donde podéis complementar el aprendizaje y los usos de DBI que veamos aquí.
![\begin{figure}
\begin{center}
\includegraphics[width=0.8\textwidth]{DBIflux}
\end{center}
\end{figure}](img7.png) |
El primer paso para consultar una base de datos por medio de DBI será conectarse a su
servidor desde nuestro programa Perl, normalmente a través de la red. Para hacerlo debemos crear
una instancia de la clase DBI y usar su método
connect .
La información necesaria para conectarse es:
![\begin{figure}
\begin{center}
\includegraphics[width=0.8\textwidth]{BDoperaciones}
\end{center}
\end{figure}](img8.png)
dbi:mysql:test:servidor.ccg.unam.mx
Si vas a conectarte a un servidor de base de datos Oracle necesitarás el controlador DBD::Oracle,
de la misma manera que necesitarás DBD::mysql para conectarte a una máquina que corra mysql.
Puedes descargar de CPAN
los controladores que necesites,
pero antes comprueba qué controladores hay disponibles en tu sistema con:
$ perl -e 'use DBI; print DBI->available_drivers();'
Aquí muestro un ejemplo de código de conexion a la base de datos
test
del servidor MySQL
servidor.ccg.unam.mx :
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Agosto de 2005
use strict;
use DBI;
my $conexion = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave");
if(!$conexion)
{
print "No me puedo conectar a mysql:test:servidor.ccg.unam.mx\n$DBI::errstr\n";
exit(-2);
}
$conexion->disconnect() || warn "Fallo al desconectarme de mysql:test:servidor.ccg.unam.mx\n$DBI::errstr\n";
En principio puedes mantener varias conexiones a la misma base de datos, y cada una de ellas se comportará de forma independiente:
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Agosto de 2005
use strict;
use DBI;
my $conexion1 = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave");
if(!$conexion1)
{
print "Falló conexión1 a mysql:test:servidor.ccg.unam.mx\n$DBI::errstr\n";
exit(-2);
}
my $conexion2 = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave");
if(!$conexion2)
{
print "Falló conexión2 a mysql:test:servidor.ccg.unam.mx\n$DBI::errstr\n";
exit(-3);
}
$conexion1->disconnect() || warn "Falla desconexión1 de mysql:test:servidor.ccg.unam.mx\n$DBI::errstr\n";
$conexion2->disconnect() || warn "Falla desconexión2 de mysql:test:servidor.ccg.unam.mx\n$DBI::errstr\n";
A veces es difícil conectarse a servidores muy cargados de trabajo, por lo que podemos probar a conectarnos varias veces, hasta que lo consigamos:
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Agosto de 2005
use strict;
use DBI;
my $conexion;
until( $conexion = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave") )
{
warn "No consigo conectarme a a mysql:test:servidor.ccg.unam.mx\nLo intento de nuevo en 1 minuto\n";
sleep(60);
}
$conexion->disconnect() || warn "Fallo al desconectarme de mysql:test:servidor.ccg.unam.mx\n$DBI::errstr\n";
Una vez que hemos establecido una conexión con una base de datos podemos utilizar el objeto DBI creado para consultarla y manipularla por medio de comandos SQL, siempre y cuando tengamos los permisos necesarios.
Os recomiendo que miréis el libro Programming the Perl DBI para ver un repertorio más amplio de formas de interacción con bases de datos.
Usaremos el método
DBI::prepare()
para crear un objeto manipulador
(statement handler) y luego el método
execute()
de este nuevo objeto para efectivamente manipular la base de datos.
#!/usr/bin/perl -w
# Ejemplo escrito por Bruno Contreras en Agosto de 2005
use strict;
use DBI;
## por Dios, trata de conectarlo, Carlos
my $conexion = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave",{
PrintError => 0, ### avisa si hay errores mediante warn
RaiseError => 1 ### avisa y termina si hay errores mediante die
}
);
## ahora manipula la base de datos mediante consultas SQL
my $tabla = "blast";
my $consulta="SELECT hit FROM $tabla WHERE evalue < 1e-03";
my $manipulador = $conexion->prepare($consulta); # normalmente, sólo devuelve un objeto
# manipulador si la consulta es válida
$manipulador->execute();
while (my @datos=$manipulador->fetchrow_array) # puedes usar fetchrow_hash con columnas como claves
{
# haz algo con estos datos, como imprimirlos
foreach my $dato (@datos)
{
print "$dato,";
}
}
$manipulador->finish();
$conexion->disconnect();
Otra manera útil de capturar la salida de una consulta sería:
use strict;
use DBI;
my $conexion = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave",{PrintError,0,RaiseError,1});
my $consulta="SELECT * FROM blast WHERE evalue < 1e-03";
my $manipulador = $conexion->prepare($consulta);
$manipulador->execute();
my $ref_arreglo = $manipulador->fetchall_arrayref(); # referencia a un arreglo que contiene
# los resultados
$conexion->disconnect();
Otra manera alternativa de mostrar los resultados de una consulta SQL SELECT podría ser:
use strict;
use DBI;
my $conexion = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave",{PrintError,0,RaiseError,1});
## ahora manipula la base de datos mediante consultas SQL
my $tabla = "blast";
my $consulta="SELECT hit FROM $tabla WHERE evalue < 1e-03";
my $manipulador = $conexion->prepare($consulta);
$manipulador->execute();
# imprime todos los resultados obtenidos de la consulta
$manipulador->dump_results();
# guarda los resultados en un archivo en formato de 80 columnas
$manipulador->execute();
open(CONSULTA,">resultados.txt") || die "no puedo crear resultados.txt\n";
$manipulador->dump_results(80, "\n", ':', \*CONSULTA);
$conexion->disconnect();
Finalmente os muestro una forma de hacer tus programas con DBI más comprensibles, usando el
método bind_columns:
use strict;
use DBI;
my $conexion = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave",{PrintError,0,RaiseError,1});
## ahora manipula la base de datos mediante consultas SQL
my ($hit,$evalue); # columnas de $tabla
my $tabla = "blast";
my $consulta="SELECT hit,evalue FROM $tabla WHERE evalue < 1e-03";
my $manipulador = $conexion->prepare($consulta);
$manipulador->execute();
# asocia columnas de salida con sus nombres
$manipulador->bind_col( 1, \$hit );
$manipulador->bind_col( 2, \$evalue );
# procesa la salida
while ($manipulador->fetch())
{
print "$hit $evalue\n";
}
$conexion->disconnect();
Para manipular variables de tamaño ilimitado en tablas (tipos LONG y LOB) hay que hacer consultas SELECT modificadas, que podéis ver aquí.
Para otros comandos SQL, como INSERT, UPDATE, DROP o CREATE, podemos seguir usando la combinación
prepare() +
execute() ,
podemos hacerlo también en un sólo paso, con el método
do() .
Podéis usarlo para sentencias como:
DELETE FROM blast WHERE hit = 'ECOLI_CRP';
UPDATE tabla SET var1 = var1+1 WHERE id = 247;
INSERT INTO variables VALUES ('grafo', 0);
Éste es un ejemplo de uso de
do() :
use strict;
use DBI;
my $conexion = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave",{PrintError,0,RaiseError,1});
# filas contiene el número de filas afectadas por la sentencia,
# en caso de error devuelve falso
my $filas = $conexion->do("DELETE FROM blast");
if(!$filas)
{
print "no pude ejecutar la expresión SQL\n";
}
$conexion->disconnect();
Para las bases de datos basadas en transacciones estas sentencias sólo tendrán efecto si tras cursarlas se ejecuta el comando COMMIT, que por defecto DBI ejecuta automáticamente por nosotros. Pero también podemos hacerlo explícitamente:
use strict;
use DBI;
my $conexion = DBI->connect("dbi:mysql:test:servidor.ccg.unam.mx","usuario","clave",{PrintError,0,RaiseError,1,AutoCommit,0});
# filas contiene el número de filas afectadas por la sentencia,
# en caso de error devuelve falso
my $filas = $conexion->do("INSERT INTO variables VALUES ('grafo', 0)");
if(!$filas)
{
print "no pude ejecutar la expresión SQL\n";
}
else
{
# si quieres fijar estos cambios:
$conexion->commit();
# si quisieras descartalos:
$conexion->rollback();
}
$dbh->rollback();
$conexion->disconnect();
Ahora os toca trabajar a vosotros. Podréis usar los ejemplos vistos hasta ahora y esta fuente adicional como inspiración para estos ejercicios:
Con la excusa de ir aprendiendo diferentes aspectos de Perl, hemos ido viendo algunas de sus aplicaciones en bioinformática, pero nos hemos dejado algunas muy importantes fuera. Creo que este curso se quedaría cojo si no viéramos al menos brevemente temas como:
Vamos a ver ahora ejemplos de cómo Perl puede ser muy útil en estos campos. Os pido que me escribáis (mirar sección 1) si creéis que me dejado fuera alguna aplicación de importancia, trataré de incluirla para la siguienre versión.
Dentro del mundo Open Source, la mayoría de los paquetes para tratamiento estadístico de datos de microarreglos se basan en la plataforma R y la biblioteca de paquetes Bioconductor, que no tienen mucho que ver con Perl.
Pero sí hay otros paquetes escritos en Perl para estos fines, como GP3, AMAD: Another Microarray Database, Statistics tools for gene expression o cellwall.stanford.edu/scripts/. Google os dará muchos más, así que, como otras veces, vale la pena ver por ahí que hay hecho antes de lanzarnos a escribir nosotros programas que ya existen.
Cuando hablamos de microarreglos normalmente hablamos de matrices muy grandes donde cada fila normalmente se corresponde con un gen. En Perl es muy sencillo manejar este tipo de datos, como veremos en el siguiente ejemplo.
Vamos a utilizar los datos de Gutiérrez-Ríos et al. (2004), accesibles desde aquí, para este ejercicio. La idea es que os copiéis el archivo plano con los datos a vuestra máquina local y os pido que escribáis un programa que:
menE, cyoA y flgL, seleccione mediante coeficientes
de correlación los genes cuyo patrón de expresión se parezca. Os sugiero que uséis el
módulo
Stats.
Ésta es una de las aplicaciones más prometedoras y quizás menos conocidas de la bioinformática. ¿Pero qué es la minería de datos? En pocas palabras es la búsqueda de conocimiento enterrado en la literatura de forma automatizada. Es una herramienta muy útil para lograr un mejor aprovechamiento de la literatura científica, permitiéndonos acceder a conceptos y contenidos que están ya descritos en las publicaciones científicas, sin ser necesariamente especialistas en ese área. Os recomiendo esta presentación para profundizar sobre este tema. Para haceros una idea de qué se está haciendo ahora mismo en el campo podéis mirar este artículo de Müller et al. (2004).
Hay multitud de herramientas en la web para hacer minería de datos, de alcance general o especializadas en biología, pero aquí vamos a ver cómo aproximarnos a estos problemas desde Perl.
En este ejercicio vamos a trabajar con dos conjuntos de datos separados.
El primero consta de 426 resúmenes o abstracts de artículos científicos devueltos por PubMed con la consulta 'CRP + coli + regulation', que podéis encontrar en el archivo CRP_coli_regulation_PubMed.txt.
El segundo es este artículo completo gb-2000-1-2-research0004.xml en formato XML, estraído de BioMedCentral, una de las editoriales que están publicando artículos científicos de libre acceso.
Vuestro trabajo consiste en escribir código Perl que:
La detección de genes es quizás la tarea más importante del procesamiento de los datos generados por un proyecto de secuenciación. De hecho, las predicciones de genes, correctas o equivocadas, son la base para casi todo el trabajo experimental que se haga después de alcanzar una versión medianamente estable del genoma.
La tarea de localizar genes es complicada y muy costosa, sobre todo por el tamaño de los datos implicados y por la dificultad de distinguir la señal del ruido acumulado en los genomas. En eucariotas se complica todavía más por los enormes genomas que tienen, por la cantidad de secuencias basura que aparentemente no tienen una función biológica y por la presencia de intrones. En otras palabras, es todo un mundo al que podéis iniciaros en www.genefinding.org. Por lo que he visto, Perl no se utiliza demasiado en este campo por razones de eficiencia. Dado el coste computacional del problema, se prefieren lenguajes como C/C++.
Sin embargo, para efectos de este curso, podemos hacer una aproximación un tanto basta al problema desde Perl, aprovechando sus fortalezas en el tratamiento de cadenas. El algoritmo que os propongo consiste en 4 pasos:
Hemos llegado ya a la parte donde os toca aplicar todo lo aprendido previamente para realizar las tareas que se han asignado a cada grupo. Esto será trabajo de equipo y así será evaluado.
Los datos crudos con los que tendréis que trabajar podéis encontrarlos
en el directorio del curso
en /home/compu2/.
Se os pide que diseñéis e implementéis programas en lenguaje Perl para interaccionar con vuestra base de datos, la de vuestro grupo o proyecto, e ir poblándola con los datos necesarios, según determinásteis en el módulo anterior de vuestro proyecto, tras todo el procesamiento que sea necesario.
Como véis no son instrucciones muy concretas, por lo que tendréis que ir justificando por escrito vuestro trabajo como grupo. Os recomiendo que consultéis readyset para redactar vuestros documentos. Aquí hallaréis plantillas.
La idea de todo esto es enseñaros un concepto nuevo, que no hemos visto antes en el curso: que normalmente vale la pena hacer un poco de planificación antes de ponerse a escribir código. Muchos programadores hacen todo esto mentalmente, sin documentos, si acaso algunas notas. Pero cuando se trata de proyectos en equipo o de cierta envergadura, entonces se hace más necesario recurrir a los documentos.
Por último, cuando vayáis avanzando en vuestro proyecto, debéis tratar de que todo, documentos y código, sean lo más sencillo posible. Que no os importe si los requisitos o las especificaciones ocupan media hoja, si tan pequeño es el proyecto.
Nuestro ciclo de desarrollo particular seguirá este flujo:
Aquí iré poniendo enlaces relacionados con este curso y los que me vayáis enviando vosotros:
Certamen de poesía en Perl
An introduction to bioinformatics algorithms
This document was generated using the LaTeX2HTML translator Version 2002-2-1 (1.70)
Copyright © 1993, 1994, 1995, 1996,
Nikos Drakos,
Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999,
Ross Moore,
Mathematics Department, Macquarie University, Sydney.
The command line arguments were:
latex2html bioinfoPerl -split 0 -dir bioinfoPerl1 -no_navigation
The translation was initiated by Bruno Contreras Moreira on 2007-06-15
![\begin{figure}
\begin{center}
\includegraphics[width=0.8\textwidth]{cicloD}
\end{center}
\end{figure}](img9.png)
