26 de diciembre de 2020

Pangenomas para dummies: ejemplos en plantas, aplicaciones y retos

 Hola, el 17 de diciembre di una charla virtual invitado por mis compañeros de la Estación Experimental de Aula Dei-CSIC titulada "Pangenomas para dummies: ejemplos en plantas, aplicaciones y retos". 

https://chilmedia.org/v2/media/c046b9fc-9684-49f5-834d-42cf6828d7a0.jpg

Un pangenoma se define como la unión de todos los genomas de una especie. El análisis de pangenomas es una herramienta habitual en microbiología, sobre todo en bacterias para analizar su ecología y patogenicidad. En cambio, nuestro conocimiento de los pangenomas de plantas es todavía limitado. En la charla repaso resultados recientes en mono y dicotiledóneas, incluyendo nuestro trabajo reciente sobre el híbrido Brachypodium hybridum y sus progenitores, que demuestran la utilidad de esta aproximación para explorar la diversidad genética en poblaciones y bancos de
germoplasma. En cuanto a la mejora, el reto es cómo distinguir los genes accesorios que contribuyen a la adaptación de aquellos que son reliquias evolutivas.

La charla quedó grabada y puedes escucharla de nuevo en 

https://balanbbb.corp.csic.es/playback/presentation/2.0/playback.html?meetingId=9ea1a71e86c265e77c16b8078be749094142699f-1608199170685

Hasta pronto,

Bruno


30 de noviembre de 2020

AlphaFold resuelve el plegamiento de proteínas (en CASP14)

Hola, 

estos días está transcurriendo CASP14, la edición 14 del certámen de predicción  de estructura de proteínas. La última vez que hablamos de CASP en este blog fue en esta entrada del año pasado. Por recordar un poco, en CASP participan grupos de investigación de todo el mundo que tratan de modelar un conjunto de secuencias de proteínas cuyas estructuras se han resuelto experimentalmente, pero que solamente se publican despúes de la temporada de predicción. Por tanto, los grupos y su algoritmos trabajan relativamente a ciegas en esas predicciones.

Digo relativamente porque en realidad se apoyan en la creciente colección de estructuras conocidas del PDB, del orden 10E5,  y en las millones de secuencias de proteínas conocidas (del orden 10E8). Por esa razón unas secuencias son más fáciles, porque se parecen a otras conocidas, y otras más difíciles, porque no se parecen demasiado a nada conocido.

En la edición CASP14 había un total de 92 secuencias de aminoácidos, cada una correspondiente a un dominio. La siguiente figura, obtenida de https://predictioncenter.org/casp14/zscores_final.cgi resume los resultados,  mostrando que los dos mejores grupos de la última década (Baker y Zhang) han sido ampliamente superados por AlphaFold2 (columna de la izquierda, del que ya habíamos hablado aquí):

 


En definitiva, la combinación de estrategias de aprendizaje automático de AlphaFold2, descritas en https://deepmind.com/blog/article/alphafold-a-solution-to-a-50-year-old-grand-challenge-in-biology, han superado con mucho a todos los expertos que más saben de este problema tan difícil. 

Qué opinan los expertos? Aquí tenéis por ejemplo a Torsten Schwede , Mohammed AlQuraishi o a Alfonso Valencia

Supongo que no estará de más esperar a CASP15 para comprobar que este progreso se mantiene en el tiempo, pero por ahora parece que lo han resuelto. Solamente nos queda exigir a DeepMind, la matriz de AlphaFold2, que libere su predictor para fines académicos y de esa manera devuelvan a la comunidad lo que de ella han extraído en estos últimos años.

Un saludo,

Bruno

PD Nota importante: como recuerda Alfonso Valencia en https://twitter.com/Alfons_Valencia/status/1333682759366303745, no es lo mismo predecir la estructura que la reacción de plegamiento


4 de noviembre de 2020

Course on scripting with the Linux shell

 Hi,  Carlos Cantalapiedra and me recently put together teaching material about scripting in the linux terminal.

The material can be found at repository https://github.com/eead-csic-compbio/scripting_linux_shell 

There are five sessions and the goal is for you to learn the basics of the Linux shell and scripting for data sciences such as genomics and plant breeding:

session title required time URL
0 Setup prior to course session 0
1 Linux basics and files 2h session 1
2 Processes and scripts 2h session 2
3 Parsing with regular expressions 2h session 3
4 Perl one-liners 2h session 4
5 Advanced scripts 2h session 5 


