26 de enero de 2012

Apuntes sobre las XI Jornadas Bioinformáticas

Hola,
tras el anuncio de hace unas semanas, voy a comentar un poco mis impresiones sobre las XI Jornadas de Bioinformática. Después de la sesión inaugural del lunes, a su vez precedida por el simposio de estudiantes, a las que no pude asistir, el día 24 realmente empezaron mis jornadas. Como valoración general, creo que los asistentes hemos tenido ocasión de aprender y de discutir sobre los problemas del campo, y hemos tenido la oportunidad de escuchar charlas muy buenas. Para los que no pudieron venir, ahí va mi resumen de las charlas a las que asistí y notas sobre algunos pósters.
Modelo 3D de la region cromosómica ENCODE ENm008, que contiene el locus de la α-globina, adaptado de www.ncbi.nlm.nih.gov/pmc/articles/PMC3056208.




24 de Enero
Hoy hemos podido escuchar un montón de charlas sobre temas variopintos, y me ha llamado la atención la abrumadura presencia de la palabra mágica NGS (Next Generation Sequencing) en muchas de ellas. También en los pósters expuestos ha habido muchos ejemplos de la aplicación de estas herramientas, sobre todo de RNAseq.

De lo que he visto muy poco hoy ha sido de ChIPseq, tan sólo el póster de Ionas Erb con el software Pro-Coffee para alinear secuencias de DNA de promotores, publicado en NAR.

Otro póster que me llamó la atención fue el trabajo de Minoche sobre la evaluación sistemática de los errores típicos de la plataforma de secuenciación Solexa, publicada en Genome Biology.

Sonia Tarazona me explicó su póster sobre RNAseq y me invitó a probar el  software Qualimap, que permite evaluar el efecto del coverage sobre la interpretación de las diferencias de expresión medidas en RNAseq.

Leo Mirny nos explicó en detalle la técnica de Hi-C para localizar regiones de cromatina cercanas en el espacio celular y cómo en su grupo han usado este tipo de datos para entender el empaquetamiento del núcleo de levaduras y de Homo sapiens, construyendo una matriz que se parece mucho a un mapa de contactos de proteínas, algo que también por la tarde explicó Davide Bau en su trabajo sobre la transcripción y la estructura de la cromatina en un locus de alfa globina.

Juan Ramón González nos comentó los métodos predominantes actuales para la normalización de conteos de lecturas RNAseq (TMM, EDAseq y CQN) y presentó los problemas que tienen las distribuciones de Poisson (con un parámetro libre) y la binomial negativa (con dos) para modelar algunos datos reales, y mostró ejemplos convincentes del uso de la de Poisson-Tweedie como lo mejor de los dos mundos, con un tercer parámetro para elegir según el caso el mejor modelo estadístico, dada la dispersión de los datos reales de esta teconología. Propone su paquete de Bioconductor/R tweeDEseq como herramienta para esta tarea.

Eva María Novoa nos dió una clase magistral sobre el uso de codones en procariotas y eucariotas, presentando evidencia de la importancia de las enzimas UMS (en proka) y hetADATS (adenosine deaminasas de euka), que modifican terceras bases de los tRNAs, para explicar las diferentes frecuencias de codones en todos los bichos conocidos. Su trabajo se publica en Cell.

Nacho Medina nos deslumbró con la capacidad de su equipo del CIPF para crear una tubería de análisis de datos de NGS que nos permitirá como usuarios hacerlo todo en sus servidores en "tiempos de minutos", aprovechando la optimización que han hecho de los distintos algoritmos y del hardware subyacente, que incluye, si no recuerdo mal, CPUs, GPUs dedicadas y discos de estado sólido. Me llamó mucho la atención su navegador genómico HTML5, que se puede probar en http://genomemaps.org.

Tomas Marques nos volvió a hablar de primates, ya lo había hecho en Málaga en el 2010. Esta vez trató de convencernos de la importancia de mirar con lupa los datos de NGS, en su caso de ensamblaje genómico, de usar el software adecuado para nuestros objetivos, y de hacer el control de calidad en casa, sin delegarlo. Menciona que en su labo usan GATK como software para variant calling tras haberlo comparado con otros.

