22 de diciembre de 2015

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]

25 de noviembre de 2015

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

3 de noviembre de 2015

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:
 O si preferimos el resultado en 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)  
   



14 de octubre de 2015

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

28 de septiembre de 2015

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.

14 de septiembre de 2015

Secuenciación de amplicones y genotipado de alto rendimiento

Secuenciación de amplicones (SA)  es una traducción aproximada al español de la técnica de "Amplicon sequencing" que junto con las tecnologías de secuenciación masiva (del inglés new generation sequencing, NGS) permite genotipar cientos/miles de individuos en un único experimento.

La secuenciación de amplicones (SA) consiste en secuenciar los productos de múltiples PCRs. Un amplicón se define como el conjunto de secuencias obtenidas de cada PCR individual.

Antiguamente se realizaban PCRs individuales y se secuenciaban uno a uno los productos. Con las nuevas técnicas de NGS, podemos incluir etiquetas de DNA (por ej. una secuencia única de 6 nucleótidos) diferentes para cientos de muestras o individuos y clasificar más tarde las secuencias o reads resultantes de una única secuenciación (Binladen et al. 2007; Meyer et al. 2007).
Esquema de etiquetado y amplificación para la secuenciación de amplicones.

Mediante esta técnica podremos genotipar individuos y distinguir los diferentes alelos (con la secuenciación tradicional a veces es complicado separar alelos de un mismo gen). El principal problema de las técnicas de NGS es su alta tasa de error, que a su vez puede ser compensada incrementando la profundidad de secuenciación (el número de reads). Otros problemas pueden ser los errores de la polimerasa o la generación de quimeras (una secuencia mezcla de otras).

En el siguiente enlace podemos ver un vídeo explicativo del proceso de secuenciación de amplicones:
http://www.jove.com/video/51709/la-secuenciacin-de-prxima-generacin-de-16s-arn-ribosomal-genes?language=Spanish

Básicamente existen 4 etapas en el análisis por AS con NGS:
  1. Diseño experimental de los primers usados para amplificar los genes de interés (marcadores) y las etiquetas a usar para distinguir los diferentes individuos o muestras.
  2. Amplificación por PCR de los marcadores en el laboratorio, generalmente se realiza una PCR por cada muestra.
  3. Secuenciación de los productos de amplificación. Las tecnologías de NGS más usadas para SA son: Illumina, 454 e Ion Torrent.
  4. Análisis bioinformático de los datos de secuenciación. El análisis incluye separación de las reads en amplicones, corrección de errores de secuenciación, filtrado de reads minoritarias/contaminantes y generación de genotipados.
Etapas de la técnica de secuenciación de amplicones mediante NGS.

Aplicaciones 


SA es utilizado para realizar clasificaciones taxonómicas usando genes como: cytochrome c oxidase subunit 1 (CO1), genes rRNA (16S/18S/28S), genes específicos de plantas (rbcL, matK, and trnH-psbA) y espaciadores internos nucleares (ITSs) (Kress et al. 2014; Joly et al. 2014). Los genes anteriores se distinguen por una tasa de mutación suficientemente rápida como para distinguir especies cercanas y a la vez suficientemente estables como para distinguir congéneres.

Lista de genes habitualmente usados como marcadores taxonómicos (fuente: Kress et al. 2014)

Un experimento pionero de SA fue la determinación de la diversidad microbiana en aguas marinas profundas (Sogin et al. 2006), usando primers flanquenado la región V6 hipervariable de la subunidad 16S rRNA bacteriana. Dicho estudio descubrió miles de poblaciones minoritarias de organismos no conocidos con anterioridad.

Otro gran campo de aplicación es el genotipado de familias de genes de alta complejidad, como el complejo mayor de histocompatibilidad, que poseen múltiples loci y diferente número de copias entre individuos, incluso de la misma especie (Babik et al. 2010; Lighten et al. 2014). El complejo mayor de histocompatibilidad (MHC) de clase I y II codifica receptores celulares que presentan antígenos a las células del sistema inmune y son los genes más polimórficos conocidos en vertebrados . El MHC humano también se conoce como HLA (Human Leukocyte Antigen) y juega un papel clave en la compatibilidad en el transplante de órganos. Los loci del MHC son tan polimórficos que no hay dos individuos en una población no endogámica que posean el mismo conjunto de alelos (excepto gemelos).

