31 de marzo de 2013

Procesar argumentos con múltiples opciones en un script en Python

En los últimos tiempos estoy vendiendo mi alma al diablo y he empezado a usar Python. Personalmente prefiero Perl, Python me parece un lenguaje más anárquico y desorganizado. Sin embargo, Python está de moda en Bioinformática y muchos de los módulos más actualizados se están escribiendo en este lenguaje.

Hoy voy a hablar del problema que me he encontrado para encontrar un módulo que procese las opciones de línea de comandos que normalmente introducimos en los scripts, por ej.: ./script.py -i archivo1 -o archivo2 Estas opciones suelen ser archivos de entrada, salida, o simplemente parámetros que cambian la ejecución del programa. También sirve para pedir ayuda sobre como usar el script: ./script.py -h 

En Perl existe el módulo Getopt::Long para procesar dichos parámetros de línea de comando de una forma sencilla y robusta.

En Python el módulo equivalente es getopt. A continuación se muestra un ejemplo de uso. Sin embargo este módulo no contempla la posibilidad de procesar múltiples opciones para un mismo argumento, por ejemplo, varios archivos de entrada.

 #!/usr/bin/python  
 # -*- coding: utf-8 -*-   
 """  
      Script que usa el módulo 'getopt' para parsear los argumentos y opciones de línea de comandos  
 """  
 import sys, re, os, copy  
 # Importar el módulo 'getopt'  
 import getopt  
 # Imprimir en pantalla el script y sus argumentos y opciones  
 print 'ARGV:', " ".join(sys.argv[0:]),"\n"  
 # Información de ayuda de uso del script  
 def help() :  
      """Ayuda sobre opciones del script en línea de comandos"""  
      print "usage:",sys.argv[0], "[options]\n"  
      print " -h this message\n"  
      print " -i input file/s\n"  
      print " -o output file\n"  
 # Definir los argumentos obligatorios (:) y opcionales, en sus formas abreviada y completa  
 try:  
      options, args = getopt.getopt(sys.argv[1:], "hi:o:", ["help", "input", "output"])  
 except getopt.GetoptError as err:  
      # Imprimir ayuda y salir si hay algún error o argumento incorrecto  
      print str(err)  
      help()  
      sys.exit(2)  
 output = None  
 verbose = False  
 # Definir el tipo de las variables que guardarán las opciones  
 INP_input = ''  
 INP_output = ''  
 # Asignar las opciones a cada argumento  
 for _opt, _arg in options:  
      if _opt in ("-i", "--input"):  
           INP_input = _arg  
      elif _opt in ("-o", "--output"):  
           INP_output = _arg  
      elif _opt in ("-h", "--help"):  
           help()  
           sys.exit()  
      else:  
           assert False, "unhandled option"  
 # Imprimir en pantalla todo lo anotado  
 print "Script:",sys.argv[0],"\n"  
 print "Argumentos y opciones:\n",  
 for _opt,_arg in options :  
      print " -"+_opt+": "+_arg  

Para solucionar este problema de Python de procesar argumentos de línea de comandos con múltiples opciones he creado mi propio código, puesto que me parece fundamental para cualquier script bioinformático. En el siguiente ejemplo se pueden definir listas para guardar múltiples opciones para un único argumento (ej. varios archivos). El código es provisional y le faltarían muchas opciones que sí tiene el módulo getopt, sin embargo, ofrece una solución al problema y espero ir mejorándolo con el tiempo.

 #!/usr/bin/python  
 # -*- coding: utf-8 -*-   
 """  
      Script que parsea los argumentos y opciones de línea de comandos  
      permitiendo usar listas para guardar argumentos con múltiples opciones  
 """  
 import sys, re, os, copy  
 from types import *  
 # Imprimir en pantalla el script y sus argumentos y opciones  
 print 'ARGV:', " ".join(sys.argv[0:]),"\n"  
 # Información de ayuda de uso del script  
 def help() :  
      """Ayuda sobre opciones del script en línea de comandos"""  
      print "usage:",sys.argv[0], "[options]\n"  
      print " -h this message\n"  
      print " -i input file/s\n"  
      print " -o output file\n"  
      exit()  
 # Definir el tipo de las variables que guardarán las opciones  
 INP_input = []  
 INP_output = ''  
 # Asignar a cada argumento una de las anteriores variables  
 options = {'h': 'help', 'help': 'help', 'i': 'INP_input', 'input': 'INP_input', 'o': 'INP_output', 'output': 'INP_output'}  
 argv_var = str()  
 # Leer los argumentos y opciones  
 for _argv in sys.argv[1:] :  
      # Leer los argumentos que comienzan con guión  
      if re.match("^-", _argv) :  
           for _opt,_var in options.iteritems() :  
                if re.compile("^-%s$" % (_opt)).match(_argv) :  
                     argv_var = copy.deepcopy(_var)  
                     if type(vars()[argv_var]) is FunctionType :  
                          vars()[argv_var]()  
                     break  
      # Anotar las opciones de cada argumento  
      else :  
           if argv_var is not None :  
                if type(vars()[argv_var]) is ListType :  
                     vars()[argv_var].append(_argv)  
                else:  
                     vars()[argv_var] += _argv  
 # Imprimir en pantalla todo lo anotado  
 print "Script:",sys.argv[0],"\n"  
 print "Argumentos y opciones:\n",  
 for _opt,_var in options.iteritems() :  
      if vars()[_var] is not None and type(vars()[_var]) is not FunctionType:  
           #print _opt,pprint(type(vars()[_var]))  
           if type(vars()[_var]) is ListType:  
                print " -"+_opt+": "+str(", ".join(vars()[_var]))  
           else :  
                print " -"+_opt+": "+str(vars()[_var])  

