Alphafold2: estructuras de proteínas para todos los públicos

Hola,

el jueves pasado escribía aquí sobre RoseTTAFold, que ponía al alcance de los programadores una potente herramienta para predecir la estructura de cualquier proteína que tenga suficientes secuencias parecidas en los archivos. Unas horas más tarde me enteré de que su principal competidor, AlphaFold2, del que hemos escrito aquí, aquí y aquí , había hecho también sus deberes y publicaba sus modelos neuronales para predecir el plegamiento de secuencias de aminoácidos. Hay muchos hilos en Twitter, como éstos: 1 y 2 .

En un artículo en Nature repasan sus resultados en CASP14, como su notable precisión de 1.5 Å de RMSD para todos los átomos modelados, y muestran que esos resultados excelentes se sostienen en gran medida en un conjunto de 3144 estructuras recientemente publicadas en el Protein Data Bank tras descartar las estructuras conocidas con % identidad de secuencia > 40. El siguiente vídeo resume el proceso de predicción de una secuencia con múltiples dominios y 863 aminoácidos. Recuerdo que la única información de partida es la secuencia de aminoácidos y un alineamiento múltiple de esa secuencia con otras secuencias no redundantes de proteínas homólogas:

 

Fuente: https://www.nature.com/articles/s41586-021-03819-2

Alphafold2 no es perfecto, ya que en la validación contra el PDB hay un en torno a un 25% de predicciones con valores de RMSD elevados, pero en la mayor parte de los casos es excelente:

Fuente: https://www.nature.com/articles/s41586-021-03819-2

Una de las fortalezas del método es que es capaz de estimar bien la calidad los modelos que genera (ver Figura 2b de https://www.nature.com/articles/s41586-021-03828-1). Entre las limitaciones reseñadas por los autores del trabajo, destacan:

1) La calidad de las predicciones cae de manera significativa si el alineamiento múltiple de secuencias homólogas a la de interés tiene una profundidad < 30 (leer más aquí). 

2) Las estructuas de Alphafold son de menor calidad para dominios poco compactos, con pocos contactos, como los dominios puente. Esto contrasta con las buenas prestaciones observadas (3) para predecir estructuras cuaternarias de homómeros.

Para los que queráis probarlo hay varias opciones:

i) El contenedor Docker descrito en https://github.com/deepmind/alphafold que ojo, requiere 2.2TB de espacio si instalas todas las bases de datos. 

ii) Un cuaderno Colab con un predictor simplificado en https://colab.research.google.com/github/deepmind/alphafold/blob/main/notebooks/AlphaFold.ipynb

iii) Las predicciones ya disponibles en UniProt para un total de 21 especies (https://alphafold.ebi.ac.uk/download), incluyendo 4 plantas (Arabidopsis thaliana, soja, arroz y maíz). Puedes leer más en 4, 5 y 6 y 7.

Hasta pronto,

Bruno

Introducción a AlphaFold

Hola,
vengo de escuchar a Andrew W Senior en el auditorio Kendrew del EMBL-EBI hablar de cómo funciona AlphaFold para predecir la estructura de proteínas. Como titulaba en Diciembre, todavía no sabemos plegar proteínas, pero cada vez tenemos más información en el Protein Data Bank (del orden de 150K estructuras) para poder predecir distancias entre residuos y ángulos diedros.

Fuente: https://deepmind.com/blog/alphafold

Trataré de resumir lo que nos ha contado antes de que salga el artículo, pero ojo, ha advertido que el código no piensan publicarlo por ahora. Su sistema comprende varios módulos que se ejecutan secuencialmente:

1) Extración de correlaciones evolutivas entre residuos de una secuencia de aminoácidos en forma de modelo de Potts (CCMpred) a partir de perfiles de secuencias homólogas obtenidas con PSIBLAST o HHblits, como hace HHpred.

2) Red neuronal profunda, con al menos dos variantes:
2.1) predicción de distancias reales entre C-betas, no contactos, a partir de histogramas precalculados en el rango de 2 a 22 Ansgtrom. En esto se parecen a RaptorX.
2.2) predicción de ángulos diedros phi y psi

NOTA1: Las predicciones de distancias les permiten asignar estructura secundaria con una precisión Q3  del 84% usando las estructuras de CASP11.

NOTA2: Parten la matriz de contactos en submatrices (crops) que comprenden siempre un trozo de la diagonal principal y esa manera, y con la estructura secundaria, propagan localmente las restricciones de distancias entre residuos cercanos.

3) Diferenciación del potencial de distancias/ángulos por métodos de minimización de gradientes. Lo hacen partiendo la secuencia en dominios o con la secuencia entre y observan que la segunda manera funciona mejor.

4) Relajan el esqueleto obtenido y le añaden cadenas laterales completas con Rosetta. En sus manos esto no mejora el modelo de manera significativa.

Hasta luego,
Bruno

Tabla de estabilidad de sustituciones de aminoácidos