Estadísticas del número de alelos conocidos para la familia de genes del HLA (fuente: base de datos IMGT-HLA)
Hasta hace poco era necesario clonar y secuenciar uno por uno los diferentes alelos de este tipo de genes para conseguir una secuencia fiable. Actualmente tan tediosa tarea puede ser simplificada mediante un único experimento de NGS que incluya múltiples individuos y múliples genes. El secuenciador de nueva generación (ej.: Illumina, 454 o Ion Torrent) leerá las sequencias individuales de cada uno de los alelos. Actualmente existen incluso kits comerciales para simplificar el proceso: Illumina TruSeq Custom Amplicon, Roche 454 Fluidigm Access Array or Life Technologies Ion Torrent Ion AmpliSeq.

13 de agosto de 2015

rperl: de perl a c++

Hola,
el proyecto RPerl acaba de liberar la primera versión en CPAN tras más de dos años de desarrollo en https://github.com/wbraswell/rperl . Su mascota es un correcaminos y a su autor principal (Will Braswell), a tenor de lo que escribe en la página del proyecto, parece que le gusta la literatura fantástica o las películas de semana santa.
 
Volviendo a lo importante, la filosofía del proyecto es hacer un compilador que permita convertir código Perl, siguiendo una versión limitado de su sintaxis, en código C++ que se compila y enlaza con Inline::CPP, del que ya hablamos en otra entrada.

Apenas he hecho algunas pruebas, os dejo mis notas:
  1.  Para instalar RPerl hay instrucciones detalladas en: https://github.com/wbraswell/rperl/blob/master/INSTALL. En mi caso no funcionaron a la primera, pero tras actualizar g++ a la versión 4.7 Io logré con:

    $ cpan -i  Inline::C
    $ cpan -if Inline::CPP
    $ cpan -i RPerl 
    

  2. Los ejemplos más sencillos de uso los encontré en: https://github.com/wbraswell/rperl/tree/master/lib/RPerl/Learning
  3. Para encontrar otros ejemplos más completos tuve que rebuscar en: https://github.com/wbraswell/rperl/tree/master/lib/RPerl/Test. Por ejemplo, para ver la sintaxis RPerl para escribir en un archivo encontré esto, o esto otro para expresiones regulares.
Como ejemplo de la ganancia en velocidad por usar tipos de datos estáticos y algoritmos compilados, los autores muestran un one-liner que ordena 5000 enteros:

$ perl -e 'use RPerl::Algorithm::Sort::Bubble; my $a = [reverse 0 .. 
5000]; use Time::HiRes qw(time); my $start = time; my $s = 
integer_bubblesort($a); my $elapsed = time - $start; print Dumper($s); 
print "elapsed: " . $elapsed . "\n";'

Este experimento, en sus manos, tarda 15s con Perl interpretado y 0.045s con RPerl precompilado. Por tanto parece que para algunas tareas numéricas vale pena complicar un poco el código Perl para obtener estas ganancias, no?

Me queda por ver si también vale la pena en otras tareas habituales en nuestro campo como leer y procesar archivos de gran tamaño y consultar luego datos extraído por medio de hashes, pero eso queda para otro día,
hasta luego,
Bruno



6 de julio de 2015

Alineamiento de secuencias de proteína y filogenias

Hola,
como anunciábamos hace unas semanas, esta semana participamos en un curso de verano de la Universidad de Zaragoza que empieza hoy en la ciudad de Jaca, al pie mismo del Pirineo. Por esta razón hemos preparado un material sobre Alineamiento de secuencias de proteína y filogenias que hoy colgamos en la Red para que podáis consultarlo si a alguien le interesa:

http://eead.csic.es/compbio/material/alineafilog

y que también podés descargar en formato PDF en:

http://digital.csic.es/handle/10261/117608


El índice del curso es el siguiente:
  •  Análisis jerárquico de la estructura de proteínas
    • Estructura primaria
    • Estructura secundaria
    • Estructura terciaria y cuaternaria
      • Demarcación de dominios en proteínas

Un saludo,
hasta pronto,
Bruno

13 de mayo de 2015

Curso de Python para biólogos - Lección 7. Lectura y escritura de archivos

Con Python podemos leer, crear y modificar archivos de texto. Por ejemplo podemos leer un fichero con miles/millones de secuencias de DNA y extraer información del mismo, lo cual sería muy complicado de forma manual. Se verá un ejemplo al final de la lección.

Trabajando con directorios

Antes de empezar a crear y leer archivos, vamos a consultar el directorio de trabajo por defecto de Python con el comando 'os.getcwd', crear un nuevo directorio para trabajar en esta lección con 'os.makedirs' y cambiar el directorio de trabajo por este nuevo directorio ('os.chdir'). Para realizar todo ello primero tenemos que importar el módulo 'import os' (Operating System) que nos proporcionará los métodos mencionados. Un módulo es una extensión de Python y se debe importar para poder trabajar con sus herramientas.

Creación de archivos