Figure of the standard streams, taken from https://en.wikipedia.org/wiki/Standard_streams

 

If you spot errors please send pull requests, hope this helps some of you out there,

Bruno 

PD si prefieres aprender en español echa un vistazo a https://github.com/vinuesa/intro2linux


































21 de octubre de 2020

Proteínas como autómatas finitos

 Hola,

hoy quería comentar un artículo reciente de Bozovic et al,  donde observan el alosterismo de un dominio PDZ al unirse a su péptido ligando (ver más aquí). Lo que me ha llamado más la atención de este trabajo es que combinando medidas experimentales y simulaciones de dinámica molecular, definen poblaciones de conformaciones y cronometran las transiciones, de manera que al final proponen un modelo que parece un autómata finito, como los que se usan en teoría de la computación o para explicar expresiones regulares.

Figura 1. Cambios en la conformación del dominio proteico PDZ2.
Izq: unido a su ligando (trans). Dcha: tras liberarse (cis). Las distancias 
entre los residuos 20,71 y 4,55 sirven para hacer un seguimiento de
las conformaciones en experimentos y simulaciones de dinámica 
molecular. Figura tomada de https://www.pnas.org/content/117/42/26031

Figura 2. Autómata que describe los 5 estados/conformaciones de PDZ2.
El diámetro de cada estado es proporcional a su frecuencia en las
poblaciones moleculares. Las flechas son las transiciones entre estados en ms.
 

Al final, resulta que las herramientas desarrolladas hace décadas en las ciencias de la computación sirven hoy para describir poblaciones de proteínas. O dicho de otro modo, esta proteína computa,

hasta luego,

Bruno


15 de octubre de 2020

Plant Genomes in a Changing environment 2020

Hi, this year's Plant Genomes in a Changing Environment was virtual with hashtag #PlantGenomes20.  Check the full program here. I post here my notes on the talks I was able to attend, either live or on the recorded videos, in case they're useful for anyone out there.

Giles Oldroyd (Sainbury lab, UK) talked about a pathway conserved in plants forming intracellular endosymbiotic mycorrhizaephytes. The pathway is not found in species forming exclusively extracellular symbioses, such as ectomycorrhizae, and those forming associations with cyanobacteria. This is described at https://www.nature.com/articles/s41477-020-0613-7. Receptors used with potential symbionts are the same used to recognize chitin and peptide-glicans of pathogens. Symbiotic potential decreases when there's enough nutriends at hand.

 

Figure from https://www.biorxiv.org/content/10.1101/804591v1.full

 

N Marmiroli (SITEIA.PARMA, Italy) presented a poster about the https://simbaproject.eu, where they are testing in plots in Italy and Germany the application of Plant Growth-Promoting Microbes to wheat, maize, tomato and potato under abiotic stress conditions.

Trevor Nolan (Duke U, USA) talked about their efforts in harmonizing single cell mRNA sequencing data from 110K cells to construct an atlas for the developing Arabidopsis thaliana root. Read more at https://www.biorxiv.org/content/10.1101/2020.06.29.178863v1


Francesco Licausi (Oxford U, UK) described the strategies of flooding tolerance in plants, which can be used to breed flooding varieties of crops. Flooding is a composite stress, he focus on two components: oxygen sensing and responses to oxidative stress, which are mediated among others by NAC transcription factors such as ANAC017, read more at https://onlinelibrary.wiley.com/doi/full/10.1111/pce.13037 and https://onlinelibrary.wiley.com/doi/abs/10.1111/tpj.14975


 Ana Caño Delgado (CRAG, Spain)
delivered a lecture on physiology and molecular biology of Arabidopsis plants and their responses to drought mediated by brassinosteroids receptors in the root tip. Among many experiments, they measure response to drought as function of root branching angles. Read more at https://www.nature.com/articles/s41467-018-06861-3 and https://science.sciencemag.org/content/368/6488/266.abstract

Steve Long (Lancaster U, UK) shared his knowledge on breeding for incresing photosynthesis efficienty. This figure is from a 2010 review https://www.annualreviews.org/doi/10.1146/annurev-arplant-042809-112206, but still quite impressive:



Alexander Borowsky (U California, Riverside, USA) uses their own coexpression network pipeline (Wisecaver 2019) and the Taiji pipeline to analyze their ATACseq data in order to find key transcription factors regulating responses to water logging and water deficit in rice. They plan to do validation work using CRISPR to confirm their results.