Hola,
a menudo nos preguntamos qué efecto puede tener un cambio de aminoácido sobre la estabilidad de una proteína. Una pregunta relacionada es si tendrá fenotipo una sustitución no sinónima (missense). De hecho hay toda una colección de herramientas que se han desarrollado para esto (ver por ejemplo el curso de bioinformática estructural). El trabajo reciente de Rocklin et al (2017), donde estudian la estabilidad de miles de miniproteínas (hasta 50aa) sintéticas y unas 500 naturales del PDB, nos proporciona la mejor respuesta hasta la fecha. Aunque es un artículo complejo, en esencia lo que hacen es medir la estabilidad de miles de proteínas expresándolas y exponiéndolas a proteasas en levadura.
De esa manera pueden estimar el efecto que tienen mutaciones individuales dependiendo de su contexto de estructura secundaria, ya sea una alfa-hélice, una lámina beta o un lazo o loop:

Efecto sobre la estabilidad de mutaciones en diferentes contextos estructurales. Los valores negativos, como los de la prolina en general, son desestabilizadores. Adaptada de http://science.sciencemag.org/content/357/6347/168.full.

En la siguiente figura se muestran ejemplos de mutaciones desestabilizadoras en amarillo:

Adaptada de http://science.sciencemag.org/content/357/6347/168.full.

Hasta pronto,
Bruno

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

Leyendo el diagrama de Ramachandran

Hola,
la semana pasada les planteé a mis alumnos de la Licenciatura en Ciencias Genómicas el problema de asignar estructura secundaria a los residuos de una proteína, una vez medidos los ángulos dihedros del esqueleto peptídico psi y phi. Para esta tarea lo lógico es usar el diagrama de Ramachandran, como éste tomado de la wikipedia:

La solución naive consiste en delimitar las regiones del diagrama por medio de condiciones tales como 'si phi > X  AND psi < Y' entonces un residuo está en estado Z. Sin embargo, para tener mucha precisión sería necesario tener muchas condiciones como ésta, dada las formas irregulares del diagrama.


Una solución interesante, propuesta por F.Peñaloza, consiste en realmente leer el mapa, como lo hacemos nosotros. Lo que él hizo fue obtener un mapa de Ramachadran cuadrado con unos cuantos colores conocidos,

Diamagra de Ramachandran RGB de F.Peñaloza.

Diagrama anterior con los colores RGB empleados

para luego convertirlo a texto ASCII (con algo como text-image). Con este mapa en memoria, es sencillo obtener la estructura secundaria para unas coordenadas phi,psi dadas. 

Aprovechando el módulo GD de perl podemos leer la imagen directamente, como muestra el siguiente código:

 use strict;  
 use GD;  
   
 my $RAMAPLOT = 'ramachandran.png';  
 my $plot = GD::Image->newFromPng($RAMAPLOT);  
 my ($maxpsi,$maxphi) = checkPlot($plot);  
   
 my $inputPSI = 160;   
 my $inputPHI = -45;   
 printf("El estado de estructura secundaria para %s,%s es: %s\n",  
    $inputPSI,$inputPHI, getSimpleSS($maxpsi,$maxphi, $inputPSI,$inputPHI));  
   
 sub getSimpleSS  
 {  
   my ($maxpsi,$maxphi,$psi,$phi) = @_;  
   my %colorSS = (  
     '236,30,36','E',  
     '164,74,164','E',  
     '156,218,236','H',  
     '4,162,236','H',  
     #'180,230,28','L'   
     );  
   
   if($psi > 180 || $psi < -180 || $phi > 180 || $phi < -180)  
   {  
     die "# valid psi/phi are [-180,+180]\n";  
   }  
   $psi = (0.5 - (0.5 * $psi/180)) * $maxpsi;  
   $phi = (0.5 + (0.5 * $phi/180)) * $maxphi;  
     
   my @RGB = $plot->rgb( $plot->getPixel($phi,$psi) );  
   return $colorSS{"$RGB[0],$RGB[1],$RGB[2]"} || 'C';  
 }  
   
 sub checkPlot  
 {  
   my ($plot,$print_pixels) = @_;  
   my ($maxphi,$maxpsi) = $plot->getBounds();  
   if(!$print_pixels){ return ($maxphi,$maxpsi) }  
     
   foreach my $phi ( 0 .. $maxphi )   
   {  
     foreach my $psi ( 0 .. $maxpsi )   
    {  
       my @RGB = $plot->rgb( $plot->getPixel($phi,$psi));  
       print "$RGB[0] $RGB[1] $RGB[2]\n";  
    }  
   }  
   return ($maxphi,$maxpsi);  
 }  

Un saludo, Bruno

'Fold it' o como el juego agudiza el ingenio