Una vez que estamos en el directorio de trabajo deseado, vamos a crear un archivo llamado 'example.txt', para ello usaremos la función 'open' con la opción 'w' que indica escritura (write). El nuevo archivo será un objeto guardado en la variable 'f' y con el método 'write' podremos insertar un texto en el mismo. No debemos olvidar cerrar siempre el archivo después de leer o modificar sus contenidos con el método 'close'. Si abrimos el archivo con nuestro explorador de archivos, veremos el texto que contiene.

Usando un bucle 'for' podemos insertar automáticamente múltiples líneas de texto:

Lectura de archivos

Para leer un archivo, primero debemos abrirlo con la función 'open' y la opción 'r' que indica lectura (read). El método 'read' leerá todo el archivo de una vez si no se especifica ningún argumento. Importante, después de leerlo debemos cerrarlo con la función close.

Una  forma más adecuada de leer un fichero puede ser línea a línea con un bucle 'for', cada línea será almacenada en la variable 'line' en cada iteración.

Existen varias formas alternativas de leer los ficheros, entre ellas los métodos 'read' y 'readline', ambos aceptan especificar el número de caracteres que serán leídos. En próximas lecturas, Python continuará en la posición que terminó la anterior.

Ejercicio. Contar el número de genes codificantes que tiene el genoma de Escherichia coli

Para ello, primero descargaremos el fichero del proteoma de E. coli del siguiente enlace. Después, copiaremos el fichero descargado a nuestro directorio de trabajo. Y por último escribiremos un código en Python que cuente el número de veces que aparece el símbolo '>' al comienzo de una línea, dicho símbolo indica el comienzo de una proteína.

Si hemos hecho todo bien, el número de genes debería ser '4140', o similar (las nuevas versiones del genoma pueden variar ligeramente dicho número). Podemos confirmarlo en la web de KEGG.

Próxima lección

 En la próxima lección se hará una introducción a las expresiones regulares.

Curso de Python para biólogos - Lección 6. Bucles 'for'

En la lección de hoy presentaremos los bucles 'for', muy similares a los 'while' explicados en la lección 4, pero con un código más sencillo.

Bucles 'for':

Un bucle 'for', al igual que 'while', repite la ejecución de un bloque de código un número determinado de veces. Si invocamos el bucle 'for' con un nombre de variable más 'in' y una lista, el número de repeticiones vendrá determinado por el número de elementos de la lista, que serán pasados de uno en uno a la variable durante cada iteración. Veamos algunos ejemplos:

Con la función 'range' podemos especificar un listado de números que serán pasados a la variable especificada en el bucle. El primer parámetro de la función será el número inicial en el bucle, el segundo parámetro indicará el último número y el tercer parámetro proporcina el incremento a aplicar en cada iteración, por defecto el incremento es 1.

Como ya se ha comentado 'for' es una simplificación de 'while' cuando se trabaja con listas o números. El mismo código con 'while' es posible, pero será más largo y complejo.

Veamos como podemos usar 2 bucles 'for' anidados para recorrer los elementos de una matriz:

Sentencias de control de bucles: 'break', 'continue' y 'pass'

A veces nos interesará salir de un bucle (ya sea 'while' o 'for') antes de terminar todas las iteraciones, la forma de conseguirlo es mediante las sentencias de control 'break' y 'continue'. 'Break' terminará totalmente el bucle y continuará la ejecución del código del resto de programa. 'Continue' terminará la ejecución de una iteración y pasará directamente a la siguiente iteración del bucle, sin salir del bucle. 'Pass' no hace nada, es simplemente una sentencia que puede ser usada cuando no se requiere ejecutar ninguna acción pero es requerido escribir algo.

Bucles 'for' y diccionarios:

Los bucles 'for' también pueden ser útiles para procesar una a una las claves de un diccionario, sus valores o ambas cosas a la vez.


Ejercicios: 

Volvamos a los ejercicios propuestos en la lección anterior y resueltos con bucles 'while', veamos cómo pueden resolverse de una forma más sencilla mediante bucles 'for'.



Próxima lección

En la próxima lección se explicará como leer y escribir archivos con Python.



24 de abril de 2015

Diferentes formas de poblar una lista en R

Hola,

trabajando con Gene Ontologies (GOs) he empezado a hacer unas pruebas con una biblioteca de Bioconductor llamada topGO. Esta sirve para comprobar si los términos de GOs de un conjunto de genes son diferentes de los que encontramos en otro conjunto.

Por ejemplo, si tenemos los genes de un organismo y nos interesan los que en nuestro tratamiento aumentan su expresión, podemos comprobar si los términos GO asociados a estos últimos son diferentes de los que tenemos para el conjunto de todos los genes del organismo.

