27 de diciembre de 2010

Enviar y recibir datos de un formulario con PHP

Por mucho que nos guste programar en Perl, a veces hay que recurrir a otros lenguajes más adecuados para propósitos concretos. Este es el caso de PHP cuando se trata de realizar interfaces web para nuestros scripts de Perl.
El primer problema que nos surge al pasar a PHP es cómo enviar y recibir datos a través de un formulario web. Los datos enviados los podemos usar para ejecutar en el servidor nuestro script de Perl y generar resultados.
A continuación muestro el ejemplo de un formulario de búsqueda de secuencias DNA o Proteínas:

 <?php  
 // Shows a search input data formulary  
 function showSearchForm(){  
      ?>  
      <h3>Search form:</h3></br>  
      <form enctype="multipart/form-data" action="<?php print $_SERVER['PHP_SELF'] ?>" method="POST">  
      Search name: <input type="text" name="search_name" size="40" maxlength="127"></br></br>  
      Email: <input type="text" name="user_email" size="40" maxlength="127"></br></br>  
      Search input type:  
      <input type="radio" name="search_type" value="stamp"> DNA  
      <input type="radio" name="search_type" value="blastp"> Protein  
      </br></br>  
      Exclude query from results <input type="checkbox" name="search_exclude" value="yes"></br></br>  
      Limit number of results per query: <input type="text" name="search_limit" size="2" maxlength="2"></br></br>  
      Query data or file:</br></br>  
      <textarea name="query_data" rows="10" cols="67"></textarea>  
      </br></br>  
      <input type="file" name="query_file" size="20"></br></br>  
      Genomes:   
      <select name="specie" class="selectbutton">  
      <option value="human">Homo sapiens</option>  
      <option value="mouse">Mus musculus</option>  
      <option value="arabidopsis">Arabidopsis thaliana</option>  
      <option value="rice">Oryza sativa</option>  
      </select>  
      </br></br>  
      <input type="submit" name="search_submitted" value="Search">  
      </form>  
      <?php  
 }  
 ?>  

El siguiente paso consistirá en validar los datos introducidos para volver a mostrar el formulario si los datos fueran erróneos. Usaremos la variable especial $_SESSION para guardar los datos enviados con el formulario (almacenados en la variable $_POST). La variable $_SESSION nos permitirá acceder a su contenido si ejecutamos 'session_start()' al comienzo del código PHP. Antes de guardar los datos enviados, borraremos los contenidos anteriores de $_SESSION. Si en vez de enviar un archivo, hemos enviado datos en formato texto, guardaremos estos datos en un archivo temporal ($_SESSION['query_file']). Si faltan datos en algún campo del formulario, crearemos un mensaje de error en $_SESSION['errors'].

 <?php  
 // Validates search form data  
 function validateSearchData(){  
      // Stores in $_SESSION variable all the post data and in $_SESSION['errors'] if some important data is missing  
      // Clear previous $_SESSION values  
      unset($_SESSION['errors']);  
      unset($_SESSION['search_name']);  
      unset($_SESSION['user_email']);  
      unset($_SESSION['search_type']);  
      unset($_SESSION['search_exclude']);  
      unset($_SESSION['search_limit']);  
      unset($_SESSION['specie']);  
      unset($_SESSION['output_dir']);  
      unset($_SESSION['query_file']);  
      unset($_SESSION['query_data']);  
      (!empty($_POST['search_name'])) ? $_SESSION['search_name'] = trim($_POST['search_name']) : $_SESSION['errors'][] = 'name for the search';  
      (!empty($_POST['user_email'])) ? $_SESSION['user_email'] = trim($_POST['user_email']) : $_SESSION['errors'][] = 'email';;  
      (!empty($_POST['search_type'])) ? $_SESSION['search_type'] = $_POST['search_type'] : $_SESSION['errors'][] = 'search input type';  
      if (!empty($_POST['search_exclude'])){  
           ($_POST['search_exclude'] == 'yes') ? $_SESSION['search_exclude'] = 1 : $_SESSION['search_exclude'] = 0;  
      }  
      (!empty($_POST['search_limit'])) ? $_SESSION['search_limit'] = trim($_POST['search_limit']) : '';  
      (!empty($_POST['specie'])) ? $_SESSION['specie'] = $_POST['specie'] : '';  
      // Store the query data into a file inside the /tmp folder  
      if (!empty($_FILES['query_file']['name']) && !empty($_FILES['query_file']['tmp_name']) && $_FILES['query_file']['size']>0){  
           $_SESSION['query_file'] = $_FILES['query_file']['tmp_name'];  
      }elseif (!empty($_POST['query_data'])){  
           $_SESSION['query_data'] = $_POST['query_data'];  
           $_SESSION['query_file'] = '/tmp/'.$_SESSION['search_name'].".data";  
           $file_handle = fopen($_SESSION['query_file'], "w");  
           fwrite($file_handle, $_POST['query_data']);  
           fclose($file_handle);  
      } else{  
           $_SESSION['errors'][] = 'data or file to process';  
      }  
      if (isset($_SESSION['errors'])){  
           return;  
      }  
 }  
 ?>  

Finalmente, falta introducir las dos funciones anteriores en un archivo (por ejemplo 'index.php') para que sean leídas e interpretadas por el servidor web. Si los datos del formulario son válidos se ejecutará un pequeño código de Perl, si no se volverá a mostrar el formulario indicando los datos que faltan.

 <?php  
 // Example file of a PHP search formulary with input data validation and script execution  
 session_start();  
 ?>  
 <html>  
 <head>  
 <title>Search Form</title>  
 </head>  
 <body>  
 <?php  
 if(!isset($_POST['search_submitted'])){  
      ?><h3>Search form:</h3><?php  
      showSearchForm();  
 } else {  
      validateSearchData();  
      if (!isset($_SESSION['errors'])){  
           executeScript();  
      } else {  
           ?><h3>Search form:</h3><?php  
           print '<font color="red"><b>Please specify: '.join(', ',$_SESSION['errors']).'.</b></font></br></br>';  
           unset($_SESSION['errors']);  
           showSearchForm();  
      }  
 }  
 ?>  
 </body>  
 </html>  
 <?php  
 // Shows an input data formulary filled with previous values  
 function showSearchForm(){  
      ?>  
      <form enctype="multipart/form-data" action="<?php print $_SERVER['PHP_SELF'] ?>" method="POST">  
      Search name: <input type="text" name="search_name" size="40" maxlength="127" <?php ($_SESSION['search_name']) ? print 'value="'.$_SESSION['search_name'].'"' : '' ?> ></br></br>  
      Email: <input type="text" name="user_email" size="40" maxlength="127" <?php ($_SESSION['user_email']) ? print 'value="'.$_SESSION['user_email'].'"' : '' ?> ></br></br>  
      Search input type:  
      <input type="radio" name="search_type" value="stamp" <?php ($_SESSION['search_type'] == 'stamp') ? print 'checked' : '' ?>> DNA  
      <input type="radio" name="search_type" value="blastp" <?php ($_SESSION['search_type'] == 'blastp') ? print 'checked' : '' ?>> Protein  
      </br></br>  
      Exclude query from results <input type="checkbox" name="search_exclude" value="yes" <?php ($_SESSION['search_exclude'] == 'yes') ? print 'checked' : print 'checked' ?>></br></br>  
      Limit number of results per query: <input type="text" name="search_limit" size="2" maxlength="2" <?php ($_SESSION['search_limit']) ? print 'value="'.$_SESSION['search_limit'].'"' : print 'value="10"'; ?> ></br></br>  
      Query data or file:</br></br>  
      <textarea name="query_data" rows="10" cols="67"><?php ($_SESSION['query_data']) ? print $_SESSION['query_data'] : '' ?></textarea>  
      </br></br>  
      <input type="file" name="query_file" size="20"></br></br>  
      Genomes:   
      <select name="specie" class="selectbutton">  
      <option value="human" <?php (!isset($_SESSION['specie']) || $_SESSION['specie'] == 'human') ? print 'selected' : '' ?>>Homo sapiens</option>  
      <option value="mouse" <?php ($_SESSION['specie'] == 'mouse') ? print 'selected' : '' ?>>Mus musculus</option>  
      <option value="arabidopsis" <?php ($_SESSION['specie'] == 'arabidopsis') ? print 'selected' : '' ?>>Arabidopsis thaliana</option>  
      <option value="rice" <?php ($_SESSION['specie'] == 'rice') ? print 'selected' : '' ?>>Oryza sativa</option>  
      </select>  
      </br></br>  
      <input type="submit" name="search_submitted" value="Search">  
      </form>  
      <?php  
 }  
 // Validates search form data  
 function validateSearchData(){  
      // Stores in $_SESSION variable all the post data and in $_SESSION['errors'] if some important data is missing  
      // Clear previous $_SESSION values  
      unset($_SESSION['errors']);  
      unset($_SESSION['search_name']);  
      unset($_SESSION['user_email']);  
      unset($_SESSION['search_type']);  
      unset($_SESSION['search_exclude']);  
      unset($_SESSION['search_limit']);  
      unset($_SESSION['specie']);  
      unset($_SESSION['output_dir']);  
      unset($_SESSION['query_file']);  
      unset($_SESSION['query_data']);  
      (!empty($_POST['search_name'])) ? $_SESSION['search_name'] = trim($_POST['search_name']) : $_SESSION['errors'][] = 'name for the search';  
      (!empty($_POST['user_email'])) ? $_SESSION['user_email'] = trim($_POST['user_email']) : $_SESSION['errors'][] = 'email';;  
      (!empty($_POST['search_type'])) ? $_SESSION['search_type'] = $_POST['search_type'] : $_SESSION['errors'][] = 'search input type';  
      if (!empty($_POST['search_exclude'])){  
           ($_POST['search_exclude'] == 'yes') ? $_SESSION['search_exclude'] = 1 : $_SESSION['search_exclude'] = 0;  
      }  
      (!empty($_POST['search_limit'])) ? $_SESSION['search_limit'] = trim($_POST['search_limit']) : '';  
      (!empty($_POST['specie'])) ? $_SESSION['specie'] = $_POST['specie'] : '';  
      // Store the query data into a file inside the /tmp folder  
      if (!empty($_FILES['query_file']['name']) && !empty($_FILES['query_file']['tmp_name']) && $_FILES['query_file']['size']>0){  
           $_SESSION['query_file'] = $_FILES['query_file']['tmp_name'];  
      }elseif (!empty($_POST['query_data'])){  
           $_SESSION['query_data'] = $_POST['query_data'];  
           $_SESSION['query_file'] = '/tmp/'.$_SESSION['search_name'].".data";  
           $file_handle = fopen($_SESSION['query_file'], "w");  
           fwrite($file_handle, $_POST['query_data']);  
           fclose($file_handle);  
      } else{  
           $_SESSION['errors'][] = 'data or file to process';  
      }  
      if (isset($_SESSION['errors'])){  
           return;  
      }  
 }  
 // Executes a Perl script and prints the output  
 function executeScript(){  
      $output = `perl -e 'print "Hello world\n"'`;  
      print '<pre>'.$output.'</pre>';  
 }  
 ?>  