8 de marzo de 2013

Beca de verano en regulación del cloroplasto

Hola,
nuestro laboratorio en la Estación Experimental de Aula Dei (EEAD-CSIC) tiene experiencia en el desarrollo de herramientas bioinformáticas para el estudio de las redes de regulación transcripcional [1,2] y ha contribuído recientemente al análisis de los genomas cloroplásticos [3].  Ahora ofertamos una beca de verano, de la convocatoria JAE-IntroCP,  con dos objetivos:

1) Sentar las bases para un análisis evolutivo de las huellas filogenéticas de
regulación transcripcional en el cloroplasto, teniendo en cuenta que hay ya
disponibles más de 200 genomas plastídicos en las bases de datos públicas. Este
estudio pretende identificar secuencias reguladoras en los promotores de genes
cloroplásticos ya sean codificados en el propio cloroplasto o en el genoma
nuclear, con el fin de ampliar nuestro conocimiento de la regulación de este
orgánulo esencial para la viabilidad de las plantas [4,5].

2) El aprendizaje por parte del estudiante de herramientas bioinformáticas para el análisis funcional de genomas. A lo largo de la estancia la persona seleccionada se familiarizará con el trabajo cotidiano de la investigación en Bioinformática en el campo de la genómica. Esperamos que los alumnos interesados en este proyecto tengan sólidos conocimientos en Biología Molecular y/o Bioquímica y se valorará especialmente experiencia previa o interés en programar y aplicar herramientas bioinformáticas.

Si estás interesado por favor contacta con:
Inmaculada Yruela Guerrero ( i.yruela at eead.csic.es  )
Bruno Contreras Moreira ( bcontreras at eead.csic.es )

Actualización: Convocatoria publicada el 8 de Abril, el plazo acaba el 23 de Abril


Ejemplos de recursos y publicaciones relevantes para este proyecto son:

[1] http://floresta.eead.csic.es/footprintdb

[2] Lozada-Chávez I, Espinosa Angarica V, Collado-Vides J. and Contreras-Moreira
B.(2008) The role of DNA-binding specificity in the evolution of bacterial regulatory networks. Journal of Molecular Biology, 379(3): 627-43.
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2491489/

[3] Yruela,I. and Contreras-Moreira,B. (2012) Protein disorder in plants: a view from the chloroplast. BMC Plant Biology, 12:165
http://www.biomedcentral.com/1471-2229/12/165

[4] Liere K, Weihe A and Börner T. The transcription machineries of plant
mitochondria and chloroplasts: Composition, function, and regulation. Journal of Plant Physiology 168 (2011) 1345–1360.
http://www.ncbi.nlm.nih.gov/pubmed/21316793

[5] Wang Y, Ding J, Daniell H, Hu H and Li X. Motif analysis unveils the possible co-
regulation of chloroplast genes and nuclear genes encoding chloroplast proteins.Plant Mol Biol. 2012 Sep;80(2):177-87.
http://link.springer.com/article/10.1007%2Fs11103-012-9938-6



4 de marzo de 2013

PubReader, el guiño de PubMed a los tablets con HTML5 y CSS3

PubMed, la mejor base de datos de artículos científicos online, en concreto su subsección de consulta gratuita PubMed Central (PMC), ha estrenado el año con novedades tecnológicas. Una de las cosas que más me gusta de PubMed y en general de las bases de datos del NCBI, además de que su consulta es gratuita, es su sobriedad, su buen funcionamiento y su resistencia a añadir florituras a su interfaz online. Creo que la máxima de 'si algo funciona bien no hay que intentar cambiarlo' se aplica a PubMed.

Sin embargo, el otro día al consultar nuestro nuevo artículo en PubMed me sorprendí al ver a la derecha una imagen de un iPad, similar a la típica que muestra Amazon con su Kindle. La imagen llevaba el título 'Click here to read article using PubReader' y por supuesto probé qué era eso. Muy gratamente descubrí que se trata de PubReader, una nueva herramienta para leer cómodamente los artículos en pantallas táctiles y tablets.