Pero uno de los pasos previos, antes de hacer tests estadísticos con topGO, es asociar nuestros genes a términos GO. En Bioconductor y otras bases de datos hay listas ya precalculadas que podemos utilizar. Sin embargo, por distintos motivos (por ejemplo, al trabajar con un organismo para el cual no está tan estudiada la anotación funcional), podemos querer cargar nuestras propias asociaciones para que topGO trabaje sobre ellas.

Ya que un gen puede tener asociado más de un término GO, los desarrolladores de topGO han decidido que la estructura de datos a utilizar sea una lista de vectores de caracteres (técnicamente una "named list of character vectors"). Dicho de otra forma, una lista de genes, con cada una de las entradas de la lista asociada a un vector de nombres de términos GO. Una opción, si queremos trabajar con nuestra propia anotación, es utilizar la función 'readMappings' de topGO, que crea esta estructura si le damos un formato de fichero adecuado:

"It consists of one line for each gene with the following syntax:
gene_ID/TAB/GO_ID1, GO_ID2, GO_ID3, ...."

La otra opción es que creemos nosotros mismos la lista de vectores a partir de los datos que queramos. A veces, cuesta menos transformar un fichero de datos a otro formato (y usar la función 'readMappings', por ejemplo), y otras es más conveniente mantener nuestro formato de fichero y manipularlo como estructura de datos, sin crear otro fichero. En este caso, es además una buena oportunidad para probar algunas estructuras de datos y sentencias de R, que es lo que haremos nosotros.

¿Cómo es el formato de datos que tenemos en el fichero inicial? Obviamente no es el formato requerido por topGO (pues usaríamos la función 'readMappings' sin más). Nosotros vamos a partir de una tabla con 2 columnas (separadas por tabulador): una para el identificador del gen, la otra para un identificador de término GO. De esta forma, los distintos términos GO de cada gen quedan realmente en varias líneas:

gene_1    GO:000001
gene_1    GO:000004
gene_2    GO:000003
...

Esto obviamente puede cambiar según el origen de los datos, con qué programa se hayan obtenido las anotaciones, etc.

No vamos a entrar en detalles aquí sobre cómo leer el fichero con R o en dar ejemplos concretos de los datos. Es muy fácil generarlos y no es el objetivo de esta entrada. Lo importante para entender el código es que vamos a trabajar con la salida de read.table, en la variable 'go_tab', que será de tipo data.frame. En 'go_tab', los campos de las dos columnas se llaman 'identifier' (la columna con genes) y 'GOterm' (la columna con términos GO).

Vamos a ver varias alternativas para transformar éste 'data.frame' en una lista de vectores. En primer lugar, probaremos creando directamente una lista que iremos poblando con un bucle 'for'. Éste código es muy legible e intuitivo y programadores de muchos lenguajes diferentes pueden entenderlo fácilmente con solo echarle un vistazo.

El código completo de éste fragmento:

 gene2GO = list()  
 for (i in seq(1:nrow(go_tab))) {  
  identifier = go_tab$identifier[i]  
  go_term = go_tab$GOterm[i]  
  if (identifier %in% names(gene2GO)){  
   gene2GO[[identifier]] = c(gene2GO[[identifier]], as.character(go_term))  
  } else {  
   gene2GO[[identifier]] = as.character(go_term)  
  }  
 }  

En éste ejemplo, por tanto, vamos a recorrer la tabla de datos (go_tab) importada del fichero, e iremos construyendo nuestra lista final (gene2GO) concatenando nuevos términos GO para cada gen:

 gene2GO[[identifier]] = c(gene2GO[[identifier]], as.character(go_term))  

Vemos aquí una de las peculiaridades de R en el manejo de listas, que es el uso de corchetes dobles '[['. Si hay dudas de por qué se usa aquí '[[' en lugar de '[', la explicación rápida es que en R para acceder a datos en una lista hay que usar '[['. Para ampliar esto, está detallado en http://adv-r.had.co.nz/Subsetting.html

Sin embargo, éste método es bastante lento en mi equipo (Dell Optiplex 9010) y con mis datos (unas 90.000 filas en el fichero de entrada). Pongamos, por referencia y para comparar con métodos que veremos a continuación, que éste ejemplo en el que usamos una lista y un bucle 'for' tarda 1' (54'' en promedio, realmente).

Parte de por qué es así puede tener que ver con el bucle 'for'. Más adelante veremos alguna alternativa al bucle. Sin embargo, también tiene que ver el rendimiento de la búsqueda de los elementos en la lista.