Tanya Vavouri nos explicó, si lo entendí bien, que los espermatozoides humanos maduros conservan sólo un 4% de los nucleosomas, pero que justo esos pueden ser muy importantes para pasar información epigenéticas al nuevo cigoto, porque se correlacionan con picos de %GC muy cercanos a promotores.

Javier Macia nos mostró ejemplos de cómo modelar circuitos electrónicos a base de puertas lógicas implementadas con células de levadura modificadas.

Ya hacia el final del día Toni Giorgino nos mostró dinámicas moleculares espectaculares de un dominio SH2 y su ligando, y Pablo Minguez publicó un montón de resultados sobre la conservación de sitios en proteínas eucariotas que pueden sufrir modificaciones postraduccionales. Lo siento, a las dos últimas charlas no me pude quedar.

25 de Enero
El tercer y último día del congreso arrancó con una conferencia plenaria de Luis Serrano donde nos resumió los estudios de su grupo sobre Mycoplasma pneumoniae, un parásito bacteriano con genoma extremadamente reducido que se puede cultivar en el laboratorio. Su charla se puede resumir  como la aplicación de todas las herramientas bioinformática, genómicas, proteómicas y metabolómicas disponibles para tratar de caracterizar la biología de este bicho, que tiene solamente unos 10 factores de transcripción homólogos de otros conocidos en otras especies. Por destacar algo de una charla muy densa pero amena, Luis habló que encontraban una correlación <0.50 entre los niveles de expresión génica y las cantidades de proteína detectadas por MS en distintas condiciones. Esta observación no es muy novedosa, pero sí la explicación que proponía, basada en la divergencia de las secuencias Shine-Dalgarno en los mensajeros, que motivarían menores afinidades de los ribosomas y por tanto menores tasas de traducción.

Antonio Mérida presentó el software Sma3s para la anotación de genomas, que comparó con otros programas como Blast2GO. Roderic Guigó apuntaba  tras la charla que la anotación de un genoma actualmente debe incluir no sólo genes codificantes sino también RNAs reguladores, por ejemplo.

Paolo Ribeca dió la primera charla del día dedicada al tema estrella (NGS). En esta conferencia Paolo repasó los principales problemas del mapeo de lecturas (reads) sobre un genoma de referencia y las limitaciones de los principales programas (Bowtie,BWA,SOAP,MrFAST,MrsFAST) a la hora de hacer búsquedas exhaustivas y flexibles sobre posibles dianas genómicas, algo que su propia plataforma GEM parece haber resuelto y acelerado considerablemente. Su mensaje de precaución es que el usuario de este tipo de software debe conocer con cierta precisión cómo funciona el programa que va a usar y qué limitaciones tiene, en vez de confiar en el programa a ciegas y dejar que tome decisiones, no siempre transparentes, por ti.

Darío Guerrero mostró datos sobre la validación de una tubería de preprocesamiento de lecturas NGS y el ensamblaje de datos de RNAseq, con seqtrimnext y fulllengther, respectivamente. Uno de los programas con los que comparó sus resultados fue Mira.

Beatriz García nos contó el desarrollo de software de aprendizaje automático para la asignación de secuencia de proteínas sin anotar a rutas metabólicas.

Ya en la tarde, Patrick Aloy nos contó un proyecto de su grupo que reconstruye una red de interacciones de proteínas implicadas en la enfermedad de Alzheimer y su uso, junto con una base de datos de fármacos y sus efectos terapéuticos, para predecir el efecto de nuevos compuestos para su tratamiento, así como para replantearse el uso de otros. Parte de este trabajo está publicado aquí.

Ana Rojas nos explicó como su grupo había reconstruido la filogenia de la superfamilia de proteínas RAS, encontrando por el camino los residuos funcionales responsables de las diferencias funcionales de las diferentes familias.

