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

Genotipado HLA y software disponible

En la presente entrada intentaré hacer una recopilación de software para el genotipado de HLA (MHC humano), en la sección de comentarios podéis añadir otras herramientas que intentaré incorporar al texto. Para escribir la entrada me resultó muy útil el siguiente post en Biostar forum y Omictools.

Antes de empezar es necesario mencionar que existe una base de datos muncial de alelos HLA denominada IMGT-HLA que recopila las miles de secuencias conocidas (y públicas) de genes y transcritos para esta familia. Todas las herramientas de genotipado emplearán las secuencias de esta base de datos para realizar sus predicciones.

Número de alelos registrados para cada tipo de HLA en la base de datos IMGT-HLA (Septiembre 2015).

Las herramientas de genotipado de HLA alinean (map) las reads de NGS a las secuencias de referencia de los principales loci HLA humanos presentes en la mencionada base de datos (clase I: HLA-A, HLA-B y HLA-C, clase II: HLA-DR y HLA-DQ).

Dependiendo del software, se puede procesar reads de secuenciación genómica, exómica o transcriptómica, aunque cuando se quieren analizar cientos/miles de individuos en un único experimento se suelen enriquecer por la técnica de secuenciación de amplicones (diseñando primers que amplifican las regiones menos conservadas, ver entrada anterior) o captura con sondas específicas para HLA.

El mapeo de reads a las secuencias de referencia puede realizarse directamente o tras realizar un ensamblaje de novo previo. Realizando un ensamblaje previo será más fácil encontrar alelos únicos de HLA puesto que los contigs resultantes darán menos mapeos ambiguos. Sin embargo el ensamblaje de esta familia de genes parálogos generará también un gran número de ensamblajes erróneos o quimeras (falsos contigs mezcla de dos secuencias análogas). A su vez el mapeo directo de reads puede generar ambigüedades puesto que numerosas reads alinearán con múltiples referencias a la vez.

Típica estrategia de genotipado por mapeado (alineamiento) de reads a sequencias de referencia. A la izquierda las reads son ensambladas de novo antes de ser alineadas. A la derecha las reads son directamente alineadas. Imagen modificada de Warren et al. (2012).

Listado de software para el genotipado de HLA

Sólo se listan herramientas libres para uso académico ordenadas por orden cronológico de la última versión del software:

• seq2HLA (Jun 2015): diseñado para procesar reads de RNA-Seq, mapea las mismas a las secuencias alélicas de referencia (IMGT-HLA) generando genotipos con una puntuación de probabilidad para los mismos y los niveles de expresión de los alelos predichos.
• HLAreporter (May 2015): primero filtra las reads que mapean a los diversos alelos de un único gen, las ensambla de novo y los contigs resultantes son de nuevo alineados a los alelos de referencia iniciales para asignar genotipos.
• HLAminer (Feb 2015): realiza un ensamblaje de novo de las reads (de casi cualquier procedencia) para después alinear los contigs resultantes contra los alelos de referencia.
• Optitype (Apr 2014): otro método que acepta diversos tipos de datos y también se basa en el mapeo a secuencias exónicas de referencia. Los resultados del mapeo son representados en forma matricial, las reads en filas y los alelos en columnas. En la matriz se identifican como máximo 2 alelos que explican el mayor número de reads mapeadas.
• PHLAT (Feb 2014): además de analizar datos genómicos, transcriptómicos y exómicos, ha sido también testado con datos de amplicones. Mapea reads a las secuencias de referencia seleccionando múltiples alelos candidatos y selecciona la pareja de alelos con la mayor probabilidad de acontecer juntos.
• HLAforest (Jan 2013): similar a seq2HLA, analiza reads de RNA-Seq, aunque puede ser usada con otro tipo de datos reduciendo su precisión.
• ATHLATES (Jun 2012): similar a HLAminer, filtra y ensambla las reads para después identificar exones de IMGT-HLA en los contigs ensamblados. Está diseñado para reads de sequenciación de exoma.
• GATK-HLA Caller (Dec 2011): similar a seq2HLA, alinea, filtra y calcula probabilidades para cada genotipo.

Por último explicaré el software diseñado en mi laboratorio...

• AmpliHLA (Sep 2015), no es el mejor, simplemente es diferente. Está únicamente enfocado al análisis de datos de secuenciación de amplicones usando primers que amplifiquen diferentes regiones de los genes HLA de interés y etiquetas de DNA que diferencien las muestras.

AmpliHLA requiere un pre-procesado online de los datos de NGS con la herramienta AmpliSAS. AmpliSAS clasifica las reads por muestra/individuo, corrige errores de secuenciación y filtra artefactos de secuenciación y PCR. AmpliSAS está diseñado para el genotipado de cualquier tipo de gen, especialmente si no tenemos alelos de referencia previos (como generalmente ocurre con muchos organismos cuyo genoma no ha sido secuenciado o regiones complejas del genoma como los genes que codifican las moléculas de MHC).

Un archivo Excel generado tras el análisis con AmpliSAS ha de ser introducido en el formulario de AmpliHLA y el programa automáticamente unificará marcadores (diversas regiones amplificadas de un mismo gen) y buscará sus secuencias en la base de datos humana para genotipar con la máxima precisión posible cada individuo. El principal inconveniente es el requerimiento de múltiples PCRs y diversos marcadores por gen para conseguir un genotipado de calidad. La principal ventaja es la obtención de un genotipado de-novo que permite descubrir alelos incluso si no están presentes en la base de datos humana.

Esquema de funcionamiento de la herramienta de genotipado de novo mediante secuenciación de amplicones: AmpliSAS. Primero las reads son separadas por muestras y marcadores. Después se realiza un clustering de los errores de secuenciación con sus secuencias de origen. Por último se filtran las reads minoritarias y se asignan alelos por cada muestra y marcador.