Una alternativa al operador '%in%' podría ser:

 if (exists(identifier, gene2GO)){  

Pero en principio esta opción es más lenta, en mi caso unos 2' (1'57'').

En ambos casos, la comprobación la hacemos para acceder después al gen en cuestión. Por tanto, se está buscando 2 veces el mismo elemento en la lista: una en la comprobación, otra al obtener el elemento. Una forma de evitar esto sería cogiendo directamente el objeto (1 sola búsqueda) y comprobando si existe o no (NULL):

 current = gene2GO[[identifier]]  
 if (is.null(current)){  

(Nota: ojo que ahora el if comprueba si NO existe, no si existe como en los ejemplos anteriores).

Vemos que con esta otra forma tenemos un tiempo promedio de unos 26''. Aproximadamente la mitad que con el uso de '%in%', como esperaríamos.

Sin embargo, ya hemos comentado antes que el uso de un bucle 'for' en R también puede tener algo que ver con la velocidad de nuestro código (por ejemplo ver: http://faculty.nps.edu/sebuttre/home/R/apply.html). En R se suele preferir utilizar funciones que van recorriendo las estructuras de datos y aplicando una función sobre cada elemento de la estructura. Es el caso de la familia de funciones 'apply' (ver por ejemplo https://nsaunders.wordpress.com/2010/08/20/a-brief-introduction-to-apply-in-r/). Veamos el equivalente al código anterior usando la función 'apply':

 gene2GO = list()  
 add_to_list <- function(go_tab){   
  identifier = go_tab[[1]]  
  go_term = go_tab[[2]]  
  current = gene2GO[[identifier]]  
  if (is.null(current)){  
   gene2GO[[identifier]] <<- as.character(go_term)  
  } else {  
   gene2GO[[identifier]] <<- c(gene2GO[[identifier]], as.character(go_term))  
  }  
 }  
 invisible(apply(go_tab, 1, add_to_list))  

Ahora está tardando unos 21'' en mi máquina. Una ligera mejora respecto al uso del bucle 'for', que podría tener impacto en conjuntos de datos muy grandes o cuando queremos correr muchos análisis.

Vemos varias peculiaridades en el código anterior. La función 'apply' va a recorrer nuestra lista y llamará a la función 'add_to_list' para cada elemento. Sin embargo, en nuestro caso no queremos simplemente transformar los elementos de la lista (por ejemplo, formateando el nombre del gen), si no que queremos poblar una estructura de datos diferente a partir de los valores de la lista recorrida. Por eso se utiliza el operador <<-:

"The operators are normally only used in functions, and cause a search to made through parent environments for an existing definition of the variable being assigned."

Otra posible opción sería utilizar un parámetro opcional en la función, pero no me queda claro que vaya a funcionar para escribir en él, ya que el parámetro (la lista de vectores que estamos creando) necesitaría pasarse por referencia, y me suena que en R no he pasado parámetros así anteriormente ¿Os suena? ¿Algún consejo sobre esto? ¿Alguna otra opción?

La otra peculiaridad es la función 'invisible'. Ya que las funciones en R devuelven un valor y 'apply' lo propaga para obtener la posible transformación de 'go_tab' en otra cosa, si no asignamos el resultado de 'apply' a una variable tendremos el resultado en nuestra salida. Para evitar esto, podemos asignar el resultado a una variable "dummy" o utilizar el método 'invisible', que a mi me ha parecido más adecuado y legible.

Otra diferencia, más básica, es cómo accedemos a nuestros datos originales. En el primer ejemplo, lo que hacíamos es acceder a los campos del 'data.frame' y, en un campo dado, obtener el dato correspondiente según el bucle for:

identifier = go_tab$identifier[i]

Con 'apply' en cambio estamos recibiendo cada una de las líneas del 'data.frame' como un 'vector', por lo que accedemos directamente a los datos en éste vector:

identifier = go_tab[[1]]

Por ahora hemos mejorado nuestro algoritmo cambiando la forma de comprobar la existencia de un elemento en la lista para luego obtenerlo, y pasando de un bucle 'for' a usar la función 'apply'. Pero aún podemos ir más allá. ¿Por qué utilizar una lista? Quizás podamos mejorar el rendimiento utilizando un entorno con un mapeado "hash" (https://stat.ethz.ch/R-manual/R-devel/library/base/html/environment.html).

 gene2GO = new.env()  
 add_to_list <- function(go_tab){  
  identifier = go_tab[[1]]  
  go_term = go_tab[[2]]  
  current = gene2GO[[identifier]]  
  if (is.null(current)){  
   gene2GO[[identifier]] <<- as.character(go_term)  
  } else {  
   gene2GO[[identifier]] <<- c(gene2GO[[identifier]], as.character(go_term))  
  }  
 }  
 invisible(apply(go_tab, 1, add_to_list))  
 gene2GO = as.list(gene2GO)  

Aquí, en lugar de una lista creamos un entorno ('environment'), que por defecto utiliza una función "hash" para mapear su contenido. Por lo demás, como vemos, lo estamos manejando aquí de forma muy similar a la lista, aunque hemos añadido al final una conversión a formato de lista (en una sola línea pasamos de nuestro "hash" a una lista de vectores).

Éste código tarda unos 7'', incluída la conversión. Hay que tener en cuenta que el uso de un "hash" puede ser más o menos aconsejable según el volumen de datos con el que se trabaje. Desde luego, con nuestros datos, hemos mejorado bastante el rendimiento prestando atención a distintas alternativas que ofrece R y cuál sería la mejor combinación de las que hemos probado.

Y ya tendríamos nuestra lista creada. Nuestra intención era utilizar esta lista para informar a topGO de qué términos de GOs se relacionan con nuestros genes ¿recordáis? Le pasaremos esta lista a topGO como parámetro 'gene2GO' en la función en la creamos el primer objeto necesario con éste paquete, e indicaremos que la función de anotación es del tipo "annFun.gene2GO". Para eso, y terminar así la preparación de los datos, tendríamos que importar las listas de genes a comparar, o una sola lista con todos pero con algún criterio que permita diferenciar un subconjunto de otro. Crearíamos entonces el objeto con:

 GOdata <- new("topGOdata", description="whatever", ontology = "BP", allGenes = background, geneSel = test,   
          annot = annFUN.gene2GO, gene2GO = gene2GO, nodeSize = 10)  

Donde:
'background' es nuestro conjunto de genes de referencia.
'test' es el grupo de genes que estamos interesados en comparar con la referencia.
'gene2GO' la lista que hemos generado.

Seguiríamos entonces con las siguientes fases de topGO: los test de enriquecimiento y el análisis de los resultados. Si te ha entrado el gusanillo de probar, topGO está perfectamente documentado:

http://www.bioconductor.org/packages/release/bioc/vignettes/topGO/inst/doc/topGO.pdf

Un saludo!


20 de abril de 2015

curso de verano "Estructura y Función de Proteínas"


Cuándo y dónde: Del 6 al 10 de Julio en el Palacio de Congresos de Jaca, Huesca.

Información completa: https://cursosextraordinarios.unizar.es/curso/2015/estructura-y-funcion-de-proteinas-v-edicion


El curso está dirigido a:
  • Estudiantes de los últimos cursos de los grados en Química, Biología, Física, Farmacia, Medicina, Bioquímica, Biotecnología y Veterinaria.
  • Estudiantes de másteres de Bioquímica, Biología Molecular y Celular, Biotecnología y Tecnología de Alimentos.
  • Estudiantes de doctorado.
Profesionales en activo en cualquiera de las áreas de conocimiento y especialidades arriba indicadas, o de aquellas afines con interés por la Biofísica.

El curso tiene como objetivo mostrar y evaluar diferentes metodologías de uso habitual en los laboratorios de bioquímica y biofísica de proteínas con objeto de mejorar el conocimiento sobre la relación estructura y función en estas macromoléculas, sin olvidar sus posibles aplicaciones para mejorar nuestra sociedad. Como ponentes participan en el curso Profesores e Investigadores especialistas en cada una de las áreas de conocimiento.

Como novedad, respecto a ediciones anteriores, se propone el último día un taller de “Cómo escribir un artículo científico” y “de los entresijos, a nivel editorial, detrás de algunas revistas del grupo Elsevier” (uno de los gigantes en la actualidad en la publicación de textos científicos). Este taller es de interés para cualquier profesional en investigación y docencia.

Organizado por: Milagros Medina y José Luis Neira


14 de abril de 2015

HOWTO install Grid Engine on multi-core Linux box to run GET_HOMOLOGUES

Hi,
I post this in English so that all users of GET_HOMOLOGUES can benefit. It is inspired on a previous tutorial posted in scidom and the expertise of David Ramírez from our provider SIE.


The aim of this entry is to demonstrate how to set up a Grid Engine queue manager on your local multi-core Linux box or server so that you can get the most of that computing power during your GET_HOMOLOGUES jobs. We will install Grid Engine on its default path, /opt/, which requires superuser permissions. As I'll do this un Ubuntu, I will do 'sudo su' to temporarily get superuser privileges; this should be replaces simply with 'su' in other Linux flavours. Otherwise you can install it elsewhere if you don't have admin rights.

1) Visit http://gridscheduler.sourceforge.net , create a new user and download the latest 64bit binary to /opt:

$ cd /opt
$ sudo su 
$ useradd sgeadmin
$ wget -c http://dl.dropbox.com/u/47200624/respin/ge2011.11.tar.gz $ tar xvfz ge2011.11.tar.gz
$ chown -R sgeadmin ge2011.11/
$ chgrp -R sgeadmin ge2011.11/

$ ln -s ge2011.11 sge

2) Set relevant environment variables in /etc/bash.bashrc  [system-wide, can also be named /etc/basrhc] or alternatively in ~/.bashrc for a given user:

export arch=x86_64
export SGE_ROOT=/opt/sge
export PATH=$PATH:"$SGE_ROOT/bin/linux-x64"
export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:$SGE_ROOT/lib"

And make this changes live:

$ source /etc/bash.bashrc

3) Set your host name to anything but localhost by editing /etc/hosts so that the first line is something like this (localhost or 127.0.x.x IP addresses are not valid):