Julia Bailey-Serres (U California, Riverside, USA) gave a keynote with two parts. While the second was related to A Borowsky0s talk, on the first part, she and her team looked at the transcriptional regulation of responses to submergence/flooding by comparing 4 different species (rice, Medicago truncatula, domesticated tomato and its dryland-adapted wild relative Solanum pennellii). By combining RNA-seq and ATACseq data and different computational analyses they found 68 Submerged Upregulated oRthologous Families (SURFs) upregulated in both monocots and dicots, with hypoxia-responsive promoter element (HRPE) motifs being abundant in rice but scarce in promoters of SURFs in dryland-adapted S. pennellii. This piece of work is described at https://science.sciencemag.org/content/365/6459/1291

Figure from https://science.sciencemag.org/content/365/6459/1291

They have really nive figure, shown below, where DNA motifs instead of transcription factors are the vertices of the graph. Despite the differences observed in the network across species, 4 motifs seem to be highly conserved:

Figure from https://science.sciencemag.org/content/365/6459/1291

 

Klaas Vandepoele (VIB-Ugent Center for Plant Systems Biology, Belgium) and his team have been doing a conceptually similar analysis, with RNA-seq from A. thaliana, Populus trichocarpa and maize, in order to identify gene triplets with conserved leaf expression pattern which are conserved across monocots and dicots. They validate their results by phenotyping A. thaliana mutants in the glass house.  Their work has been published at https://onlinelibrary.wiley.com/doi/full/10.1111/pbi.13223. He also mentions this paper https://pubmed.ncbi.nlm.nih.gov/31748815 which describes modules of genes involved in cell proliferation in the leaf.

Jeff Habben (Corteva, USA) describes their extensive field trials in low and high-yielding locations, mostly across the US, to check whether different transgenic lines that overexpress maize transcription factor zmm28 (MADS-box) have higher and more stable yields. Their work is described in detail at https://www.pnas.org/content/116/47/23850

Facundo Romani (CONICET, Argentina) explains how liverwort Marchantia polymorpha defends itself against hervibores by accumulating terpenes in oil bodies, a process controlled by HD-ZIP transcription factor MpC1HDZ. This work is described at https://www.sciencedirect.com/science/article/abs/pii/S0960982220307685

Sandra Knapp (Natural History Museum-London, UK) , a plant taxonomist, talks about diversity of Solanaceae (nightshades), starting with potatoes

Source: Sandra Knapp