9 de diciembre de 2010

Compresión de secuencias de ADN

Un problema que me ha surgido recientemente es el de manejar en memoria RAM grandes volúmenes de secuencias de ADN, en mi caso genes, pero que podrían ser también lecturas de secuenciación masiva. Cuando el tamaño total de las secuencias supera la RAM disponible el sistema operativo se verá forzado a hacer swapping/intercambio con el disco y el proceso se ralentiza considerablemente. Una de las posibles soluciones, la que nos ocupa hoy, es la compresión de las secuencias, de manera que ocupen menos espacio físico en RAM.
El siguiente programa muestra cómo podemos realizar esto en Perl, y nos sugiere que esta estrategia vale realmente la pena a medida que el tamaño de las secuencias se aproxima a las 200 bases:




 use strict;  
 use Compress::Zlib;  
   
 # http://perldoc.perl.org/Compress/Zlib.html  
   
 my (@S,$source,$compressed,$uncompressed);  
   
 print "size(source)\tsize(compressed)\tratio\n";  
 for(my $size=25;$size<=500;$size += 25)  
 {  
    # secuencia aleatoria de longitud $size  
    @S = map { ('A','C','G','T')[rand(4)] } 1 .. $size;  
    $source = join(/ /,@S); #print "$source\n";  
   
    $compressed = compress($source);  
    $uncompressed = uncompress($compressed) ;  
      
    printf("%d\t%d\t%g\n",  
       $size,length($compressed),length($compressed)/$size);  
 }  
   

Un saludo,
Bruno

26 de noviembre de 2010

Interfaz SOAP de 3D-footprint

En las recientes Jornadas de Bioinformática 2010 (en Málaga) di una pequeña charla
sobre nuestra base de datos 3D-footprint, con el mismo título del artículo donde la publicamos en NAR "3D-footprint: a database for the structural analysis of protein-DNA complexes" (el  PDF del material presentado está disponible aquí, el vídeo  de la disección de interfaces proteína-DNA aquí).

Logo footprint de 2bam_AB

Tras la presentación me preguntaron si 3D-footprint tenía una interfaz de servicios web. La respuesta es afirmativa, en particular es una interfaz SOAP que se puede usar fácilmente desde Perl, como se muestra en el siguiente ejemplo:

 #!/usr/bin/perl -w  
   
 # prog 10 Bruno Contreras-Moreira  
 # web services client for 3D-footprint  
 # ejemplo basado en http://floresta.eead.csic.es/3dfootprint/tutorial.html#ws  
   
 use strict;  
 use SOAP::Lite;  
   
 my $URI  = 'http://floresta.eead.csic.es/footprint';  
 my $WSURL = 'http://floresta.eead.csic.es/3dpwm/scripts/server/ws.cgi';  
   
 my $soap = SOAP::Lite  
 -> uri($URI)  
 -> proxy($WSURL);  
   
 # 1) consulta de secuencia proteica (de un posible factor de transcripción)  
 my $secuencia =   
 "IYNLSRRFAQRGFSPREFRLTMTRGDIGNYLGLTVETISRLLGRFQKSGMLAVKGKYITIENNDALAQLAGHTRNVA";  
   
 my $result = $soap->protein_query($secuencia); # 1 param: 1(string) sequence  
   
 unless($result->fault){ print $result->result(); }   
 else{ print 'error: ' . join(', ',$result->faultcode,$result->faultstring); }  
   
 # los resultados muestran los posibles residuos de interfaz proteína-DNA  
 # en mayúsculas  
   
   
 # 2) consulta de motivo de DNA  
 my $motivo = 'tntga';  
 $result = $soap->motif_query($motivo,'local',1);   
 # 3 params: 1(string) motif in IUPAC|CONSENSUS|TRANSFAC format,   
 # 2(string) local|global,   
 # 3(float) evalue cutoff  
   
 unless($result->fault){ print $result->result(); }   
 else{ print 'error: ' . join(', ',$result->faultcode,$result->faultstring); }  
   
 # los resultados muestran complejos proteína-DNA cuyos motivos de DNA reconocidos  
 # se parecen a éste  

19 de noviembre de 2010

Escapar caracteres especiales en un texto antes de utilizar expresiones regulares

Hoy me he encontrado con un pequeño bug del script 'pfam_scan.pl' que sirve para buscar diferentes secuencias de proteínas en una base de datos de dominos de  Pfam (HMMs), todo ello usando la herramienta 'hmmscan' de HMMER3.

El error simplemente consistía en que una subrutina insertaba variables dentro de una expresión regular sin haber escapado los caracteres especiales previamente, con lo cual se insertaban directamente caracteres del tipo '(', '/', ... directamente dentro de la expresión regular y no funcionaba como era deseado.

Para evitarlo basta con escapar dichos caracteres especiales insertando el texto entre '\Q' y '\E' como bien nos ha sugerido Joaquín (más información en perlop):

 my $text_with_special_chars = 'abc$%()/def';  
 my $start_of_text_to_search = '$%';  
   
 $text_with_special_chars =~ /(\Q$start_of_text_to_search\E.+)/;  
   
 print $1; # $1 = '$%()/def'  

15 de noviembre de 2010

Mezclando datos en un solo archivo Excel

A menudo me ha tocado intercambiar datos con colegas en formato Excel y uno de los problemas que he padecido es el de tener que ir juntando a mano diferentes conjuntos de datos en diferentes hojas del mismo libro de cálculo. Una molestia añadida es la necesidad de importar archivos en formato texto, como columnas separadas por tabuladores, por ejemplo, en el mismo libro. Para automatizar este tipo de tareas y usar Excel/OpenOffice realmente sólo para el análisis de los datos, podemos echar mano por ejemplo del módulo Spreadsheet del CPAN. El mismo módulo puede servir además para leer y procesar datos directamente desde archivos Excel.

Volviendo al problema descrito arriba, supongamos por ejemplo que queremos unir en un sólo archivo excel tres archivos (datos1.xls, datos2.xls y datos3.tab) en uno solo que llamaremos todos_datos.xls . Pues bien, por medio del programa descrito debajo sería tan sencillo como ejecutar en el terminal:

 $ ./my_xls_merge.pl todos_datos.xls datos1.xls datos2.xls datos3.tab

Yo lo he probado en Linux y los archivos generados se leen perfectamente con el MS Excel, con la única pega de que se pierden las gráficas que hubiera ya en los archivos de entrada.

El código es el siguiente:

 #!/usr/bin/perl -w  
   
 # programa 9 my_xls_merge.pl   
 # programita para unir hojas de diferentes archivos .xls o .tab en un solo libro excel
 # OJO: conserva el valor de cada celda de una hoja de calculo, no las formulas,   
 # y se pierden los diagramas por ejemplo  
 # inspirado por http://www.perlmonks.org/?node_id=743574  
   
 use strict;  
 use Spreadsheet::ParseExcel;  
 use Spreadsheet::WriteExcel;   
 use File::Basename;  
   
 my $MAXLENGTHSHEETNAME = 31;  
   
 die "# usage: $0 <outfile.xls> <file1.xls> ... <fileN.xls>\n" if(scalar(@ARGV) < 2);  
   
 my @infiles = @ARGV;  
 my $outfile = shift(@infiles);  
   
 if(-e $outfile)  
 {  
    print "# a file named '$outfile' already exists, overwrite? (Y/N)\n";  
    my $overwrite = <STDIN>;  
    if($overwrite ne "Y\n"){ exit }  
    else{ unlink($outfile) }  
 }  
   
 # instantiate resulting excel object  
 my $outbook = Spreadsheet::WriteExcel->new($outfile) ||   
    die "# $0: cannot create '$outfile': $!";  
   
 for my $file (@infiles)   
 {  
    if(!-e $file)  
    {  
       print "# skipping non-existing file '$file'\n";  
       next;  
    }  
   
    # remove path   
    my ($filename,$copysheet_name) = (fileparse($file));   
      
    my $excel = Spreadsheet::ParseExcel::Workbook->Parse($file);  
    if(defined($excel->{'File'})) # excel format infiles  
    {  
       foreach my $sheet (@{$excel->{Worksheet}})   
       {   
          $copysheet_name = $sheet->{'Name'} . '-' . $filename;  
          if(length($copysheet_name) > $MAXLENGTHSHEETNAME)  
          {  
             print "# abbreviating $copysheet_name to ";  
             $copysheet_name = substr($copysheet_name,0,$MAXLENGTHSHEETNAME-3).'...';  
             print "$copysheet_name\n";  
          }  
         print "# adding '$file' sheet '$sheet->{'Name'}' ($copysheet_name)\n";  
   
        my $copysheet = $outbook->add_worksheet($copysheet_name);  
            
          $sheet->{'MaxRow'} ||= $sheet->{'MinRow'};  
          foreach my $row ($sheet->{'MinRow'} .. $sheet->{'MaxRow'})   
          {  
               my @rowdata = map { $sheet->{'Cells'}->[$row]->[$_]->{'Val'} }   
                      $sheet->{'MinCol'} .. $sheet->{'MaxCol'};  
          $copysheet->write($row,0,\@rowdata);   
          }  
       }  
    }  
    elsif($file =~ /\.tab/)# .tab infiles, seria trivial manejar .csv por ejemplo  
    {  
       $copysheet_name = $filename;  
       print "# adding '$file' no_sheet ($copysheet_name)\n";  
         
       my $copysheet = $outbook->add_worksheet($copysheet_name);  
         
       open(INTAB,$file) || die "# $0 : cannot read '$file': $!";  
       while(<INTAB>)  
       {  
          chomp($_);  
          my @tabdata = split(/\t/,$_);  
          $copysheet->write($.-1,0,\@tabdata);  
       }  
       close(INTAB);  
    }  
 }  
   
 $outbook->close();  
 print "# outfile: $outfile\n";   
   

Un saludo

9 de noviembre de 2010

Semana de la Ciencia 2010

Ayer comenzó la Semana de la Ciencia 2010, el mayor evento de comunicación social de la ciencia y la tecnología que se celebra en España, organizada por el Ministerio de Ciencia e Innovación.

Museos, Universidades, centros de investigación, parques tecnológicos o empresas organizan exposiciones, cursos, visitas, talleres, mesas redondas, excursiones o conferencias, acercando al público en general su quehacer diario, tanto sus aspectos más llamativos como los más desconocidos.

