Bioinformática en C (y perl6 async)

Hola, y feliz año!
tras la última entrada, y antes de que escriba algo más completo sobre perl6,
parece que la gestión de tareas asíncronas tiene mucho potencial...

Pero vamos a lo que vamos. Hoy quería compartir un documento sobre cómo programar en C con los compiladores actuales, porque algunos todavía seguimos usando sintaxis anticuada y esto es relevante en nuestro campo, dado que hay muchas aplicaciones bioinformáticas escritas en C, normalmente por cuestiones de eficiencia.

El documento en cuestión es: https://matt.sh/howto-c

y ofrece un montón de trucos y recomendaciones para modernizar y simplificar la escritura de código en C, con ejemplos sencillos como éste para vectores dinámicos:

 uintmax_t arrayLength = strtoumax(argv[1], NULL, 10);
    void *array[arrayLength];

 /* Hace innecesario liberar el puntero al final.

    Ojo: debe ser un tamaño razonable para el sistema donde se ejecute. */

Ha sido traducido al español aquí y comentado/criticado aquí y en menéame,
un saludo,
Bruno

Despedimos el 2015, viene Perl 6

Buenas,
quedan ya pocos días para terminar el año, o lo que es lo mismo, pocos días para que podamos probar la primera versión oficial de Perl 6, tras demasiados (15) años de desarrollo, de la mano del compilador Rakudo. Si no os suena nada de esto, podéis empezar por leer las preguntas frecuentes o directamente empezar tutoriales como perl6intro.

Ha habido mucho morbo con la muerte de Perl, y seguirá habiendo a lo largo de 2016 discusiones sobre la compatibilidad de perl 6 y 5, el que usamos ahora, sobre la eficiencia del nuevo lenguaje, o sobre la portabilidad de CPAN a perl 6. Todo esto habrá que ir viendo como se desarrolla muy pronto, pero desde luego es una noticia importante para los usuarios de Perl, esperemos que buena!

Un saludo y los mejores deseos para el próximo año, long live perl!
Bruno

Post-datas del 29 de Diciembre: 

1) Podéis ampliar esto en [http://perlweekly.com/archive/231.html]

2) Para instalar rakudo vía git debes tener abierto el puerto 9418; en caso contrario es posible hacerlo vía https añadiendo estas líneas a tu archivo ~/.gitconfig: 

[url "https://"]
    insteadOf = git://
[url "https://github.com"]
    insteadOf = git@github.com

Post-data del 1 de Enero: 