She discusses how we usually think of aridification as something bad but according tho their work in Australia (https://onlinelibrary.wiley.com/doi/full/10.1111/jse.12638) it seems to be certainly a driver of diversification. She suggests we should be looking at plant lineages that have this pattern to prepare for climate change. One possible way would be to learn from resistant weeds in order to create resistant crops (https://www.sciencedirect.com/science/article/abs/pii/S0304423818307994). This is aligned with the work of Dana MacGregor, which also had a talk in this session.

She terminates by acknowledging the collective effort of people studying nightshades around the world (http://solanaceaesource.org) and a warning about the excess of P/N consum by agriculture and the current rates of extinction, that cause an unsustainable loss of diversity.

Michael Purugganan (New York University, USA) talks about the domestication of rice (which tokk 4K years) and molecular studies to understand the geographic spread of primarily japonica rice in Asia. His two main figures are:

 Figure from https://pubmed.ncbi.nlm.nih.gov/28087768/


 

Tal Dahan Meir (Weizmann Institute of Science, Israel) spoke about a long term observation seris of a wild Emmer wheat population in Isreal. Unfortunately I could not see the most relevant slides in the saved video, but het talk reminded me of the Eviatar Nevo's talk last year at the Brachypodium conference (see https://www.youtube.com/watch?v=UhNzH39zLbk).

Guy Naamati (Ensembl Plants, EMBL-EBI, UK) presented a poster about the wheat Pan Genome and status at Ensembl Plant. The programmatic access recipes mentioned in the poster can be found at https://github.com/Ensembl/plant_tools/tree/master/demo_scripts

Click on the image to download and see full size

6 de octubre de 2020

Violines desde varios ficheros en R

 Hola,

la semana pasada necesitaba hacer una gráfica para mostrar lado a lado varias distribuciones de miles de observaciones. Para ello me decanté por las llamadas gráficas violín, que mi colega Pablo Vinuesa ya había usado en esta figura (panel D) hace unos años.


 Los violin plots permiten resumir distribuciones de manera concisa, mostrando la mediana (punto blanco), el rango intercuartil (el segmento negro grueso) y la densidad de puntos como una curva a cada lado. Por esta última característica son superiores a los diagramas de caja (boxplots). 

El siguiente código en R muestra como calculé esta gráfica horizontal usando la librería vioplot de R a partir de múltiples ficheros de texto, de un valor por línea y extensión .tsv, uno por cada serie de datos.

library(vioplot)

# setwd("path") # if required
filedir="./"

# parse input TSV file names
repeat_files = list.files(path=filedir, pattern="\\.tsv")
series_names = gsub("\\.tsv", "", repeat_files)

# actually read files into data frames
repeats = lapply(repeat_files, function(i){
  log10(read.table(i, header=FALSE))
})
names(repeats) <- series_names

# increase left and bottom margins to make room for axis labels
par(mar = c(4, 11, 1, 1)) 

plot("", 
  ylim = c(0.5, length(repeats)+0.5), 
  # low X values truncated in example   
  xlim = c(2, max(unlist(repeats))), 
  yaxt = "n",  
  ylab = "", 
  xlab = "log10 repeat length")
axis(2, labels = series_names, at = c(1:length(repeats)), las=1)

# add violins one by one
lapply(seq_along(repeats), function(x)
  vioplot(repeats[[x]], at = x,  add = T, box = F, horizontal=T)
)

Hasta pronto,

Bruno

12 de septiembre de 2020

Alineamiento global con el algoritmo Wavefront (WFA)

Hola, en esta entrada vuelvo a una de las piedras angulares de la biología computacional, el alineamiento de secuencias. Este problema tiene múltiples caras, algunas ya discutidas en este blog, pero desde el algoritmo de Needleman-Wunsch su formulación fundamental es el alineamiento global de dos secuencias q y t de longitudes n y m, desde el principio hasta el final. 

La última vuelta de tuerca la acaban de publicar Santiago Marco Sola y colaboradores en Bioinformatics, donde describen su algoritmo de alineamiento llamado wavefront (WFA). En el artículo los autores muestran como WFA y su variante heurística WFA-Adapt son más eficientes tanto en tiempo de cálculo como en consumo de memoria que las alternativas del estado del arte, incluyendo los algoritmos de la librería KSW2, empleada por  minimap2 (visto en este blog).

Recomiendo la lectura del artículo porque es muy didáctico y está escrito en un estilo sencillo, dada la naturaleza del problema. A continuación resumo aquí las ideas principales. 

El problema que resuelve WFA es un alineamiento global entre dos secuencias usando el modelo affine-gap como función coste. Esto significa que para cada par de posiciones alineadas el coste del alineamiento aumenta en 0 puntos en caso de ser idénticas (a), en x puntos en caso de ser diferentes (x=6 en BWA-MEM) o con un coste lineal para las inserciones que se calcula como g(n) = o + e n  donde o es el coste de abrir un indel y e el de extenderlo n bases. WFA es un algoritmo exacto que calcula el alineamiento óptimo como aquel que termina en la celda (n,m) de la matriz de programación dinámica (Figura 1, izquierda) con un coste total más pequeño. Hasta aquí nada nuevo, es un algoritmo que busca diagonales que alinean las dos secuencias.

La novedad de WFA es que define furthest-reaching points (fr), vectores  Fs,k que indican para una diagonal k el punto más lejano donde se alcanza un coste s (ver Figura 1 izquierda, vectores M0, M4 y M8 desde el origen, donde M=matches, de la misms manera que I=indel y D=deletion). En su algoritmo reformulan el alineamiento por programación dinámica calculando vectores fr para un coste s en base a los vectores fr calculados para costes menores, pero de manera que sólo una fracción de los fr se llegan a calcular. El algoritmo descrito en la Figura 2 recibe su nombre porque para cada coste s se define el frente de onda WFs como el conjunto de todos los vectores fr con coste total s. El alineamiento optimo se corresponde a la secuencia de frentes de onda desde WF0 a WFs que alcanzan la coordenada (n, m) con el menor coste s. En el ejemplo de la Figura 1 solamente es necesario guardar en memoria 3 WF (0, 4 y 8) para calcular el alineamiento óptimo y luego reconstruirlo. A diferencia de otros algoritmos de alineamiento, WFA es más eficiente cuánto más se parecen las secuencias a alinear y sus operaciones son fácilmente paralelizables de manera portable con instrucciones SIMD.

 
Figura 1. Alineamiento global de dos secuencias en una matriz de programación dinámica (izq) y su representación en forma de vectores fr (furthest-reaching points, dcha). En la matriz las celdas (0,0) y (6,6) marcan respectivamente el principio y el final del alineamiento. Tomada de Bioinformatics


Figura 2. Pseudocódigo para el algoritmo WFA, tomado de Bioinformatics.

Para terminar esta entrada, el código de WFA está escrito en C y se compila fácilmente con gcc si lo clonas desde el repositorio https://github.com/smarco/WFA . Allí encontrarás ejemplos sencillos de cómo llamar a las funciones de alineamiento, un saludo,

Bruno

 

7 de septiembre de 2020

Intento explicar Perl 7

Buenas,

a finales de Junio participé en la más reciente conferencia sobre Perl, https://perlconference.us/tpc-2020-cloud . La principal novedad, sobre la que ya twiteé en aquellas fechas, fue el anuncio de una nueva versión de Perl, en concreto Perl7. Han pasado ya dos meses y no hay mucho más que contar, pero parece que podemos ir haciéndonos a la idea de qué significa esto. 

Perl 7 tendrá como punto de partida la última versión de Perl5, en concreto v5.32, publicada en Enero de este año. La principal diferencia es que tendrá por defecto varias características que eran opcionales hasta ahora pero que la comunidad de desarrolladores considera son ya necesarias para modernizar el lenguaje. Éstas siguien en discusión, pero parece que van a incluir las siguientes (mira también estas FAQ):

  • enable strict by default
    • % perl -Mstrict program.pl
  • enable warnings by default
    • % perl -Mwarnings program.pl
  • disable bareword filehandles
    • % perl -M-bareword::filehandles
  • disable multidimensional array emulation (a Perl 4 trick)
    • % perl -M-multidimensional program.pl
  • enable subroutine signatures
    • % perl -Mfeature=signatures program.pl
  • change prototypes to use the :prototype attribute

 

La ley del mínimo esfuerzo

Sólo tendrás que poner use v5.32 al principio de tus scripts.

Qué pasa si no quiero cambiar nada?

Al parecer las versiones v5.30 y v5.32 tendrán soporte unos diez años más, esa es la promesa.

Qué pasa con Perl6 o raku?

Se considera que es tan distinto a Perl que es un lenguaje independiente, que tendrá su propia vida dentro de la familia.

Cuando haya más novedades iré contándolas por aquí, hasta luego,

Bruno

31 de agosto de 2020

Ubuntu nativo en windows 10

Hola, 

durante años he utilizado MobaXterm como herramienta para poder conectarme a servidores Linux por SSH y SFTP desde mi portátil con Windows 10. Esta herramienta te ofrece un terminal Linux, con casi todas las características que esperas, dentro de un sistema Windows. Además, te permite instalar componentes de Linux que van más allá del terminal BASH. Por ejemplo, te permite instalar paquetes de python o perl de manera sencilla. 

Pego aquí una imagen de mi MobaXterm en la carpeta /mnt/c, que corresponde al disco C:\


Pero de lo que  quería hablar hoy en realidad es de que Windows 10 ya soporta Ubuntu de manera nativa y de esa manera podemos tener un sistema Linux completo, no solamente un terminal, sin salir de nuestra sesión Windows. Para instalar debes seguir estos pasos:

1) Asegúrate de que tienes una versión de Windows10 con la versión 1709 o superior, como se explica en https://www.protocols.io/view/ubuntu-on-windows-for-computational-biology-sfuebnw 

2) En la tienda de Microsoft (Microsoft Store) busca e instala el software "Ubuntu". Hay otras opciones, yo elegí Ubuntu a secas.

3) Busca en el sistema la aplicación "Activar o desactivar las características de Windows" (en inglés Turn Windows Features on or off) y selecciona el subsistema de Windows para Linux.

4) Reinicia y arranca la aplicación Ubuntu. Tras unos minutos de instalación deberás elegir tu usuario y seña para Linux. Después deberás ver algo como:

Puedes ver que las particiones de mi sistema Windows (C y Q) se montan de manera automática en /mnt, como en MobaXterm, y que mi home de Linux está en una partición diferente: /home.

Ahora ya puedes hacer algo tan importante, por ejemplo, como instalar un compilador: sudo apt-get install gcc

Hasta pronto,

Bruno