En concreto, la Agencia Estatal Consejo Superior de Investigaciones Científicas (CSIC) a la que pertenece nuestro centro de investigación, la Estación Experimental de Aula Dei (EEAD), organiza diversas actividades de divulgación en todo el país durante esta semana que se pueden consultar en su página web: http://www.semanadelaciencia.csic.es/

En concreto, en Aragón, una de las actividades más importantes se realiza en la Delegación del CSIC en Aragón del 2 al 10 de Noviembre.  Participan la Estación Experimental de Aula Dei, Instituto Pirenaico de Ecología, Instituto de Carboquímica, Instituto de Ciencias de Materiales de Aragón, y el Laboratorio de Investigación en Tecnologías de la Combustión. Cada participante dispone de un Stand con capacidad suficiente para preparación de experimentos atractivos, talleres didácticos, exposición de productos vistosos y utilización de demostradores para explicar conceptos de ciencia en los que se basan sus líneas de investigación.

8 de noviembre de 2010

Buscar un elemento en un array de perl

En Perl no tenemos la útil función "in_array()" que posee PHP, por lo cual cuando queremos saber si un array contiene un determinado elemento tenemos que escribir un poquito más de código.

Aunque buscando en google "perl array search element" obtenemos muchas soluciones a este problema, aquí presentaré 2 muy sencillas.

Para saber únicamente si un elemento está o no está en el array:

 # Search an element inside an array and return 1 or 0 if the element is present or not, example: in_array(\@my_array, $element);  
 sub in_array {  
      my ($array,$search_for) = @_;  
      my %items = map {$_ => 1} @$array; # create a hash out of the array values  
      my $exist = (exists($items{$search_for}))?1:0;  
      return $exist;
 }  

Si queremos además poder buscar expresiones regulares en todo el array y
conocer la o las posiciones de los resultados:

 # Search a regular expression inside an array and retrieve the positions of the hits, example: in_array(\@my_array, "/$expression/");  
 sub in_array {  
      my ($array,$search_for) = @_;  
      my @pos;  
      if ($search_for =~ /^\/.+\/$/){  
           $search_for =~ s/^\/|\/$//g;  
           @pos = grep { $_ =~ /$search_for/ } @{$array};  
      } else {  
           @pos = grep { $_ eq $search_for } @{$array};  
      }  
      (scalar @pos)?return (1,\@pos):return (0);  
 }  

Otra opción sería usar el módulo Tie::IxHash que permite crear hashes manteniendo el orden de los elementos y buscar los índices de los mismos:

 # Search a regular expression inside an array and retrieve the positions of the hits, example: in_array(\@my_array, "/$expression/");  
 sub in_array {  
      my ($array,$search_for) = @_;  
      my @pos;  
      if ($search_for =~ /^\/.+\/$/){  
           $search_for =~ s/^\/|\/$//g;  
           @pos = grep { $_ =~ /$search_for/ } @{$array};  
      } else {  
           @pos = grep { $_ eq $search_for } @{$array};  
      }  
      (scalar @pos)?return (1,\@pos):return (0);  
 }  

2 de noviembre de 2010

Formación en Bioinformática en Europa


 
La semana pasada estuvimos en las X Jornadas de Bioinformática, que desde hace ya un par de ediciones aglutinan a la mayoría de grupos de España y Portugal que se dedican a esto. Además de un montón de presentaciones científicas y paneles de pósters,tuvimos una sesión dedicada a discutir opciones de formación en Bioinformática. De entre todo lo que se dijo me gustaría destacar un calendario, mantenido por Pedro Fernandes, del Gulbenkian Training Program in Bioinformatics, donde se irán publicando las eventos de formación que se organicen en Europa. El enlace al calendario es: http://sites.google.com/site/bioinftraining/calendarofevents

Un saludo,
Bruno

22 de octubre de 2010

Alineamiento de perfiles con HHalign

Hola,
hoy discutiré un problema habitual cuando trabajamos con familias de proteínas, el de alinear entre si dos grupos de secuencias previamente alineadas (MSAs=multiple sequence alignments) , como se ilustra en el esquema siguiente:


Dadas las secuencias del grupo 1, que hemos alineado entre si con ayuda de programas como muscle, clustal o MAFFT, y las secuencias del grupo 2, que hemos procesado de manera similar, ahora queremos ver cómo se alinean todas juntas, por ejemplo para inferir diferentes funciones o estructuras.

En el caso trivial sería tan sencillo como poner todas las secuencias en el mismo archivo y hacer un alineamiento múltiple de una vez. Sin embargo, el caso que discutimos aquí es más complicado, el que se da cuando no es sencillo alinear una secuencia del grupo 1 con otra del 2, porque han divergido bastante. Entonces podemos probar a convertir el alineamiento MSA1 a un perfil y a alinearlo contra el perfil correspondiente a MSA2, como si fuera un alineamiento pareado. Podemos incorporar, además de la información evolutiva capturada en cada perfil, información de estructura secundaria para ayudar a guiar el alineamiento.

Todo esto es sencillo de realizar con el paquete de programas HHsearch, que podemos descargar de ftp://toolkit.lmb.uni-muenchen.de/HHsearch.

Los pasos a seguir son 3:

1) Obtener archivos FASTA de MSA1.faa (30 secuencias) y MSA2.faa (35 secuencias). Para mayor precisión en el alineamiento es recomendable añadir en ambos archivos una predicción de estructura secundaria, en el formato de PSIPRED. Por ejemplo, el contenido de MSA1.faa comenzaría con las siguientes líneas:
>aa_pred
MVVSIGVFDGVHI--GHQKVL...
>ss_pred
CEEEEECCCCEEH--HHHHHH...
>ss_conf
8888617775007--889999...
>seq1 (y todas las 29 siguintes a continuación)
MVVSIGVFDGVHI--GHQKVL...
2) Convertir MSA1.faa y MSA2.faa en perfiles de Markov (MSA1.hhm y MSA2.hhm), con el ejecutable hhmake:

$ ~/soft/HHsearch/hhmake -i MSA1.faa -seq 30
$ ~/soft/HHsearch/hhmake -i MSA2.faa -seq 35


3) Alinear MSA1 contra MSA2, ya sea de manera global o local:

$ ~/soft/HHsearch/hhalign -i MSA1.hhm -t MSA2.hhm -glob -seq 65 -ofas globalMSA.faa

La salida obtenida es algo como:
REMARK: in -mac -global mode -mact is forced to 0
Query file is in HHM format
Read in HMM seq1 with 467 match states and effective number of sequences = 3.1
Query file is in HHM format
Read in HMM seq31 with 576 match states and effective number of sequences = 1.8
Using maximum accuracy (MAC) alignment algorithm ...
Printing alignments in FASTA format to globalMSA.faa 
Aligned seq1 with seq31: Score = 218.04   P-value = 1     
Done
Un ejemplo de aplicación real de esta técnica de alineamiento lo tenéis en el artículo http://www.biomedcentral.com/1471-2148/10/311 .

15 de octubre de 2010

Cómo crear un blog de laboratorio

Hoy en día todo laboratorio que se precie tiene que tener una presencia mínima en internet. Eso supone tener una página web donde se puedan consultar los datos de contacto, líneas de investigación, proyectos, publicaciones, ofertas de trabajo, etc. referentes al laboratorio. Toda esta información será una buena carta de presentación para otros investigadores que busquen colaboraciones, estudiantes y doctores que quieran entrar en el laboratorio, y toda aquella gente que esté interesada en los temas que se investigan en el laboratorio. Por ejemplo, ésta es nuestra web: http://www.eead.csic.es/compbio/

 La información de una página web de laboratorio es más o menos estática, esto es, no cambian sus contenidos frecuentemente, normalmente sólo se actualiza 2 o 3 veces al año. Sin embargo la actual tendencia de Web 2.0 hace que la gente busque en internet información cada vez más actualizada, dinámica y cambiante. Por ello un Blog de Laboratorio puede servirnos para ampliar la información de nuestra página web de una forma más flexible, personal y frecuente. A continuación se darán unas pautas y consejos para ayudar a cualquier investigador que quiera crear su propio blog o un blog para su laboratorio.

¿Por qué crear un Blog?
Un blog nos permitirá publicar de forma rápida y sencilla contenidos en internet. De esta forma mejorará la visibilidad de nuestras investigaciones en Google. Además, el vocabulario que se suele usar en los blogs es más directo y sencillo, con lo cual se amplía también el rango de gente a la que le pueda interesar lo que publiquemos.

¿Cómo y dónde crear un Blog?
Actualmente existen dos grandes plataformas gratuitas para crear un blog en internet: Blogger y Wordpress. Las dos son muy buenas opciones, y decidirse por una u otra a veces es cuestión de gustos, este artículo os puede ayudar a decidir: http://entretecnologia.blogspot.com/2010/04/blogger-y-wordpress-como-elegir.html En mi humilde opinión, Blogger es más sencillo y Wordpress es más completo y personalizable, así que para un laboratorio puede servirnos muy bien la primera opción.
 
¿Qué escribir en el Blog?
No tiene sentido crear un Blog para el laboratorio y escribir en él cosas de política, así que mi recomendación es que sea un blog con temática cercana a la de las investigaciones del laboratorio. No por ello han de ser contenidos técnicos y complicados, si no no los leería nadie. Y tampoco hay que ser muy estrictos en respetar siempre la temática del blog y de vez en cuando se puede uno permitirse el lujo de hablar de algo diferente, como por ejemplo este post. Un buen consejo es escribir párrafos cortos y poner de vez en cuando alguna imagen.

¿Cada cuánto debo escribir en el Blog?
No hay ninguna obligación de escribir todos los días en el blog, pero sí que sería conveniente publicar algo todas las semanas, o al menos un artículo todos los meses. Cuanto más cosas escribamos y más interesantes sean, más gente nos leerá y conseguiremos nuestro objetivo de informar y dar a conocer nuestro laboratorio.

¿Y si leo este artículo y creo un Blog?
Si es así publica un comentario en este post con el título, temática y dirección de tu Blog de Laboratorio, así empezarás a darlo a conocer...


Espero que algún investigador se anime después de leer esto, y empiecen los comentarios :D

Publica tu código