1) Tutorial exprés en [https://learnxinyminutes.com/docs/perl6]
2) Si ya usas perl5 [http://doc.perl6.org/language/5to6-nutshell]

Post-data de Junio de 2016: 

Mientras no sea pueda instalar perl6 directamente en Ubuntu prefiero usar rakudobrew para instalar perl6 como un ecosistema aislado, como se explica en: [http://rakudo.org/how-to-get-rakudo/#Installing-Rakudo-Star-Source-Rakudobrew]

Si los médicos tuvieran que trabajar de bioinformáticos...

Hoy os publico una viñeta que me ha hecho gracia y se adapta muy bien a la bioinformática... desgraciadamente se cumple totalmente, sobretodo lo del conserje...

El original aquí: http://sinergiasincontrol.blogspot.com/2008/10/33-suposiciones.html

Anotando datos del NCBI de forma automática con un script de Perl

Esta semana me di cuenta que un antiguo script de Perl que leía datos bibliográficos de Pubmed automáticamente había dejado de funcionar.

La explicación es que el NCBI ha dejado de dar soporte a su servicio web con el protocolo SOAP desde el 1 de julio de 2015. No soy un experto informático, pero hace unos años ese protocolo estaba de moda y ahora parece que ya no (imagino que por razones de seguridad).

He aquí el código OBSOLETO que da el error:

 use SOAP::Lite;  
   
 my $pubmed_id = '24234003';  
   
 my $WSDL = 'http://www.ncbi.nlm.nih.gov/soap/v2.0/efetch_pubmed.wsdl';  
 my $soap = SOAP::Lite->service($WSDL)  
                ->on_fault( sub { my($self, $res) = @_;  
                     die  
                     "faultcode:", ref $res ? $res->faultcode : $self->transport->status, "\n" ,  
                     "faultstring:", ref $res ? $res->faultstring : $self->transport->status, "\n";  
                });  
 my $response = $soap->run_eFetch(     SOAP::Data->name(id => $pubmed_id),  
                          SOAP::Data->name(db => "pubmed")  
                     );  
   
 my %results;  
   
 if (ref($response) eq "HASH") {  
      $results{'title'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'ArticleTitle'};  
      $results{'journal'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Journal'}{'Title'};  
      $results{'pubmed'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'PMID'};  
      $results{'volume'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Journal'}{'JournalIssue'}{'Volume'};  
      $results{'issue'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Journal'}{'JournalIssue'}{'Issue'};  
      $results{'year'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Journal'}{'JournalIssue'}{'PubDate'}{'Year'};  
      $results{'pages'} = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'Pagination'}{'MedlinePgn'};  
      $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$results{'pubmed'};  
      if (ref($response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'AuthorList'}{'Author'}) eq "ARRAY") {  
           my @authors;  
           foreach my $author (@{$response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'AuthorList'}{'Author'}}){  
                my $initials = join('.',split('',$author->{'Initials'}));  
                push(@authors, $author->{'LastName'}." ".$initials);  
           }  
           $results{'authors'} = join(", ",@authors);  
      } elsif (ref($response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'AuthorList'}{'Author'}) eq "HASH") {  
           my $author = $response->{'PubmedArticle'}{'MedlineCitation'}{'Article'}{'AuthorList'}{'Author'};  
           my $initials = join('.',split('',$author->{'Initials'}));  
           $results{'authors'} = $author->{'LastName'}." ".$initials;  
      }  
      foreach my $key (keys %{$results}){  
           if (!defined($results{$key})) {  
                $results{$key} = '';  
           }  
      }  
 }  
   
 print %results;  
   

Ahora las consultas automáticas al NCBI se han de realizar insertando nuestra consulta en una URL e indicando el formato deseado de respuesta. Por ejemplo, la consulta equivalente al código SOAP anterior se puede realizar con la siguiente URL:

• http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=pubmed&id=24234003 

 O si preferimos el resultado en XML:

• http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=pubmed&id=24234003&retmode=xml

De esta forma, el código anterior queda reescrito de la siguiente forma:

 use LWP::Simple;  
   
 my $pubmed_id = '24234003';  
   
 my $eutils_url = 'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/';  
 my $eutil = 'efetch';  
 my $db = 'pubmed';  
 my $type = 'docsum';  
 my $mode = 'text';  
   
 my $url = sprintf('%s%s.fcgi?db=%s&id=%s&rettype=%s&retmode=%s', $eutils_url, $eutil, $db, $pubmed_id, $type, $mode);  
   
 my $data = get($url);  
   
 my %results;  
   
 if (defined($data)){  
      $data =~ s/\012\015?|\015\012?/ /g;  
      if ($data =~ /\d+\:\s*(.+?)\.(.+?)\.(.+?)\.(.+?)\s.+?\;(\d+)\((\d+)\)\:([\d\-]+)./){  
           @{$results{'authors'}} = split(/,\s?/,$1);  
           $results{'title'} = $2;  
           $results{'journal'} = $3;  
           $results{'year'} = $4;  
           $results{'volume'} = $5;  
           $results{'issue'} = $6;  
           $results{'pages'} = $7;  
           $results{'pubmed'} = $pubmed_id;  
           $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$pubmed_id;  
      } elsif ($data =~ /\d+\:\s*(.+?)\.(.+?)\.(.+?)\.(.+?)\s/){  
           @{$results{'authors'}} = split(/,\s?/,$1);  
           $results{'title'} = $2;  
           $results{'journal'} = $3;  
           $results{'year'} = $4;  
           $results{'pubmed'} = $pubmed_id;  
           $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$pubmed_id;  
      }  
 }  
   
 while (my ($key, $value) = each %results) {  
      if (ref($value) ne ARRAY){  
           printf("%s: %s\n", $key, $value);  
      } else {  
           printf("%s: %s\n", $key, join(', ',@$value));  
      }  
 }  
   

Para otro tipo de consultas automáticas (secuencias, artículos científicos...) recomiendo leer la documentación oficial del NCBI, donde se explica la utilidad de las diferentes Entrez Programming Utilities (E-utilities): ESearch, EFetch, ESummary, EPost y ELink.

Seguro que alguno de vosotros está pensando porqué no usar BioPerl o BioPython para estas consultas. La respuesta es porque quería explicar el fin de la era SOAP y la estructura de las nuevas consultas vía URL.

Aquí el código en BioPerl:

 use Bio::DB::EUtilities;  
   
 my $pubmed_id = '24234003';  
   
 my $eutil = 'efetch';  
 my $db = 'pubmed';  
 my $type = 'docsum';  
 my $mode = 'text';  
 my $email = 'mymail@foo.bar';  
    
 my $factory = Bio::DB::EUtilities->new(     -eutil  => $eutil,  
                          -db   => $db,  
                          -id   => $pubmed_id,  
                          -email  => $email,  
                          -rettype => $type,  
                          -retmode => $mode  
                     );  
    
 my $data = $factory->get_Response->content;  
   
 my %results;  
   
 if (defined($data)){  
      $data =~ s/\012\015?|\015\012?/ /g;  
      if ($data =~ /\d+\:\s*(.+?)\.(.+?)\.(.+?)\.(.+?)\s.+?\;(\d+)\((\d+)\)\:([\d\-]+)./){  
           @{$results{'authors'}} = split(/,\s?/,$1);  
           $results{'title'} = $2;  
           $results{'journal'} = $3;  
           $results{'year'} = $4;  
           $results{'volume'} = $5;  
           $results{'issue'} = $6;  
           $results{'pages'} = $7;  
           $results{'pubmed'} = $pubmed_id;  
           $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$pubmed_id;  
      } elsif ($data =~ /\d+\:\s*(.+?)\.(.+?)\.(.+?)\.(.+?)\s/){  
           @{$results{'authors'}} = split(/,\s?/,$1);  
           $results{'title'} = $2;  
           $results{'journal'} = $3;  
           $results{'year'} = $4;  
           $results{'pubmed'} = $pubmed_id;  
           $results{'url'} = 'http://www.ncbi.nlm.nih.gov/pubmed/'.$pubmed_id;  
      }  
 }  
   
 while (my ($key, $value) = each %results) {  
      if (ref($value) ne ARRAY){  
           printf("%s: %s\n", $key, $value);  
      } else {  
           printf("%s: %s\n", $key, join(', ',@$value));  
      }  
 }  
   

Y el código en BioPython:

 from Bio import Entrez  
   
 pubmed_id = '24234003';  
   
 handle = Entrez.efetch(db='pubmed', id=pubmed_id, rettype='docsum', retmode="xml", email='mymail@foo.bar')  
   
 results = Entrez.read(handle)  
   
 for key, value in results[0].iteritems():  
      print "%s: %s" % (key,value)  
   

merge TSV files as Excel sheets (multitab2xls.pl)

Buenas,
hoy me he vuelto a encontrar con una tarea reincidente, la de juntar varios ficheros con datos separados por tabuladores (formato TSV) en un sólo fichero Excel para compartirlo más fácilmente con otros colegas.  Como decía esto me pasa con cierta frecuencia, y por eso hoy dediqué 5 minutos a escribir unas líneas de Perl que permiten automatizar esta tarea, previa instalación del módulo Spreadsheet::WriteExcel .

Para instalar dicho módulo es suficiente con:

$ sudo cpan -i  Spreadsheet::WriteExcel

El código del script multitab2xls.pl es el siguiente:

 #!/usr/bin/perl -w  
 use strict;  
 use Spreadsheet::WriteExcel;  
 use File::Basename;  
 die "# usage: $0 <target.xls> <*.tab> <*.tsv>\n" if(!$ARGV[1]);  
 my ($xlsname,$files) = (shift(@ARGV),0);  
 my $xlsbook = Spreadsheet::WriteExcel->new($xlsname);  
 my $title_sheet = $xlsbook->add_worksheet('original files');  
 foreach my $tabfile (@ARGV) {  
     $title_sheet->write($files++,0,$tabfile);  
   my $sheet = $xlsbook->add_worksheet(substr(basename($tabfile,qw(.tsv .tab)),-30));  
     open(TAB,$tabfile) || die "# cannot read $tabfile, exit\n";  
     while(<TAB>) {  
         $sheet->write_row($.-1, 0, [split(/\t/)]);  
     }  
     close(TAB);  
 }  

Y se ejecutaría de esta manera:

$ perl multitab2xls.pl outfile.xls folder/*.tab 

Un saludo,
Bruno