172.1.1.1   yourhost

4) Install Grid Engine server with all defaults except cluster name, which usually will be 'yourhost':
$ ./install_qmaster

5) Install Grid Engine client with all defaults:
$ ./install_execd

6) Optionally configure default all.q queue:
$  qconf -mq all.q

7) Add your host to list of admitted hosts:
$ qconf -as yourhost

$ exit

You should now be done. A test will confirm this. Please open a new terminal and move to your GET_HOMOLOGUES installation folder:

$ cd get_homologues-x86-20150306
$  ./get_homologues.pl -d sample_buch_fasta -m cluster

If the jobs finishes successfully then you are indeed done. I hope this helps,
Bruno





23 de marzo de 2015

Curso de Python para biólogos - Lección 5. Funciones y ejercicios prácticos

Hoy descansaremos un poco y sólo explicaremos un nuevo concepto teórico, pero muy importante, las funciones en Python. Tras ello se propondrán varios ejercicios para aplicar los conocimientos adquiridos en las lecciones previas (1, 2, 3 y 4).

Funciones en Python

Una función es un bloque de código que toma ciertos datos o variables como argumentos y devuelve otros datos o variables como resultado. Aunque no nos hayamos dado cuenta, ya hemos usado funciones anteriormente: print(), input(), int(), str(), len(), type(), min(), max()...dichas funciones vienen por defecto con Python y alguien las ha escrito previamente para que nosotros las podamos usar. La importancia fundamental de las funciones es que nos permiten escribir un código una sola vez y utilizarlo todas las veces que sea necesario sin volver a escribirlo, lo cual nos ahorra tiempo y simplifica los programas.