Hola,
hoy en vez de evaluar un programa o unas líneas de código os sugiero una columna
publicada esta misma semana en la influyente revista Nature, donde el autor, Nick Barnes, director de la Climate Code Foundation, defiende que los científicos salimos ganando y de hecho contribuimos a difundir nuestro propio trabajo, cuando publicamos el código fuente de los programas que escribimos para nuestros proyectos, aunque nunca lo hubiésemos escrito con la idea de distribuirlo. En caso contrario, defiende el autor, generalmente ese código se pierde y se convierte en abandonware (fuente de la imagen: http://top-buzz.com/2009/10).

Cualquiera que lleve cierto tiempo trabajando en Bioinformática seguro que tiene programitas y scripts que llevan tiempo languideciendo en directorios o incluso en copias de seguridad condenadas al olvido. Lo que propone el autor es que muchos de esos programas se podrían publicar por ejemplo como material suplementario junto con los artículos para los que se escribieron. Yo añadiría que son todavía más visibles, y por tanto potencialmente más útiles para otros usuarios, si los colgamos de la web en blogs como éste.


El enlace al artículo original en inglés es:
http://www.nature.com/news/2010/101013/full/467753a.html

8 de octubre de 2010

Biostar, consultas bioinformáticas a la comunidad

Hola, mi colega Pedro Olivares me habla de un foro en inglés donde puedes preguntar dudas sobre tareas bioinformáticas. Dice Pedro:

Aprovecho para recomendarte este sitio de internet de preguntas y respuestas en bioinformática. Tiene un sistema raro de reputación y hay gente bien clavada ahí dando consejos y luciéndose ante la comunidad.

http://biostar.stackexchange.com

Por lo que he visto el foro tiene mucha vida, así que igual es un recurso interesante para solucionar problemas urgentes, un saludo,
Bruno

6 de octubre de 2010

Anuncio del curso "Statistical Methods for mRNA-Seq and ChIP-seq using Bioconductor"

El departamento de bioinformática del Centro de investigación Príncipe Felipe (CIPF) de Valencia organiza este curso el 10 de Noviembre, cuesta 500 euros.

The short course will provide an overview of statistical methods and software for the analysis of next generation sequencing (NGS) data, with emphasis on transcriptome analysis (mRNA-Seq) and the study of DNA-protein interactions (ChIP-Seq). The morning lectures on statistical methodology will be followed by afternoon lab sessions demonstrating the application of the methodology to mRNA-Seq and ChIP-Seq data using R software packages from the Bioconductor Project (www.bioconductor.org).

The course will be taught by Bioconductor developeprs from the Division of Biostatistics, School of Public Health, University of California, Berkeley and has been organized in collaboration with the Department of Bioinformatics of the CIPF. The course is especially suited to bioinformaticists and advanced users of bioinformatic tools interested in dealing with NGS data. It is also strongly recommended to personnel in charge of data analysis in core facilities implementing NGS.

Maś información en http://bioinfo.cipf.es/RNA-seq2010

27 de septiembre de 2010

Ordenamiento de resultados de BLAST

Un problema con el que me he encontrado recientemente es el de ordenar resultados de BLAST contenidos en diferentes ficheros. Para definir mejor el problema supongamos que tenemos 4 genomas A,B,C y D y queremos comparar, todas contra todas, sus secuencias de proteínas por medio del programa blastp. Una manera de hacerlo sería hacer 4x4 comparaciones por parejas (AA,AB,AC,AD,BA,BB,BC,BD,CA,CB,CC,CD,DA,DB,DC,DD) teniendo en cuenta que la dirección en BLAST normalmente importa.

Una vez completados estos 16 trabajos de BLAST encontraremos que cada uno de los archivos de salida están internamente ordenados en términos de E-value. Sin embargo, dada la secuencia 1 del genoma A, si queremos averiguar en qué genomas se encuentra su secuencia más similar (best hit), deberemos mirar en todos los archivos (7 de 16) donde aparece el genoma A. Otra alternativa, la que da pie a este artículo, es fusionar, intercalar y ordenar (merge-sort) esos 16 archivos internamente ordenados para crear uno sólo, que facilite posteriores consultas.
Sirvan de muestra los archivos AA.blast:
16 16 100.00 512 0 0 1 512 1 512 0.0 1036
16 18 24.88 406 261 10 57 443 34 414 2e-24  114
16 78 25.26 475 303 12 1 452 15 460 1e-19 97.8
y AB.blast:
16 582 25.97 362 232 9 95 443 76 414 5e-23  108
16 637 28.00 300 193 5 86 377 91 375 3e-21  103

que contienen, en formato tabular, los alineamientos significativos de la secuencia 16 dentro del genoma A (contra las secuencias 16,18,78) y dentro del genoma B (582,637), con valores esperados de 0.0, 2e-24,1e-19,5e-23,3e-21 respectivamente.
El problema que vamos a resolver es como combinar estos dos archivos (o en general N archivos) en uno sólo hasta obtener una lista ordenada por E-value:
16 16 100.00 512 0 0 1 512 1 512 0.0 1036
16 18 24.88 406 261 10 57 443 34 414 2e-24  114
16 582 25.97 362 232 9 95 443 76 414 5e-23  108
16 637 28.00 300 193 5 86 377 91 375 3e-21  103
16 78 25.26 475 303 12 1 452 15 460 1e-19 97.8

Sin más preámbulos, el siguiente código Perl implementa dos posibles soluciones:

 package BlastSort;   
   
 use strict;  
 use Symbol qw(gensym);  
 use sort 'stable';   
   
 require Exporter;  
 use vars qw(@ISA @EXPORT);  
 @ISA = 'Exporter';  
 @EXPORT = qw(merge_BLAST_files merge_BLAST_files_GNUsort);  
   
 sub merge_BLAST_files   
 {  
    # Adapted from File-Sort-1.01(http://search.cpan.org/perldoc?File::Sort)  
    # Assumes infiles are BLAST output files in tabular format with sequences  
    # identified by natural numbers, such as 11,12,1439,1440 in the sample:  
    #11   1439   78.24   625   136   0   4   628   6   630   0.0    993  
    #12   1440   80.88   272   52   0   1   272   1   272   7e-125    446  
    # Order of infiles is RELEVANT as merging is stable, so sequences from  
    # the first files will be given sorting priority  
      
    my ($outfile,@infiles) = @_;  
      
    my (%fh,%order,@fhorder,$filein,$id,$first,$n_of_fhs,$curr,$line);  
   
    # for Schwartzian transform (ST) see  
    # http://www.hidemail.de/blog/perl_tutor.shtml#sort_orcish_schwartzian  
    sub blastsort { $a->[2] <=> $b->[2] || $a->[3] <=> $b->[3] }  
    sub blastmap   
    {   
       my @tmp = split(/\s+/,$fh{$_});   
       [$_,$fh{$_},$tmp[0],$tmp[10]]   
       # returns anonymous array[4]: filehandle, line,query_id,E-value  
    }  
   
    ## open all input BLAST files and keep filehandles, ordered 0..N  
    $n_of_fhs = 0;  
   foreach $filein (@infiles)   
    {  
      $id = gensym(); # get valid id for filehandle  
         
       open($id,$filein) ||   
          die "# merge_BLAST_files : cannot read $filein: $!";  
   
       $order{$n_of_fhs} = $id;  
       $n_of_fhs++;  
    }  
    @fhorder = (0 .. $n_of_fhs-1);  
      
     
    ## open outfile and ensure IO buffer is used  
    $| = 0;   
   unlink($outfile) if(-s $outfile);  
    open(OUT,">$outfile") ||   
       die "# merge_BLAST_files : cannot create $outfile: $!";  
   
   ## get first BLAST line from all filehandles    
   %fh = map {  
     my $fh = $order{$_};  
     ($_ => scalar <$fh>);  
   } @fhorder;  
   
   ## start merging BLAST lines   
    while(scalar(@fhorder)>1)   
    {        
       ($first) = (map {$_->[0]} sort blastsort map &blastmap, @fhorder); #ST  
         
       print OUT $fh{$first};  
   
     $curr = $order{$first};  
     $line = scalar <$curr>;  
     if(defined($line)) # update current filehandle  
       {   
          $fh{$first} = $line;   
       }  
       else # exhausted filehandle  
       {   
          @fhorder = grep { $_ ne $first } @fhorder;  
       }  
   }  
      
    ## take care of last filehandle left and close file  
    print OUT $fh{$fhorder[0]};  
    $curr = $order{$fhorder[0]};  
    while(<$curr>){ print OUT }  
   close(OUT);    
      
    ## close all input files  
    foreach $id (0 .. $n_of_fhs-1){ close($order{$id}); }  
      
    if(!-s $outfile){ return 0 }  
    else{ return 1 }  
 }  
   
   
 sub merge_BLAST_files_GNUsort  
 {  
    my ($outfile,$tmpdirpath,$maxbufferMb,@infiles) = @_;  
   
    # local sort -k11g fails with: 0.00,0.006 and 1e-23 (coreutils v6.10)  
    # probably as LC_ALL is not set to 'POSIX'  
    # http://www.gnu.org/software/coreutils/faq/coreutils-faq.html   
   # #Sort-does-not-sort-in-normal-order_0021  
    $ENV{'LC_ALL'} = 'POSIX';  
     
   unlink($outfile) if(-s $outfile);  
   
   my $sort_command = "sort --temporary-directory=$tmpdirpath " .  
     "--buffer-size=$maxbufferMb -s -k1g -k11g -m " .  
     join(' ',@infiles)." > $outfile ";  
   
   system("$sort_command");    
   
    if(!-s $outfile){ return 0 }  
    else{ return 1 }  
 }  
   
 __END__  

El módulo contiene dos subrutinas,  merge_BLAST_files y merge_BLAST_files_GNUsort; mientras la primera muestra explícitamente cómo se desarrolla el ordenamiento externo en disco, manteniendo los N ficheros de entrada abiertos y sin guardar nada en memoria, la segunda subrutina es realmente una llamada a la utilidad sort del shell (GNU coreutils 6.10 en mi sistema), donde sí usamos la memoria para acelerar el ordenamiento, 500Mb en el siguiente ejemplo:

 #!/usr/bin/perl  
 use strict;  
 use Benchmark;  
 use BlastSort;  
   
 my @infiles = ( 'A-A.blast','A-B.blast' );    
   
 my ($outfile,$gnu_outfile,$start_time,$end_time,$sortOK) =   
    ('out.merge-sort.txt','out.gnu-sort.txt');  
   
 print "# number of BLAST files to merge-sort = ".scalar(@infiles)."\n";  
   
 $start_time = new Benchmark();  
 $sortOK = merge_BLAST_files($outfile,@infiles);  
 $end_time = new Benchmark();  
 print "\n# runtime (BlastSort): ".  
    timestr(timediff($end_time,$start_time),'all')."\n";  
   
 $start_time = new Benchmark();  
 $sortOK = merge_BLAST_files_GNUsort($gnu_outfile,'./',500,@infiles);  
 $end_time = new Benchmark();  
 print "\n# runtime    (GNU): ".  
    timestr(timediff($end_time,$start_time),'all')."\n";  

Como se observa en la gráfica, aunque la versión Perl es autoexplicativa, es exponencialmente más lenta que GNUsort. Por tanto, en general convendrá usar la segunda opción y la primera sólo tendrá interés si queremos reservar la RAM, a costa de tiempo de ejecución. La única posible complicación del GNUsort es que depende de la variable de ambiente LC_ALL, a la que deberemos dar el valor 'POSIX' para tener resultados correctos. De no ser así no ordena bien los E-valores.


17 de septiembre de 2010

¿Cuál es el mínimo valor del E-value de Blast?

Esta extraña pregunta ha surgido hoy en mi laboratorio, y analizando varios ficheros de resultados de Blast he llegado a la conclusión que el menor E-value de Blast es 1e-180, por debajo de este valor Blast devuelve 0.0 como resultado.


¿What is the minimum Blast E-value?
This strange question has an easy solution: 1e-180, it is the minimum Blast E-value, lower values are scored by Blast as 0.0.

13 de septiembre de 2010

Jornadas Bioinformáticas JBI 2010 (X Edición), nuestro laboratorio estará allí...

Las Jornadas Bioinformáticas son la cita anual obligada para los bioinformáticos españoles. Este año se celebrará su décima edición del 27 al 29 de Octubre en Torremolinos (Málaga). La organización de las mismas corre a cargo de la Universidad de Málaga, el Instituto Nacional de Bioinformática y la Red Portuguesa de Bioinformática. Este año el tema central es "La bioinformática aplicada a la medicina personalizada", sobre el cual se discutirá la integración de los campos de la biología, medicina e informática para el desarrollo de terapias más específicas y efectivas. Sin embargo, éste no será el único tema a tratar, también se compartirán resultados y experiencias en otros campos:
- Análisis de datos en técnicas de alto rendimiento como la secuenciación de nueva generación.
- Bioinformática estructural
- Algoritmos de biología computacional y técnicas de computación de alto rendimiento
- Análisis de secuencias, filogenética y evolución
- Bases de datos, herramientas y tecnologías de biología computacional
- Bioinformática en transcriptómica y proteómica
- Biología sintética y de sistemas

IN ENGLISH:


The Xth Spanish Symposium on Bioinformatics (JBI2010) will take place in October 27-29, 2010 in Torremolinos-Málaga, Spain. Co-organised by the National Institute of Bioinformatics-Spain and the Portuguese Bioinformatics Network and hosted by the University of Malaga (Spain).


This year, the reference topic is “Bioinformatics for personalized medicine” for which the conference will provide the opportunity to discuss the state of the art for the integration of the fields of biology, medicine and informatics. We invite you to submit your work and share your experiences in the following topics of interest including, but not limited to:
- Analysis of high throughput data    (NGS)
- Structural Bioinformatics
 - Algorithms for computational biology and HPC
 - Sequence analysis, phylogenetics and evolution
 - Databases, Tools and technologies for computational biology
-  Bioinformatics in  Transcriptomics and Proteomics
- System and Synthetic Biology



Nuestras aportaciones

Nuestro laboratorio va a participar en las Jornadas Bioinformáticas con tres contribuciones que presentaré a continuación:


3D-footprint: a database for the structural analysis of protein–DNA complexes
3D-footprint is a living database, updated and curated on a weekly basis, which provides estimates of binding specificity for all protein–DNA complexes available at the Protein Data Bank. The web interface allows the user to: (i) browse DNA-binding proteins by keyword; (ii) find proteins that recognize a similar DNA motif and (iii) BLAST similar DNA-binding proteins, highlighting interface residues in the resulting alignments. Each complex in the database is dissected to draw interface graphs and footprint logos, and two complementary algorithms are employed to characterize binding specificity. Moreover, oligonucleotide sequences extracted from literature abstracts are reported in order to show the range of variant sites bound by each protein and other related proteins. Benchmark experiments, including comparisons with expert-curated databases RegulonDB and TRANSFAC, support the quality of structure-based estimates of specificity. The relevant content of the database is available for download as flat files and it is also possible to use the 3D-footprint pipeline to analyze protein coordinates input by the user. 3D-footprint is available at http://floresta.eead.csic.es/3dfootprint with demo buttons and a comprehensive tutorial that illustrates the main uses of this resource.


The relation between amino-acid substitutions in the interface of transcription factors and their recognized DNA motifs

Transcription Factors (TFs) play a key role in gene regulation by binding to DNA target sequences. While there is a vast literature describing computational methods to define patterns and match DNA regulatory motifs within genomic sequences, the prediction of DNA binding motifs (DBMs) that might be recognized by a particular TF is a relatively unexplored field. Numerous DNA-binding proteins are annotated as TFs in databases; however, for many of these orphan TFs the corresponding DBMs remain uncharacterized. Standard annotation practice transfer DBMs of well known TFs to those orphan protein sequences which can be confidently aligned to them, usually by means of local alignment tools such as BLAST, but these predictions are known to be error-prone. With the aim of improving these predictions, we test whether the knowledge of protein-DNA interface architectures and existing TF-DNA binding experimental data can be used to generate family-wise interface substitution matrices (ISUMs). An experiment with 85 Drosophila melanogaster homeobox proteins demonstrate that ISUMs: i) capture information about the correlation between the substitution of a TF interface residue and the conservation of the DBM; ii) are valuable to evaluate TFs alignments and iii) are better classifiers than generic amino-acid substitution matrices and that BLAST E-value when deciding whether two aligned homeobox proteins bind to the same DNA motif.