Buenas,
desde hace al menos una década el grupo de David Baker (U.Washington) ha tenido un impacto tremendo en el estudio del plegamiento de las proteínas, yo diría que hasta el punto de ir cocinando poco a poco un premio Nobel, pero el tiempo dirá. Muy brevemente, su trabajo ha consistido en desarrollar algoritmos para la predicción de la estructura de proteínas (por comparación y ab initio) que han ido refinando en sucesivas rondas de CASP (donde siempre quedan entre los mejores) y poniendo a punto con experimentos bastante espectaculares donde comparan sus predicciones con estructuras obtenidas por cristalografía y NMR, obteniendo modelos de una calidad y complejidad sin precedentes, llegando por ejemplo a diseñar una proteína artificial y predecir su estructura terciaria con un error promedio de 1.2A. Una de las aportaciones más extravagantes de este laboratorio tan productivo fue la creación de un juego basado en el problema del plegamiento, que llamaron FoldIt.

El juego permite poner a jugadores de todo el mundo a resolver ejercicios de plegamiento de dificultad variable, tras superar un tutorial donde aprendes las reglas del juego y las herramientas disponibles. La verdad es que hasta hoy no he encontrado una manera más amena de entender qué es un puente de hidrógeno o un choque estérico. Enfin, a lo que voy, lo que parecía simplemente un intento por popularizar el 'problema del plegamiento', que a día de hoy todavía se considera sin resolver, se ha convertido también en una herramienta de innovación, como se publica en el último número de PNAS. Tras estudiar la evolución de las estrategias ganadoras encontradas por la comunidad jugadores, que se pueden guardar como recetas, se observa que se parecen mucho a los propios algoritmos que el grupo de Baker está desarrollando. O dicho de otro modo, parecería que poniendo en forma de juego un problema complejo como esto, con puntuaciones objetivas y reglas sencillas, el ingenio colectivo puede llegar a soluciones parecidas a las de los científicos especialistas, sorprendente?
Un saludo,
Bruno

Alineamiento de perfiles con HHalign

Hola,
hoy discutiré un problema habitual cuando trabajamos con familias de proteínas, el de alinear entre si dos grupos de secuencias previamente alineadas (MSAs=multiple sequence alignments) , como se ilustra en el esquema siguiente:

Dadas las secuencias del grupo 1, que hemos alineado entre si con ayuda de programas como muscle, clustal o MAFFT, y las secuencias del grupo 2, que hemos procesado de manera similar, ahora queremos ver cómo se alinean todas juntas, por ejemplo para inferir diferentes funciones o estructuras.

En el caso trivial sería tan sencillo como poner todas las secuencias en el mismo archivo y hacer un alineamiento múltiple de una vez. Sin embargo, el caso que discutimos aquí es más complicado, el que se da cuando no es sencillo alinear una secuencia del grupo 1 con otra del 2, porque han divergido bastante. Entonces podemos probar a convertir el alineamiento MSA1 a un perfil y a alinearlo contra el perfil correspondiente a MSA2, como si fuera un alineamiento pareado. Podemos incorporar, además de la información evolutiva capturada en cada perfil, información de estructura secundaria para ayudar a guiar el alineamiento.

Todo esto es sencillo de realizar con el paquete de programas HHsearch, que podemos descargar de ftp://toolkit.lmb.uni-muenchen.de/HHsearch.

Los pasos a seguir son 3:

1) Obtener archivos FASTA de MSA1.faa (30 secuencias) y MSA2.faa (35 secuencias). Para mayor precisión en el alineamiento es recomendable añadir en ambos archivos una predicción de estructura secundaria, en el formato de PSIPRED. Por ejemplo, el contenido de MSA1.faa comenzaría con las siguientes líneas:

>aa_pred
MVVSIGVFDGVHI--GHQKVL...
>ss_pred
CEEEEECCCCEEH--HHHHHH...
>ss_conf
8888617775007--889999...
>seq1 (y todas las 29 siguintes a continuación)
MVVSIGVFDGVHI--GHQKVL...

2) Convertir MSA1.faa y MSA2.faa en perfiles de Markov (MSA1.hhm y MSA2.hhm), con el ejecutable hhmake:

$ ~/soft/HHsearch/hhmake -i MSA1.faa -seq 30
$ ~/soft/HHsearch/hhmake -i MSA2.faa -seq 35

3) Alinear MSA1 contra MSA2, ya sea de manera global o local:

$ ~/soft/HHsearch/hhalign -i MSA1.hhm -t MSA2.hhm -glob -seq 65 -ofas globalMSA.faa

La salida obtenida es algo como:

REMARK: in -mac -global mode -mact is forced to 0
Query file is in HHM format
Read in HMM seq1 with 467 match states and effective number of sequences = 3.1
Query file is in HHM format
Read in HMM seq31 with 576 match states and effective number of sequences = 1.8
Using maximum accuracy (MAC) alignment algorithm ...
Printing alignments in FASTA format to globalMSA.faa 

Aligned seq1 with seq31: Score = 218.04   P-value = 1     
Done

Un ejemplo de aplicación real de esta técnica de alineamiento lo tenéis en el artículo http://www.biomedcentral.com/1471-2148/10/311 .