Las funciones se especifican con 'def' seguido del nombre de la función y entre paréntesis tantas variables como argumentos utiliza, estas variables guardarán los argumentos que demos al llamar a la función. No hay que olvidar los dos puntos al final de esta primera línea y escribir el bloque de código de la función tabulado (4 espacios). Al final de la función el resultado de la misma se devuelve con el comando 'return'. Para ejecutar o llamar una función simplemente hay que escribir su nombre y entre paréntesis sus argumentos. El resultado/s de las funciones los podemos guardar en variables o usarlos directamente.

Veamos en el primer ejemplo cómo escribir una función que realiza la suma de dos números y la devuelve como resultado. En el segundo ejemplo, la función recibe un nombre y devuelve una frase dando la bienvenida.

Ejercicio 1: Calcular el factorial de un número

El factorial de un número se define como el producto de todos los números enteros positivos desde 1. El ejercicio propuesto consiste en escribir una función que lo calcule usando un bucle 'while' que vaya restando (o sumando) número a número y guardando el resultado de la multiplicación por el anterior. No vale mirar en Google y usar funciones ya escritas.

Ejercicio 2: Contar el número de nucleótidos en una secuencia de DNA

Ahora escribiremos una función que cuente el número de pares de bases de una secuencia de DNA tomada como argumento de la función. Creo que este problema será más sencillo que el anterior, el código de la función puede ocupar tan sólo una línea, pero sugiero utilizar el bucle 'while' para leer letra a letra y aprender más.

Ejercicio 3: Encontrar un codón de inicio de la traducción en una secuencia de DNA

El ejercicio consiste en encontrar el codón de inicio de la traducción, osea la palabra 'ATG', en una secuencia de DNA. La función puede ser muy sencilla utilizando expresiones regulares, pero como no las hemos estudiado todavía, propongo reutilizar el código del Ejercicio 2 que debería leer letra a letra la secuencia e ir extrayendo palabras de 3 letras en cada iteración del bucle. Comprobar si la palabra es 'ATG' con una sentencia condicional, si lo es retornar la posición como resultado de la función, y si no seguir hasta terminar de leer toda la secuencia.