101DNA: a set of tools for Protein-DNA interface analysis

Analysis of protein-DNA interfaces has shown a great structural dependency. Despite the observation that related proteins tend to use the same pattern of amino acid and base contacting positions, no simple recognition code has been found. While protein contacts with the sugar-phosphate backbone of DNA provide stability and yield very little specificity information, contacts between amino acid side-chains and DNA bases (direct readout) apparently define specificity, in addition to some constrains defined by DNA sequence-dependent features, namely indirect readout.
Recent approaches have proposed bipartite graphs as an structural way of analysing interfaces from a protein-DNA-centric viewpoint. With this perspective in mind, we have developed a set of tools for the dissection and comparison of protein-DNA interfaces. Taking a protein-DNA complex file in PDB format as input, the software generates a 2D matrix that represents a bipartite graph of residue contacts  obtained after applying a simple distance threshold that captures all non-covalent interactions. The generated 2D matrices allow a fast and simple visual inspection of the interface and have been successfully produced for the current non-redundant set of protein-DNA complexes in the 3D-footprint database.
As a second utility to compare 2 interfaces, the 101DNA software includes an aligment tool where a dynamic programming matrix is created with the Local Affine Gap algorithm and traced back as a finite state automata. The scores between pairs of  interface amino acid residues are calculated as a function of the observed contacts with DNA nitrogen bases. This tool produces local interface alignments which are independent of the underlying protein sequence, but that faithfully represent the binding architecture. Preliminary tests show that these local alignments successfully identify binding interfaces that share striking similarity despite belonging to different protein superfamilies, and these observations support this graph-theory approach.

10 de septiembre de 2010

Cómo saber si un valor es numérico mediante una simple expresión regular