Mar Gonzàlez  nos resumió resultados de su reciente artículo sobre la variabilidad del splicing alternativo en poblaciones humanas.

Alberto Pascual García nos mostró cómo, a partir de datos de presencia/ausencia de rRNA en muestras de diferentes ambientes, se puede inferir la composición bacteriana en un ambiente y por medio de una aproximación basada en la arquitectura de las redes resultantes estudiar si los diferentes géneros tienden a agregarse o segregarse.

Hernán Dopazo nos contó resultados de un trabajo suyo que está en revisión donde sostiene que la composición de los genomas, en cuanto a elementos como genes, rRNAs, promotores, elementos repetidos, etc, se distribuye de manera parecida a la distribución de especies que se observa en ecosistemas naturales, en un proceso donde aparentemente la selección tiene poco que decir.

En la última charla que pude escuchar Jaime Huerta nos contó los progresos de su grupo para hacer filogenias anidadas, que se pueden entender como un proceso recursivo donde vamos refinando el árbol inicial recalculando de manera recursiva la topología de las ramas a medida que vamos de la raíz a las hojas, añadiendo nuevos genes ortólogos a medida que avanza el proceso.

Entre los pósters del segundo día tomé nota del servidor iLOOP para la predicción de interacciones proteína-proteína, del software TAPyR para el alineamiento de reads largos como los de 454 y del programa Pyicos para el procesamiento de datos de ChIPseq.

PD Una oferta de trabajo que se publicó en el congreso:

The Evolutionary Genomics Group in the Comparative and Computational Genomics program of the IBE (http://www.ibe.upf-csic.es/) is willing to recruit a PhD student. More information is available in the attached document. For queries, please contact Tomás Marquès-Bonet (tomas.marques@upf.edu)

13 de enero de 2012

Construir un suffix array: Manber y Myers

Hola,
hoy me gustaría volver un poco sobre el tema de los suffix array que ya ha sido tratado alguna vez en el blog.

En 1990, Manber y Myers publicaban Suffix arrays: a new method for on-line string searches, como motivo del primer ACM-SIAM symposium on Discrete algorithms. En éste fantástico e ilustrativo artículo, introducen la comparación de esta estructura con los suffix trees y proponen definiciones, notación y algoritmos para su construcción, búsqueda y optimización.

La mayor desventaja de suffix trees reside en el espacio requerido, y con los suffix arrays se consigue reducir considerablemente, siendo aún así competitivos frente a los suffix trees en términos de búsqueda y construcción.

Aquí veremos la primera parte, la creación de un suffix array ordenado (sección 3 del artículo original, sorprendentemente), con un algoritmo que requiere O(NlogN) como peor caso, si bien con optimizaciones se podría conseguir un tiempo esperado O(N); lo cual ya se consiguió de manera formal en 2003.

Quizás hoy día los suffix array van perdiendo protagonismo de primera plana, pero siguen siendo la base para muchas aplicaciones, por ejemplo implicados en métodos de obtención de la secuencia Barrows-Wheeler. Posiblemente en entradas posteriores hablemos de mejoras en el uso de suffix arrays y de otras estructuras relacionadas como FM-index y BWT.

Por ahora, aquí os dejo el código en Python (-V 2.6.6). Podéis cambiar la secuencia a indexar en la variable inicial A.

Y nos podemos preguntar ¿por qué usar éste código en lugar de la simple línea de ordenamiento MergeSort que aparecía en el código Perl de la entrada anterior sobre suffix arrays? Adjunto una tabla que lo explica en términos de tiempo requerido para la construcción del índice.


Long texto Myers MergeSort
20 0''019 0''004
100 0''023 0''005
1,000 0''055 0''013
11,000 0''717 0''401
110,000 22''200 41''075
220,000 1'36''004 2'42''189
660,000 9'31''297 25'54''598


Por supuesto estamos espectantes ante cualquier comentario, crítica o sugerencia.
¡Un saludo!

 #!/usr/bin/python

# Simple algorithm for Suffix Array creation,
# based on Manber & Myers, 1990
# by Cantalapiedra, 2011

############### FUNCTIONS ##############
## Function to print a Suffix Array
def printSA(text, sa = [], sa_name = "Suffix Array"):
  if len(sa) == 0:
       print "No records in {0}".format(sa_name)
  else:
       print sa_name
       for i in sa:
            print "{0}\t{1}".format(i, text[i:])
  print "\n"
  return

## Function to print one or several lists
def printTemp(things = [], title = "##", post_title = "#"):
  if len(things) == 0:
       print "No things to print"
  else:
       print title
       for a in things:
            print a
       print post_title
  print "\n"
  return

############### VARIABLES ##############

A = 'TTTTAGATCGATCGACTAGA$'

pos = [] # Contains the SA based on suffix numbers (SN)
prm = [] # Inverse of pos, each position is a SN containing its position in SA
bH = [] # Boolean array marking the buckets

b2H = [] # Temporary array for marking the 2H-buckets
count = [] # Counts the number of SN in a 2H-bucket

## All above variables have length = |A|

################ START ###################
# Some verbose
print "Your text is {0}\n".format(A[:len(A) - 1])

# Initial SA (POS)
pos = [(a,i) for i,a in enumerate(A)]
pos = sorted(pos)
pos = [a[1] for a in pos]

printSA(A, pos, "INITIAL SA")

# Initialize PRM
prm = [0]*len(A)

# Initial positions and buckets
prev_letter = ""
for i, pos_i in enumerate(pos):
  prm[pos_i] = i

  letter = A[pos_i]
  if letter != prev_letter:
       bH.append(1)
  else:
       bH.append(0)
  prev_letter = letter

######### H-STAGES ##########

count_changed = []

h_stage = 1
while h_stage < len(A): # max log2(|A| + 1)

  #print "stage {1} init {0}".format(time.time(), h_stage) IMPORT time module to use this

  # FIRST: re-assign prm as of buckets starting points
  #            and create count and b2H arrays
  prev_one = 0

  for i, bh_i in enumerate(bH):
       if bh_i == 1:
            prev_one = i # Is a bucket starting position
     
       if bh_i == 0:
            prm[pos[i]] = prev_one # Points to its bucket starting pos
     
       count.append(0)
       b2H.append(0)

  # Reset count[] to 0s
  for a in range(len(count)):
       count[a] = 0

  # THEN: for every bucket...
  # Here the entire SA is scanned, but in fact takes account
  # of the current bucket being covered
  prev_count = [a for a in count]
  for i, pos_i in enumerate(pos):
     
       if bH[i] == 1:
            b2H[i] = 1
     
       # We are moving the H-prefix (Ti) of suffix in i position
       ti = pos_i - h_stage
       if ti < 0: continue
     
       ### This is the heart of the algorithm
       # If its the first change of this bucket during this H-bucket run...
       if prev_count[prm[ti]] == count[prm[ti]]:
            change = True
       else:
            change = False
     
       count[prm[ti]] += 1
       count_changed.append(prm[ti])
       prm[ti] = prm[ti] + count[prm[ti]] - 1
     
       # ... that means this is a new 2H-bucket
       if change:
            b2H[prm[ti]] = 1
     
       #printTemp([prm, count, b2H], "After moving suff: {0}".format(ti), "")
     
       # If next position is another H-bucket, or is the last bucket
       if (i < len(bH) - 1 and bH[i + 1] == 1) or i == len(bH) - 1:
          
            #printTemp([prm, count, b2H], "", "NEXT H-BUCKET {0}".format(i + 1))
          
            for j in count_changed:
                 prev_count[j] = count[j]
            count_changed = []
     
       # If there are no buckets left
       if 0 not in b2H:
            break

  #print "end stage {0}".format(time.time())

  # If there are no buckets left
  if 0 not in b2H:
       break

  #printTemp([prm, count, b2H], "", "NEXT ROUND {0}".format(h_stage*2))

  # Update POS from PRM
  for j, prm_j in enumerate(prm):
       pos[prm_j] = j

  # Update BH from B2H
  bH = b2H

  # Reset vars for next H-STAGE
  count = []
  b2H = []
  h_stage *= 2

# Update POS from PRM
for j, prm_j in enumerate(prm):
  pos[prm_j] = j

printSA(A, pos, "FINAL SA")
#


La salida del ejemplo es:

Your text is TTTTAGATCGATCGACTAGA

INITIAL SA
$
AGATCGATCGACTAGA$
ATCGATCGACTAGA$
ATCGACTAGA$
ACTAGA$
AGA$
A$
CGATCGACTAGA$
CGACTAGA$
CTAGA$
GATCGATCGACTAGA$
GATCGACTAGA$
GACTAGA$
GA$
TTTTAGATCGATCGACTAGA$
TTTAGATCGATCGACTAGA$
TTAGATCGATCGACTAGA$
TAGATCGATCGACTAGA$
TCGATCGACTAGA$
TCGACTAGA$
TAGA$

FINAL SA
$
A$
ACTAGA$
AGA$
AGATCGATCGACTAGA$
ATCGACTAGA$
ATCGATCGACTAGA$
CGACTAGA$
CGATCGACTAGA$
CTAGA$
GA$
GACTAGA$
GATCGACTAGA$
GATCGATCGACTAGA$
TAGA$
TAGATCGATCGACTAGA$
TCGACTAGA$
TCGATCGACTAGA$
TTAGATCGATCGACTAGA$
TTTAGATCGATCGACTAGA$
TTTTAGATCGATCGACTAGA$

29 de diciembre de 2011

Jornadas Bioinformáticas JBI 2012 (XI Edición)

Aprovecho esta última entrada del año para dar difusión a las Jornadas de Bioinformática JBI 2012. Como muchos sabréis, este congreso es el principal punto de encuentro anual de nuestra comunidad en la península ibérica, así que nuestro laboratorio también estará en Barcelona del 23 al 25 de Enero. El programa completo de las jornadas se puede descargar en este enlace.
http://sgu.bioinfo.cipf.es/jbi2012


Este año presentaremos parte de nuestro trabajo reciente:

"Genome-wide clustering of transcription factors by comparison of predicted protein-DNA interfaces"

donde explicamos y evaluamos la anotación de interfaces de reconocimiento de DNA en secuencias de proteínas por medio de diferentes aproximaciones como BLAST, TFmodeller, DP-Bind y DISIS.






El tema principal de las jornadas será "Arquitectura genómica, anotación y diseño", sobre el cual se discutirán los diferentes avances en la integración de los campos de la Biología, Medicina e Informática en el campo de la Genómica. Además se tratarán los siguientes temas:
- Análisis de datos de secuenciación de alto rendimiento (NGS)
- Bioinformática estructural
- Algoritmos de biología computacional y computación de alto rendimiento
- Análisis de sequencias, filogenética y evolución
- Bases de datos, herramientas y tecnologías en biología computacional
- Bioinformática en transcriptómica y proteómica
- Biología de sistemas

ENGLISH:

The XIth Spanish Symposium on Bioinformatics (JBI2012) will take place in January 23-25, 2012 in Barcelona, Spain. Co-organised by the Spanish Institut of Bioinformatics and the Portuguese Bioinformatics Network and hosted by the Barcelona Biomedical Research Park (PRBB). The full program can be downloaded from this link.

This year, the reference topic is “Genome Architecture, Annotation and Design” for which the conference will provide the opportunity to discuss the state of the art for the integration of the fields of biology, medicine and informatics. We invite you to submit your work and share your experiences in the following topics of interest including, but not limited to:

- Analysis of high throughput data  (NGS)
- Structural Bioinformatics
- Algorithms for computational biology and HPC
- Sequence analysis, phylogenetics and evolution
- Databases, Tools and technologies for computational biology
-  Bioinformatics in Transcriptomics and Proteomics
- System and Synthetic Biology

Our contribution to the congress:

Genome-wide clustering of transcription factors by comparison of predicted protein-DNA interfaces

Transcription Factors (TFs) play a central role in gene regulation by binding to DNA target sequences, mostly in promoter regions. However, even for the best annotated genomes, only a fraction of these critical proteins have been experimentally characterized and linked to some of their target sites. The dimension of this problem increases in multicellular organisms, which tend to have large collections of TFs, sometimes with redundant roles, that result of whole-genome duplication events and lineage-specific expansions. In this work we set to study the repertoire of Arabidopsis thaliana TFs from the perspective of their predicted interfaces, to evaluate the degree of DNA-binding redundancy at a genome scale. First, we critically compare the performance of a variety of methods that predict the interface residues of DNA-binding proteins, those responsible for specific recognition, and measure their sensitivity and specificity. Second, we apply the best predictors to the complete A.thaliana repertoire and build clusters of transcription factors with similar interfaces. Finally, we use our in-house footprintDB to benchmark to what extent TFs in the same cluster specifically bind to similar DNA sites. Our results indicate that there is substantial overlap of DNA binding specificities in most TF families. This observation supports the use of interface predictions to construct reduced representation of TF sets with common DNA binding preferences.


23 de diciembre de 2011

Felices Navidades bioinformáticas!!!

Ya que estamos en Navidad os dejo nuestra felicitación navideña especialmente pensada para los bioinformáticos... con algo de inspiración en estructura tridimensional de proteínas y dianas terapéuticas.

Adjunto la imagen a alta resolución para que podáis intentar adivinar con qué estructura del PDB se han generado las bolas (proteínas) y adornos (ligandos) que decoran el árbol, dejar en los comentarios vuestras respuestas...

Y por supuesto... FELIZ NAVIDAD Y PRÓSPERO AÑO 2012!!!

English version:

In this Christmas time we want to share our Christmas card especially designed for bioinformaticians ... with some inspiration from three-dimensional structure of proteins and therapeutic targets.

It is a high resolution image so you can try to guess from which PDB structure have been extracted the balls (proteins) and ornaments (ligands) that decorate the tree, give your answers in the comment section ...

And of course... MERRY CHRISTMAS AND HAPPY NEW YEAR 2012!!!









19 de diciembre de 2011

metaDBSite: un metaservidor para predecir la unión de proteínas a DNA

Los meta-servidores son valiosas herramientas en bioinformática para no perder el tiempo mandando trabajos a diferentes servidores y tener que juntar después manualmente los resultados. Un meta-servidor nos ofrece utilizar varios servicios web a la vez desde una interfaz única. Ejemplos de meta-servidores son: NCBI, UniProt, GeneSilico...

Recientemente he  tenido que usar diferentes programas de predicción de interfaces de unión proteína-DNA, y he encontrado la grata sorpresa de que existía un meta-servidor que permitía realizar la predicción de una sola vez utilizando la mayoría de técnicas conocidas, este valioso meta-servidor se llama metaDBSite. No sólo eso, si no que he podido comprobar que siguen dándole soporte.

Os animo a probarlo con una sencilla secuencia proteica:
>1a02_N  
 wplssqsgsyelrievqpkphhRahYetEgsRgavkaptgghpvvqlhgymenkplglqifigtaderilkphafyqvhritgktvtttsyekivgntkvleiplepknnmratidcagilklrnadielrkgetdigRkntrvrlvfrvhipessgrivslqtasnpiecsQRsahelpmverqdtdsclvyggqqmiltgqnftseskvvftekttdgqqiwemeatvdkdksqpnmlfveipeyrnkhirtpvkvnfyvingkrkrsqpqhftyhpv  

Si comprobáis los resultados con los datos de la estructura original (cadena N de la estructura 1a02) podréis ver cómo el mejor de los métodos (BindN-RF) 'adivina' 7 de los 7 residuos que contactan el DNA entre los 18 predichos, si tenemos en cuenta que la proteína tienen 280 aminoácidos, el resultado no está nada mal...