Ejercicio 4: Calculando una secuencia consenso

El último ejercicio de hoy consiste en escribir una función que tome como argumentos 3 secuencias de DNA (preferiblemente de la misma longitud, para evitar problemas y errores, si no calcular el consenso hasta la longitud de la secuencia más corta) y devuelva una única secuencia que contenga en cada posición el nucleótido más frecuente o en su defecto la letra 'N'. Este ejercicio es el más complejo y requiere utilizar bucles 'while', sentencias condicionales...


Buena suerte y nos vemos en la próxima lección.

19 de marzo de 2015

Curso de Python para biólogos - Lección 4. Entrada de datos y bucles while

En esta nueva lección aprenderemos a pedir datos al usuario desde nuestro programa y a realizar nuestros primeros bucles. Como veremos, un bucle consiste en repetir un bloque de código varias veces sin tener que volver a escribirlo.

Entrada de datos desde la línea de comandos

A veces necesitaremos introducir datos en nuestro programa de forma interactiva, lo que significa que preguntaremos o pediremos al usuario cierta información (por ej. linux siempre nos pregunta antes de instalar un nuevo programa y debemos responder 'yes' o 'no'). En Python se utiliza la función 'input', dicha función devuelve el texto introducido por el usuario antes de pulsar Enter. El texto tendrá formato de cadena (string), por lo que no podrá ser utilizado directamente en caso de ser un número, habrá que convertir la cadena en número.

Conversión entre diferentes formatos de datos

Muchas funciones y métodos de Python sólo aceptan un tipo de datos, por ejemplo 'print' requiere cadenas o las operaciones aritméticas requieren números. Los tipos de datos que conocemos son: números, cadenas, listas y diccionarios. Simplificando, en Python existen 2 clases de números, los enteros (int) y los decimales (float). Por ello es importante conocer cómo convertir una cadena (devuleta por la función 'input' por ejemplo) en un número antes de realizar operaciones aritméticas, o convertir un número en una cadena antes de imprimirlo.
Si intentamos elevar al cuadrado un número introducido mediante 'input' nos dará error, porque antes deberemos convertir la cadena devuelta por 'input' en un tipo de número (int o float).
Si el texto devuelto por 'input' contiene decimales, al intentar convertirlo en número entero nos devolverá un error. La solución es convertirlo a decimal primero. Pero cuidado, operar con números enteros es diferente que con números decimales, y los resultados pueden variar mucho.

Validación de datos con 'try' y 'except':

Como hemos visto en los ejemplos anteriores, muchas veces tenemos errores al intentar operar con los diferentes tipos de datos. Una forma de evitar problemas es detectar los errores con la sentencia 'try:' y evitar que se ejecute el código que origina el error para ejecutar un código alternativo especificado por 'except:'. Se podría traducir cómo: "en caso de error... hacer lo siguiente...".
 Notar cómo siempre hay diferentes maneras de hacer lo mismo...
Pero 'try' y 'except' no son milagrosos, no permiten por ejemplo corregir errores en la sintaxis de nuestro código, como puede ser olvidarnos un paréntesis o unas comillas.

Bucles 'while':

Por fin vamos a adentrarnos en el fascinante mundo de los bucles. Un bucle es un bloque de código que es ejecutado varias veces, pero sólo hay que escribirlo una vez :)

La particularidad de un bucle 'while' es que es ejecutado indefinidamente (bucle infinito) mientras la condición que le sigue es cierta (True). Como en el caso de las sentencias condicionales con 'if', el código a ejecutar mientras se cumpla la condición debe ser tabulado 4 espacios para ser reconocido por Python. Veamos unos ejemplos.
Empecemos a aplicar los conocimientos adquiridos y veamos cómo usar sentencias condicionales (if...elif...else...) dentro de un bucle, la tabulación será muy importante...

Ejercicio1. Programa para adivinar un número:

Para terminar la lección haremos dos ejercicios muy interesantes. El primero será escribir un código que genera un número aleatorio entre 1 y 100, y nos preguntará sucesivante por el número hasta que lo adivinemos. Las primeras líneas (from random import randint y number=randint(1,100)) son necesarias para generar el número aleatorio guardado en la variable 'number' pero serán explicadas en próximas lecciones.
Pero, qué pasa si nos equivocamos y no escribimos un número? El programa devolverá un error a no ser que nos preparemos para ello:

Ejercicio2. Cantar una canción

En el último ejercicio de hoy, crearemos un código que nos pida escribir una canción línea a línea y después la mostrará en pantalla 2 veces, imprimiendo línea por línea con un bucle 'while'.
  
Próxima lección: bucles 'for'