Hoy no estoy muy inspirado para escribir en el blog, así que he repasado mis librerías de subrutinas de perl para ver si encontraba algo curioso para publicar y me he encontrado con esto...
¿Quién no ha tenido nunca el problema de comprobar si una expresión es numérica o no? 
Perl no posee una función que lo haga automáticamente (por lo menos que yo conozca), sin embargo con una simple línea de código podemos salir de dudas:

 # Return TRUE if an expression is numerical  
 sub is_numeric {  
      my ($exp) = @_;  
      if ($exp =~ /^(-?[\d.\-]*e[\d.\-\+]+|-?[\d.\-]\^[\d.\-]+|-?[\d.\-]+)$/){  
      # Corregida siguiendo las indicaciones de Joaquín Ferrero:
      if ($exp =~ /^-?(?:[\d.-]+*e[\d.+-]+|\d[\d.-]*\^[\d.-]+|[\d.-]+)$/){
           return 1;  
      } else {  
           return 0;  
      }  
 }  

1 de septiembre de 2010

Adaptando scripts antiguos para Blast+

Hola,
los que uséis programas de la familia BLAST habitualmente habréis notado que las últimas versiones instalables, que se pueden descargar por FTP, son de la rama nueva de desarrollo Blast+, que tiene algunas nuevas funcionalidades muy interesantes, pero que cambia los nombres de los ejecutables y la forma de pasar argumentos que sabíamos usar.


Sin embargo, los desarrolladores del NCBI ya habían previsto esta transición y acompañan a los nuevos binarios un script Perl, que se llama legacy_blast.pl, que puede ayudarnos a reconvertir código que invoca versiones antiguas de BLAST. Con este programa podemos por ejemplo traducir este comando de la versión 2.2.18, que tiene opciones de filtrado un tanto especiales:

$ blastall -F 'm S' -z 10000 -p blastp -i problema.faa -d bases_datos/seqs.faa

a esta llamada al binario de la versión 2.2.24+:

$ blastp -db bases_datos/seqs.faa -query problema.faa -dbsize 10000 -seg yes -soft_masking true

Por cierto, en mis pruebas veo que el nuevo binario puede aprovechar una base de secuencias formateada con el formatdb antiguo, que ahora se llama makeblastdb.

Un saludo,
Bruno

24 de agosto de 2010

Código en C dentro de un programa Perl

Hola, siguiendo el hilo de los comentarios a la entrada anterior de Álvaro, hoy pongo un ejemplo de como incluir una o más subrutinas escritas en C dentro de un texto en Perl, recurriendo al módulo Inline:
 #!/usr/bin/perl -w 
 # programa 7 
 use Inline C; 
 use strict; 
 my $BASE = 0.987654321; 
 my $EXP = 1000; 
 print "# $BASE^$EXP = \n"; 
 print "# ".potencia_perl($BASE,$EXP)." (perl)\n"; 
 print "# ".potencia_C($BASE,$EXP)." (C)\n"; 
 sub potencia_perl 
 { 
    my ($base,$exp) = @_; 
    my $result = $base; 
    for(my $p=1; $p<$exp; $p++) { $result *= $base; } 
    return $result; 
 } 
 __END__  
 ### alternativa en C 
 __C__ 
 #include <stdio.h> 
 float potencia_C( double base, double exp )  
 { 
    double result = base; 
    int p; 
    for(p=1; p<exp; p++) 
    {  
       result *= base; 
    }    
    return result; 
 } 

La salida obtenida pone en evidencia que diferentes lenguajes suelen tener diferentes precisiones a la hora de representar números y operar con ellos:

$ perl prog7.pl 
# 0.987654321^1000 = 
#  4.02687472213071e-06 (perl)
#  4.02687464884366e-06 (C)

19 de agosto de 2010

Generar todas las posibles combinaciones posibles de n nucleótidos o aminoácidos

Erróneamente decimos que queremos generar todas las posibles 'combinaciones' de n letras, nucleótidos o aminoácidos cuando el término correcto sería 'variaciones con repetición' (Permutations-Variations-Combinations). Por ejemplo, las 16 posibles variaciones con repetición de 2 nucleótidos serían: AA, AC, AG, AT, CA, CC, CG, CT, GA, GC, GG, GT, TA, TC, TG, TT. El número de variaciones con repetición de n nucleótidos es 4^n y de n aminoácidos 20^n.

El siguiente código, es una adaptación de un código en PHP para generar todas las variaciones con repetición de n elementos, en nuestro caso, 2 aminoácidos y 4 nucleótidos.

 #!/usr/bin/perl -w  
 # Generate all the variations with repetition of n letters  
 my @aminoacids = ('A','R','N','D','C','Q','E','G','H','L','I','K','M','F','P','S','T','W','Y','V');  
 my @nucleotides = ('A','C','G','T');  
 # Generate all the variations of 2 aminoacids  
 my $aas = variations(\@aminoacids,2);  
 # Generate all the variations of 4 nucleotides  
 my $nts = variations(\@nucleotides,4);  
 # Print the results  
 print "\nAll the variations with repetition of 2 aminoacids: ";  
 foreach my $vari (@{$aas}){  
      foreach my $elem (@{$vari}){  
           print "$elem";  
      }  
      print " ";  
 }  
 print "\n";  
 print "\nAll the variations with repetition of 4 nucleotides: ";  
 foreach my $vari (@{$nts}){  
      foreach my $elem (@{$vari}){  
           print "$elem";  
      }  
      print " ";  
 }  
 print "\n\n";  
 sub variations {  
      my ($letters,$num) = @_;  
      my $last = [map $letters->[0] , 0 .. $num-1];  
      my $result;  
      while (join('',@{$last}) ne $letters->[$#{$letters}] x $num) {  
           push(@{$result},[@{$last}]);  
           $last = char_add($letters,$last,$num-1);  
           print '';  
      }  
      push(@{$result}, $last);  
      return $result;  
 }  
 sub char_add{  
      my ($digits,$string,$char) = @_;  
      if ($string->[$char] ne $digits->[$#{$digits}]){  
           my ($match) = grep { $digits->[$_] eq $string->[$char]} 0 .. $#{$digits};  
           $string->[$char] = $digits->[$match+1];  
           return $string;  
      } else {  
           $string = changeall($string,$digits->[0],$char);  
           return char_add($digits,$string,$char-1);  
      }  
 }  
 sub changeall {  
      my ($string,$char,$start,$end) = @_;  
      if (!defined($start)){$start=0;}  
      if (!defined($end)){$end=0;}  
      if ($end == 0) {$end = $#{$string};}  
      for(my $i=$start; $i<=$end; $i++){  
           $string->[$i] = $char;  
      }  
      return $string;  
 }  

Para terminar una breve lección de estadística, fuente: Aulafacil.com

a) Combinaciones:
Determina el número de subgrupos de 1, 2, 3, etc. elementos que se pueden formar con los "n" elementos de una nuestra. Cada subgrupo se diferencia del resto en los elementos que lo componen, sin que influya el orden.
Por ejemplo, calcular las posibles combinaciones de 2 elementos que se pueden formar con los números 1, 2 y 3.
Se pueden establecer 3 parejas diferentes: (1,2), (1,3) y (2,3). En el cálculo de combinaciones las parejas (1,2) y (2,1) se consideran idénticas, por lo que sólo se cuentan una vez.
b) Variaciones:
Calcula el número de subgrupos de 1, 2, 3, etc.elementos que se pueden establecer con los "n" elementos de una muestra. Cada subgrupo se diferencia del resto en los elementos que lo componen o en el orden de dichos elementos (es lo que le diferencia de las combinaciones).
Por ejemplo, calcular las posibles variaciones de 2 elementos que se pueden establecer con los número 1, 2 y 3.
Ahora tendríamos 6 posibles parejas: (1,2), (1,3), (2,1), (2,3), (3,1) y (3,3). En este caso los subgrupos (1,2) y (2,1) se consideran distintos.
c) Permutaciones:
Cálcula las posibles agrupaciones que se pueden establecer con todos los elementos de un grupo, por lo tanto, lo que diferencia a cada subgrupo del resto es el orden de los elementos.
Por ejemplo, calcular las posibles formas en que se pueden ordenar los número 1, 2 y 3.
Hay 6 posibles agrupaciones: (1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2) y (3, 2, 1)

6 de agosto de 2010

Script para conocer el formato de un archivo

En el campo de la bioinformática se trabaja con numerosos formatos de archivos de datos biológicos (FASTA, GenBank, PDB...). Habitualmente tenemos que conocer el formato de los datos de un archivo que usamos como entrada en un script. Sin embargo, esto no es tarea fácil puesto que diferentes formatos pueden contener el mismo tipo de datos de forma poco estandarizada. En estos casos tenemos que usar nuestros conocimientos previos sobre el formato para establecer el tipo de archivo. El siguiente código de perl permite introducir información sobre los formatos más habituales en forma de expresiones regulares para posteriormente leer un fichero y asignar una puntuación a cada uno de los formatos definidos. De acuerdo a la puntuación final conoceremos el formato de fichero más probable.


 # Obtains the format of a file  
 sub obtain_file_format {  
      my ($file) = @_;  
      open(FILE,$file)|| die "# $0 : cannot read $file\n";  
      my %score;  
      while(<FILE>){  
           if (/^#/){  
                next;  
           }  
           if (/^>/){  
                $score{'fasta'}++;  
           }  
           if (/^(LOCUS|DEFINITION|ACCESSION|FEATURES)/){  
                $score{'genbank'}++;  
           }  
           if (/^(HEADER|TITLE|COMPND|SEQRES|ATOM)/){  
                $score{'pdb'}++;  
           }  
           if (/^(DE|VV|XX)/){  
                $score{'transfac'}++;  
           }  
      }  
      close(FILE);  
      my @sorted_scores;  
      if (@sorted_scores = sort { $score{$b} <=> $score{$a} } (keys %score)) {  
           return $sorted_scores[0];  
      } else {  
           return 'unknown';  
      }  
 }  

30 de julio de 2010

Bioinformática en la nube a partir de 10000 usuarios

La comunidad bioinformática tiene ya una cierta tradición en poner a disposición de los usuarios el software desarrollado a través de la red. Esto sirve para dar a conocer las últimas ideas y algoritmos que se van desarrollando en los laboratorios, de manera que se ponen a prueba en escenarios reales en manos de usuarios que no son sus creadores. Por supuesto, también sirve para ir ampliando el repertorio de herramientas que cualquiera con conexión a Internet puede usar, normalmente de forma gratuita.

Para el usuario final éstas son a menudo simplemente aplicaciones accesibles desde un navegador web, pero para sus creadores son programas más o menos complejos que viven en servidores web (como máquinas virtuales) que consumen luz y hay que administrar, que son a veces objeto de ataques informáticos, y que dependiendo del uso que tengan exigirán ampliaciones de hardware.


Qué coste real tienen estos servidores? El coste absoluto de mantener un servidor es difícil de promediar, porque dependerá de costes (luz, hardware, la carga de usuarios) y beneficios (evaluación peer to peer de los algoritmos, corrección de errores, retroalimentación de los usuarios, incluso citas bibliográficas) variables.
Otra cuestión  es si es mejor alojar estos servidores en nuestras propias máquinas o aprovecharnos de las arquitecturas de computación en la nube, que son escalables. Para ayudarnos a tomar esta decisión un reciente artículo en Bioinformatics del grupo de Cedric Notredame, donde prueban la versión paralelizada de su famoso alineador múltiple Tcoffee (basado en el principio de consistencia), sugiere que con los costes actuales sólo compensa alojar este tipo de servicios web en arquitecturas en nube a partir de 10000 usuarios al mes.

Nos vemos a la vuelta de las vacaciones.

19 de julio de 2010

Alineamiento perfil-perfil teniendo en cuenta la estructura secundaria

Con el alineamiento múltiple de secuencias de proteínas relacionadas se puede construir un perfil. A parte de su uso en técnicas de fold recognition o para predicción de estructura secundaria, este se puede utilizar para compararlo con las secuencias de otras proteínas e incluso con otro grupo de proteínas relacionadas entre sí, en este caso mediante una comparación perfil-perfil.

En el siguiente ejemplo, escrito en python y del que hay un ejemplo en perl aquí, vamos a realizar la comparación entre 2 perfiles, y apoyaremos el alineamiento con la información disponible de la estructura secundaria predicha para ambos. Utilizaremos un algoritmo de programación dinámica, por lo que si no se está familiarizado con la técnica se puede acudir a uno de los muchos libros sobre algoritmos o al excelente artículo de Sean R Eddy What is dynamic programming?, el cual se basa precisamente en ejemplos de alineamiento de secuencias biológicas.

Básicamente, necesitaremos un archivo PSSM y otro PSIPRED, perfil y estructura secundaria, respectivamente, para cada uno de los 2 grupos de proteínas. Puedes encontrar archivos de ejemplo en este curso, donde ya enlazamos para el ejemplo en perl.

Remarcar, entre los aspectos más importantes del código, que la puntuación entre perfiles se calcula como la media entre los valores de las 2 matrices y que se asigna una puntuación extra al caso de que 2 posiciones alineadas tengan el mismo estado de estructura secundaria.

 #!/usr/bin/python   

from __future__ import with_statement
import re

__title__ = "Secondary structure and pairwise DP profile-profile alignment"
__author__ = "Carlos P Cantalapiedra"
__institution__ = "CSIC AulaDei Laboratory of Computational Biology"

### MAIN STARTS UNDERNEATH THE FUNCTIONS ###

# Function that parses PSSM files
# The format of the array is:
# pssm[
# {'aa':aa1, 'weights', [w1, w2, ...]},
# {'aa':aa2, 'weights', [w1, w2, ...]},
# .
# .
# .
# ]
def read_PSSM(pssm_file):
pssm = [] # Returning dictionary
aa = ''
weights = []

# Finds the aminoacid (\w) and his weights ((\s*-?\d*.\d*)+)
reg_exp = re.compile(r"^\s*\d+\s+(\w)((\s*-?\d*.\d*)+)")

# NOTE: with the with_statement there's no need to close the file
with open(pssm_file, 'r') as file:
for line in file:
result = reg_exp.search(line)

if result:
try:
aa = result.group(1)
weights = str(result.group(2)).split()

except IndexError:
raise IOError(\
"There is a problem with the format of the PSSM file "+str(pssm_file))

pssm.append({'aa' : aa, 'weights' : weights})
return pssm

# Function that parses "Pred:" line of PSIPRED files,
# to obtain the predicted secondary structures
def read_PSIPRED(psipred_file):
psipred = []

reg_exp = re.compile(r"^\s*Pred:\s*(\w+)")

# NOTE: with the with_statement there's no need to close the file
with open(psipred_file, 'r') as file:
for line in file:
result = reg_exp.search(line)

if result:
try:
[psipred.append(a) for a in result.group(1)]

except IndexError:
raise IOError(\
"There is a problem with the format of the PSIPRED file "+str(psipred_file))

return psipred

# Function that returns the bonus score associated to secondary structure
def get_secondary_score(st1, st2):
global secondary_score
score = 0

if st1 == st2 and st2 != 'C':
score = secondary_score

return score

# Function that creates de Dynamic Programming matrix
# The format of the matrix is going to be:
# matrix[
# row(0) [('r', gap), ('r', gap), ..., ('r', gap)] # len(row) = len(prot2) + 1
# row(1) [('c', gap), (path, score), ..., (path,score)]
# row(2) [('c', gap), (path, score), ..., (path,score)]
# .
# .
# .
# row(len(prot1)) [('c', gap), (path, score), ..., (path, score)]
# # len(col) = len(prot1) + 1
# ]
#
def create_matrix_DP(pssm1, pssm2, psipred1, psipred2):
global aa_order
global gapcost

matrix_DP = [] # Returning matrix
scores = []
score = 0
score1 = 0 # score1 and score2 serve several purposes
score2 = 0
aa1_count = 0 # 2 indexes
aa2_count = 0
path = '' # 'd': diagonal, 'c': column, 'r': row

try:

# Fisrt row is set up with gap values
init_row = [('r', i * gapcost) for i in range(0, len(psipred2) + 1)]
matrix_DP.append(init_row)

for aa1 in pssm1:
aa1_count+=1

index1 = aa_order.index(aa1['aa']) # aa1 position in PSSM columns

scores = []
aa2_count = 0

for aa2 in pssm2:
aa2_count+=1

# First column is set up with gap values
if aa2_count == 1:
scores.append(('c', aa1_count * gapcost))

index2 = aa_order.index(aa2['aa'])

score1 = int(aa2['weights'][index1])
score2 = int(aa1['weights'][index2])

score = (score1 + score2) / 2.0
score += get_secondary_score(psipred1[aa1_count-1], psipred2[aa2_count-1])

# Now does apply the function for matrix positions:
# Sij = max (score(-1,-1)+score, score(-1,0)+gap, score(0,-1)+gap)
# [1] is the score in the tuple (path, score)
score += matrix_DP[aa1_count - 1][aa2_count - 1][1]
score_1 = matrix_DP[aa1_count - 1][aa2_count][1] + gapcost
score_2 = scores[aa2_count - 1][1] + gapcost

if score_1 >= score and score_1 >= score_2:
path = 'c' # Column
score = score_1

elif score_2 >= score:
path = 'r' # Row
score = score_2

else:
path = 'd' # Diagonal

# Adds the score to the list of current row
scores.append((path, score))

# Adds a row of scores
matrix_DP.append(scores)

except ValueError, value:
raise Exception("Error. Maybe there is a rare aminoacid in the PSSM.\n"+value)
except IndexError, value:
raise Exception("Error. "+str(value))

return matrix_DP

# Recursive function that recovers the alignment
def visit_alignment(i, j, output):
global matrix_DP # This could be changed to be a parameter by reference
global prot1_pssm, prot2_pssm, prot1_psipred, prot2_psipred

if i == 0 and j == 0: return # END OF RECURSIVE SEARCH

join = ' '
if prot1_pssm[i - 1]['aa'] == prot2_pssm[j - 1]['aa']: join = '|'

path = matrix_DP[i][j][0]
score = matrix_DP[i][j][1]

if path == 'd':
visit_alignment(i - 1, j - 1, output)
output.append((prot1_psipred[i - 1], prot1_pssm[i - 1]['aa'],\
join, prot2_pssm[j - 1]['aa'], prot2_psipred[j - 1]))

elif path == 'c':
visit_alignment(i - 1, j, output)
output.append((prot1_psipred[i - 1], prot1_pssm[i - 1]['aa'], join, '-', '-'))

elif path == 'r':
visit_alignment(i, j - 1, output)
output.append(('-', '-', join, prot2_pssm[j - 1]['aa'], prot2_psipred[j - 1]))

else:
raise Exception("Invalid path "+path)

return


############################ MAIN ##################################

# Global variables: program parameters
aa_order = ('ARNDCQEGHILKMFPSTWYV'); # This must copy aa order in PSSM columns
gapcost = -8
secondary_score = 2

# File locations: please, change this variables to suit your needs
prot1_pssm_file = '../pssm/1ngk_A.pssm'
prot2_pssm_file = '../pssm/1s69_A.pssm'
prot1_psipred_file = '../psipred/1ngk_A.psipred'
prot2_psipred_file = '../psipred/1s69_A.psipred'

# Printing title and author
print
print "*** "+__title__+" ***"
print "Author: "+__author__+" at "+__institution__

# Printing parameters
print
print "Aminoacids and order (PSSM files): "+aa_order
print "Gap penalty = "+str(gapcost)
print "Equal secondary structure score = "+str(secondary_score)

### 1) Opening and parsing PSSM and PSIPRED files

prot1_pssm = []
prot2_pssm = []
prot1_psipred = []
prot2_psipred = []

try:
print
prot1_pssm = read_PSSM(prot1_pssm_file)
prot2_pssm = read_PSSM(prot2_pssm_file)
print "> PSSM files successfully opened."

prot1_psipred = read_PSIPRED(prot1_psipred_file)
prot2_psipred = read_PSIPRED(prot2_psipred_file)
print "> PSIPRED files successfully opened."

except IOError, value:
print "There has been an error reading input files:"
print value

### 2) Creating the matrix of DP

matrix_DP = [[]]
len_prot1 = len(prot1_psipred)
len_prot2 = len(prot2_psipred)

try:
matrix_DP = create_matrix_DP(prot1_pssm, prot2_pssm, prot1_psipred, prot2_psipred)

except Exception, value:
print value

### 3) Printing the maximum score

print
print "Maximum score "+str(matrix_DP[len_prot1][len_prot2][1])

### 4) Recovering and printing the alignment

output = [] # [(ss1, aa1, align, aa2, ss2)]
# Lines to print
ss1 = ''
aa1 = ''
align = ''
aa2 = ''
ss2 = ''

try:
# Recursively recovering the alignment
visit_alignment(len_prot1, len_prot2, output)

# Printing the alignment
print
for nts in output:
ss1 += str(nts[0])
aa1 += str(nts[1])
align += str(nts[2])
aa2 += str(nts[3])
ss2 += str(nts[4])
print ss1
print aa1
print align
print aa2
print ss2

except Exception, value:
print "Error while tracing the alignment."
print value
else: # END OF PROGRAM
print
print "*** Alignment finished successfully ***"


Esta sería la salida del programa, donde es importante fijarse en el alineamiento obtenido (que se muestra aquí en 2 fragmentos por motivo de espacio) y cómo, a pesar de que pocos aminoácidos coinciden, los estados de estructura secundaria se muestran alineados, lo que podría interpretarse como que la estructura está más conservada que la secuencia.

 *** Secondary structure and pairwise DP profile-profile alignment ***  
Author: Carlos P Cantalapiedra at CSIC AulaDei Laboratory of Computational Biology

Aminoacids and order (PSSM files): ARNDCQEGHILKMFPSTWYV
Gap penalty = -8
Equal secondary structure score = 2

> PSSM files successfully opened.
> PSIPRED files successfully opened.

Maximum score 444.0

CCHHHHHCCHHHHHHHHHHHHHHHHCCHHHHHHCCCCCHHHHHHHHHHHHHHHHCCCCCCCCCCCCCCH
KSFYDAVGGAKTFDAIVSRFYAQVAEDEVLRRVYPEDDLAGAEERLRMFLEQYWGGPRTYSEQRGHPRL
| || | | || | | | | || || |
STLYEKLGGTTAVDLAVDKFYERVLQDDRIKHFFADVDMAKQRAHQKAFLTYAFGGTDKYDGRYMREAH
CCHHHHHCCHHHHHHHHHHHHHHHHCCHHHHHHCCCCCHHHHHHHHHHHHHHHHCCCCCCCCCCHHHHH

HHCCCCCCCCHHHHHHHHHHHHHHHHHHCCCCCCHHHHHHHHHHHHHHHH--HHHCCCC
RMRHAPFRISLIERDAWLRCMHTAVASIDSETLDDEHRRELLDYLEMAAH--SLVNSPF
|| | || |
KELVENHGLNGEHFDAVAEDLLATLKEMG---VPEDLIAEVAAVAGAPAHKRDVLNQ--
HCCCCCCCCCHHHHHHHHHHHHHHHHHCC---CCHHHHHHHHHHHHHHHHHHHHHCC--

*** Alignment finished successfully ***

13 de julio de 2010

Eliminar secuencias redundantes y transcritos repetidos (II, con CD-HIT)

Como complemento a la entrada anterior de Álvaro os paso un enlace al programa
CD-HIT, que es una opción muy cómoda y eficiente para eliminar redundancia dentro de un conjunto de secuencias de proteínas o de DNA. El artículo donde se describe la primera versión de CD-HIT se publicó en Bioinformatics en 2006.


Aunque hay una serie de servidores web que pueden ayudarnos en esta tarea, si tienes un volumen realmente grande de secuencias lo más probable es que necesites descargar el código fuente (en C++) y compilarlo con make. Una vez compilado, el programa tiene múltiples opciones, pero su uso es muy sencillo.  El algoritmo se basa en calcular la identidad de todas las secuencias (previamente ordenadas por longitud) por parejas para ir generando clusters de secuencias con una identidad superior al umbral deseado. Por defecto el umbral de identidad es del 90%, medida a lo largo de toda la secuencia (global), no sólo la alineada. Por supuesto se puede modificar este comportamiento y elegir por ejemplo umbrales más bajos o identidades locales.


Terminaré con un ejemplo de uso. Por ejemplo, dado el mismo archivo FASTA de la entrada anterior de Álvaro, para obtener el subconjunto no redundante al 80% sería necesario el siguiente comando en el terminal Linux:

$ cd-hit/cd-hit -i redundante.faa -o noredundante.faa -c 0.8 

Como resultado, además del fichero noredundante.faa, se obtienen otros dos archivos (noredundante.faa.clstr y noredundante.faa.bak.clstr) que contienen
los clusters encontrados y las secuencias redundantes observadas.

5 de julio de 2010

Eliminar secuencias redundantes y transcritos repetidos en genomas y proteomas

Una situación cada vez más habitual en genómica y proteómica es la redundancia de datos. Con las modernas técnicas de secuenciación de alto rendimiento se generan gigas y gigas de datos que es necesario filtrar y corregir para obtener resultados inteligibles.

El problema más habitual y sencillo es la redundancia de secuencias, cuando en un mismo fichero aparecen secuencias idénticas con identificadores diferentes. En otras ocasiones el problema es el contrario, identificadores idénticos apuntan a secuencias parciales de una misma secuencia (especialmente en datos de proteómica o transcriptómica).

El siguiente script filtra un fichero de sequencias en formato FASTA, dos tipos de filtro son posibles: 1) filtrar las secuencias redundantes (repetidas), 2) filtrar las secuencias de transcritos parciales (recuperando únicamente los transcritos de mayor secuencia en ficheros de proteínas)

 #!/usr/bin/perl  
   
 my ($infile) = @ARGV;  
 my $outfile = filter_fasta_file($infile,['longtrans','nonredundant']);  
   
 # Create a FASTA file with custom filtered sequences  
 sub filter_fasta_file {  
   
      my ($infile, $filter, $outfile) = @_;  
   
      if (!defined($outfile)){$outfile=$infile.'.filtered'}  
   
      open(INFILE,"$infile");  
      open(OUTFILE,">$outfile");  
   
      my ($sequences,$names,$headers) = read_FASTA_file($infile);  
   
      my (@previous_sequences,@previous_names,@previous_headers);  
   
      # my (@filtered_sequences,@filtered_names,@fitered_headers);  
   
      for (my $i=0; $i<=$#{$headers}; $i++){  
   
           my $write=1;  
   
           if (in_array($filter,'nonredundant')){  
                my ($found,$posic)=in_array($sequences, $sequences->[$i], 1);  
                foreach my $pos (@{$posic}){  
                     if ($i != $pos && $i > $pos){  
                          $write = 0; last;  
                     }  
                }  
           }  
   
           if (in_array($filter,'longtrans')){  
                $names->[$i] =~ /([^\.]+)\.?(\d?)/;  
                my $name = $1;  
                my $variant = $2;  
                my ($found,$posic)=in_array($names, "/$1/", 1);  
                foreach my $pos (@{$posic}){  
                     if ($i != $pos){  
                          if (length($sequences->[$i]) < length($sequences->[$pos])) {  
                               $write = 0; last;  
                          }  
                          if (length($sequences->[$i]) == length($sequences->[$pos]) && $sequences->[$i] eq $sequences->[$pos]){  
                               $names->[$pos] =~ /([^\.]+)\.?(\d?)/;  
                               my $name2 = $1;  
                               my $variant2 = $2;  
                               if ($variant2 < $variant){  
                                    $write = 0; last;  
                               }  
                          }  
                     }  
                }  
           }  
   
   
           if ($write) {  
                print OUTFILE $headers->[$i]."\n".$sequences->[$i]."\n";  
           }  
   
      }  
   
      close INFILE;  
      close OUTFILE;  
   
      return $outfile;  
 }  
   
   
 #################################################################################  
   
 # Create 3 array references with sequences, names and headers from a FASTA file  
 sub read_FASTA_file  
 {  
      my($FASTAfile,$full_header,$allow_empty_sequences,$allow_multi_rows) = @_;  
   
      if(!$full_header){ $full_header = $allow_empty_sequences = 0 }  
      if(!$allow_empty_sequences){ $allow_empty_sequences = 0 }  
      if(!$allow_multi_rows){ $allow_multi_rows = 0 }  
   
      my (@sequences,@names,@headers,$seq,$seqname,$header);  
      my $headerOK = 0;  
   
      open(FASTA,$FASTAfile) || die "# read_FASTA_file : cannot find $FASTAfile\n";  
      $seq = '';  
      while(<FASTA>){  
           next if(/^#/); # allow comments  
           s/\012\015?|\015\012?|^\s+|\s+$//;  
           if(/>/){  
                $headerOK++;  
                if($seq || ($allow_empty_sequences && $seqname)){  
                     if(!$allow_multi_rows){  
                          $seq =~ s/\d|\s//g;  
                     }  
                     push(@sequences,$seq); # store this FASTA sequence  
                     push(@names,$seqname);  
                     push(@headers,$header);  
                     $seq='';  
                }  
   
                if($full_header){  
                     $seqname = substr($_,1);  
                }else { # get next sequence name  
                     if (/^>\s*([^\t|^\||^;]+)/){  
                          $header = $_;  
                          $seqname = $1;  
                          if (length($seqname)>32 && $seqname =~ /^([^\s]+)+\s/){  
                                    $seqname = $1;  
                          }  
                          $seqname =~ s/^\s+|\s+$// ;  
                     }  
   
                }  
           }else{  
                if(!$allow_multi_rows){   
                     s/\s+//g;  
                }
                $seq .= $_;  
           }  
      }  
      close(FASTA);  
   
      # take care of last FASTA sequence  
      push(@sequences,$seq);  
      push(@names,$seqname);  
      push(@headers,$header);  
   
      if(!$headerOK){ @sequences = @names = @headers = (); } #Aug2007  
   
      return (\@sequences,\@names,\@headers);  
 }  
   
 #################################################################################  
   
 # Searchs a string inside an array  
 sub in_array {  
      my ($array,$search_for,$search_posic) = @_;  
 #      my %items = map {$_ => 1} @$array; # create a hash out of the array values  
 #      my $exist = (exists($items{$search_for}))?1:0;  
 #      return (exists($items{$search_for}))?1:0;  
   
      my @posic;  
   
      # If search a pattern  
      if ($search_for =~ /^\/.+\/$/){  
           $search_for =~ s/^\/|\/$//g;  
           @posic = grep { $array->[$_] =~ /$search_for/ } 0 .. $#{$array};  
      } else {  
           @posic = grep { $array->[$_] eq $search_for } 0 .. $#{$array};  
      }  
   
      if (defined($search_posic) && $search_posic){  
           (scalar @posic)?return (1,\@posic):return (0);  
      } else {  
           (scalar @posic)?return 1:return 0;  
      }  
 }  
   

An example FASTA file to test the script:

 >AT1G51370.2 | Symbols: | F-box family protein | chr1:19045615-19046748 FORWARD  
 MVGGKKKTKICDKVSHEEDRISQLPEPLISEILFHLSTKDSVRTSALSTKWRYLWQSVPG  
 LDLDPYASSNTNTIVSFVESFFDSHRDSWIRKLRLDLGYHHDKYDLMSWIDAATTRRIQH  
 LDVHCFHDNKIPLSIYTCTTLVHLRLRWAVLTNPEFVSLPCLKIMHFENVSYPNETTLQK  
 LISGSPVLEELILFSTMYPKGNVLQLRSDTLKRLDINEFIDVVIYAPLLQCLRAKMYSTK  
 NFQIISSGFPAKLDIDFVNTGGRYQKKKVIEDILIDISRVRDLVISSNTWKEFFLYSKSR  
 PLLQFRYISHLNARFYISDLEMLPTLLESCPKLESLILVMSSFNPS*  
 >AT1G70790.1 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD  
 MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS  
 IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR  
 TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP  
 TKRLD*  
 >AT1G50920.1 | Symbols: | GTP-binding protein-related | chr1:18870555-18872570 FORWARD  
 MVQYNFKRITVVPNGKEFVDIILSRTQRQTPTVVHKGYKINRLRQFYMRKVKYTQTNFHA  
 KLSAIIDEFPRLEQIHPFYGDLLHVLYNKDHYKLALGQVNTARNLISKISKDYVKLLKYG  
 DSLYRCKCLKVAALGRMCTVLKRITPSLAYLEQIRQHMARLPSIDPNTRTVLICGYPNVG  
 KSSFMNKVTRADVDVQPYAFTTKSLFVGHTDYKYLRYQVIDTPGILDRPFEDRNIIEMCS  
 ITALAHLRAAVLFFLDISGSCGYTIAQQAALFHSIKSLFMNKPLVIVCNKTDLMPMENIS  
 EEDRKLIEEMKSEAMKTAMGASEEQVLLKMSTLTDEGVMSVKNAACERLLDQRVEAKMKS  
 KKINDHLNRFHVAIPKPRDSIERLPCIPQVVLEAKAKEAAAMEKRKTEKDLEEENGGAGV  
 YSASLKKNYILQHDEWKEDIMPEILDGHNVADFIDPDILQRLAELEREEGIREAGVEEAD  
 MEMDIEKLSDEQLKQLSEIRKKKAILIKNHRLKKTVAQNRSTVPRKFDKDKKYTTKRMGR  
 ELSAMGLDPSSAMDRARSKSRGRKRDRSEDAGNDAMDVDDEQQSNKKQRVRSKSRAMSIS  
 RSQSRPPAHEVVPGEGFKDSTQKLSAIKISNKSHKKRDKNARRGEADRVIPTLRPKHLFS  
 GKRGKGKTDRR*  
 >AT1G70795.1 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD  
 MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS  
 IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR  
 TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP  
 TKRLD*  
 >AT1G70790.2 | Symbols: | C2 domain-containing protein | chr1:26700724-26702127 FORWARD  
 MEDKPLGILRVHVKRGINLAIRDATTSDPYVVITLANQKLKTRVINNNCNPVWNEQLTLS  
 IKDVNDPIRLTVFDKDRFSGDDKMGDAEIDFRPFLEAHQMELDFQKLPNGCAIKRIRPGR  
 TNCLAEESSITWSNGKIMQEMILRLKNVECGEVELMLEWTDGPGCKGLGREGSKKTPWMP  
 TKRLD*  
 >AT1G36960.1 | Symbols: | unknown protein | chr1:14014796-14015508 FORWARD  
 MTRLLPYKGGDFLGPDFLTFIDLCVQVRGIPLPYLSELTVSFIAGTLGPILEMEFNQDTS  
 TYVAFIRVKIRLVFIDRLRFFRREEAAASNTITDQTHMTSSNSSDISPASPISQPPLPAS  
 LPSHDSYFDAGIQASRLVNPRAISQHHFSSSYSDFKGKEKAKIKIGECSKRKKDKQVDSG  
 T*  
 >AT1G51370.1 | Symbols: | F-box family protein | chr1:19045615-19047141 FORWARD  
 MVGGKKKTKICDKVSHEEDRISQLPEPLISEILFHLSTKDSVRTSALSTKWRYLWQSVPG  
 LDLDPYASSNTNTIVSFVESFFDSHRDSWIRKLRLDLGYHHDKYDLMSWIDAATTRRIQH  
 LDVHCFHDNKIPLSIYTCTTLVHLRLRWAVLTNPEFVSLPCLKIMHFENVSYPNETTLQK  
 LISGSPVLEELILFSTMYPKGNVLQLRSDTLKRLDINEFIDVVIYAPLLQCLRAKMYSTK  
 NFQIISSGFPAKLDIDFVNTGGRYQKKKVIEDILIDISRVRDLVISSNTWKEFFLYSKSR  
 PLLQFRYISHLNARFYISDLEMLPTLLESCPKLESLILEMVKNQSTRRHGEKREPNVMVS  
 TVPWCLVSSLKFVELKRSIPRYEGEMELVRYVLTNSTVLKKLRLNVYYTKKAKCAFLTEL  
 VAIPRCSSTCVVLVL*