PubReader está desarrollado con HTML5 y CSS3. El esquema de visualización de PubReader es muy simple, lee el artículo de PMC en formato XML y lo transforma en HTML que junto con algo de Javascript y CSS simula un lector de libros al modo de Kindle o del nuevo Reader de Windows 8. Además podemos descargarnos el código fuente de PubReader para hacer nuestro propio lector online de libros.

22 de febrero de 2013

Probando deltablast

Hola,
hace tiempo que tenía pendiente escribir sobre una de las ramas más recientes de la familia de programas BLAST, que se llama deltablast, publicada el año pasado. La lectura del artículo original sugiere que este programa se desarrolló a partir de la publicación del algoritmo CS-BLAST (context specific BLAST).
Pero realmente la historia no comienza aquí.
Todo empezó con el algoritmo PSI-BLAST, una versión iterativa de BLASTP, que en vez comparar una secuencia de proteína contra una librería de secuencias, construye un perfil de la secuencia problema y ése es el que compara contra la librería, ganando en sensibilidad. Todo esto con un coste en tiempo de cálculo modesto.
Qué problema tiene esto? Pues que el perfil se construye en tiempo real con los homólogos encontrados en iteraciones previas de manera automática y en ocasiones se pueden contaminar y producir resultados no idóneos.

(tomada de http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2645910)
CS-BLAST, que en realidad es como un preprocesamiento de BLAST, logró superar estos problemas precalculando los perfiles, que capturan información del contexto, y de hecho mejoró la sensibilidad a la hora de detectar homólogos remotos, que han perdido claramente la similitud de secuencia.

Pues bien, deltablast es un programa del NCBI que se inspira en estas experiencias para superar a capacidad de búsqueda de CS-BLAST, apoyándose en la colección de dominios Conserved Domains Database (CDD):
(tomada de http://www.biology-direct.com/content/7/1/12)
El programa es muy sencillo de usar en la web, aquí os explico como probarlo localmente en vuestra máquina:
1) descarga la última versión de BLAST+ para tu arquitectura de
ftp://ftp.ncbi.nih.gov/blast/executables/LATEST
 
2) descomprime el software en una carpeta, por ejemplo soft/

3) descarga una copia de CDD de ftp://ftp.ncbi.nih.gov/blast/db/cdd_delta.tar.gz

4) descomprime la base de dominios en un carpeta, de nuevo por ejemplo en  soft/

5) elige una colección de secuencias protéicas contra la que buscar, por ejemplo soft/nr.faa y formatéala con:
 $ soft/ncbi-blast-2.2.27+/bin/makeblastdb nr.faa
 
6) Ya puedes buscar una secuencia de proteína, como ejemplo.faa, contra tu colección de secuencias:
$ soft/ncbi-blast-2.2.27+/bin/deltablast -query ejemplo.faa -db soft/nr.faa -rpsdb soft/cdd_delta

7) Puedes utilizar prácticamente los mismos parámetros que con BLASTP, que puedes consultar haciendo:
$ soft/ncbi-blast-2.2.27+/bin/deltablast -help

Suerte!
Bruno

8 de febrero de 2013

Filogenias universales 16S

Un obituario en el último número de Science (8 Feb 2013) me ha hecho recordar el tremendo impacto que ha tenido el trabajo del desaparecido Carl Woese en nuestra visión de la historia evolutiva, ya que a partir de sus logros pudimos comparar por primera vez con precisión molecular las historias evolutivas de linajes tan distantes como las bacterias y los animales. Además de descubrir por el camino nada menos que a las arqueobacterias, el uso de las secuencias de los genes ribosomales 16S como fuente de información filogenética todavía es rutinario, a pesar de sus limitaciones y de los avances tecnológicos, y a pesar de que ahora se hable más de redes que de árboles. Si no que se lo pregunten a los que se dedican a la metagenómica...

En las siguientes figuras se resume el impacto de la introducción de los análisis de genes 16S:

Árbol filogenético de Haeckel (tomado de http://mmbr.asm.org/content/51/2/221.long).
Árbol de distancias basado en el alineamiento de 260 secuencias de rRNA 16S (tomado de http://mmbr.asm.org/content/51/2/221.long).

Al releer la revisión de Woese del año 1987 lo que más me ha sorprendido ha sido ver que no sólo se usó secuencias, sino que se molestó en analizar en detalle la estructura de los RNAs ribosomales, e incluso localizó a nivel de estructura secundaria las zonas que permitían distinguir los grandes linajes evolutivos, como se muestra en la siguiente figura. Está ya todo inventado también aquí?

(tomada de http://mmbr.asm.org/content/51/2/221.long)

Un saludo y buen finde,
Bruno