Biología universitaria — tercer año · Bachelor Year 3

5Bioinformática y análisis de secuencias

Un biólogo que acaba de secuenciar un gen de un gusano abisal pega sus $300$ aminoácidos en un formulario web y, tres segundos después, descubre que la proteína es prima lejana de una quinasa humana, con un $31\,\%$ de identidad a lo largo de $280$ residuos y una probabilidad de $10^{-40}$ de que el parecido sea azar. Detrás de esos tres segundos hay un algoritmo de programación dinámica de 1970, una teoría estadística de los alineamientos aleatorios, matrices de sustitución destiladas de miles de familias de proteínas y una base de datos de unos cien mil millones de residuos. Este capítulo trata del razonamiento que hay dentro de la caja: cómo se alinean dos secuencias de manera que el alineamiento sea demostrablemente el mejor, cómo se logra que la puntuación signifique algo, cómo se distingue una coincidencia de una casualidad y cómo se encuentran patrones en un genoma que nadie ha mirado antes. La matemática es elemental — una recurrencia, un logaritmo, una distribución de Poisson — y conviene conocerla, porque toda conclusión extraída de una comparación de secuencias descansa en ella.

5.1 Alinear dos secuencias

Definición 5.1 (Alineamiento y puntuación)

Un alineamiento de dos secuencias las escribe una sobre la otra, con huecos (–) insertados de modo que las columnas emparejen un residuo con un residuo o un residuo con un hueco, y ninguna columna empareje dos huecos. Su puntuación es la suma, sobre las columnas, de una puntuación de sustitución $s(a,b)$ por cada par de residuos y de una penalización por hueco por cada hueco: una penalización lineal $-d$ por posición de hueco o, de forma más realista, una penalización afín $-d - (k-1)e$ por una tirada de $k$ huecos, con un coste de apertura $d$ mayor que el coste de extensión $e$ , ya que una sola inserción de varios residuos es un único suceso evolutivo. Un alineamiento global cubre ambas secuencias de extremo a extremo; un alineamiento local halla el par de subcadenas de mayor puntuación e ignora el resto, que es lo que uno quiere cuando un dominio compartido se sitúa en dos proteínas por lo demás no emparentadas.

Teorema 5.2 (Needleman–Wunsch)

Sean $x = x_{1}\dots x_{m}$ e $y = y_{1}\dots y_{n}$ , con penalización lineal por hueco $d$ . Defínase $F(i,j)$ como la mejor puntuación de un alineamiento global de los prefijos $x_{1}\dots x_{i}$ e $y_{1}\dots y_{j}$ . Entonces $F(i,0) = -id$ , $F(0,j) = -jd$ y, para $i,j \ge 1$ ,

F(i,j) = \max\bigl\{\,F(i-1,j-1) + s(x_{i},y_{j}),\; F(i-1,j) - d,\; F(i,j-1) - d\,\bigr\}.

$F(m,n)$ es la puntuación global óptima, un alineamiento óptimo se recupera siguiendo hacia atrás desde $(m,n)$ las decisiones que produjeron cada máximo, y todo el cálculo cuesta $mn$ pasos. La variante de Smith–Waterman para el alineamiento local añade $0$ como cuarta opción del máximo, pone los bordes a $0$ y lee la respuesta en la mayor entrada de la tabla.

Demostración. Considérese la última columna de cualquier alineamiento de los dos prefijos. Es una de tres cosas: $x_{i}$ sobre $y_{j}$ , $x_{i}$ sobre un hueco o un hueco sobre $y_{j}$ . Quitarla deja un alineamiento de $(x_{1}\dots x_{i-1}, y_{1}\dots y_{j-1})$ , de $(x_{1}\dots x_{i-1}, y_{1}\dots y_{j})$ o de $(x_{1}\dots x_{i}, y_{1}\dots y_{j-1})$ respectivamente, cuya puntuación es como mucho $F$ de ese par; y, a la inversa, cada uno de esos alineamientos óptimos puede extenderse con la última columna correspondiente. Así que la mejor puntuación que termina en cada tipo de columna es $F$ del par más corto más la puntuación de la columna, y el óptimo es el mayor de los tres. Los bordes están forzados (contra un prefijo vacío solo caben huecos). La inducción sobre $i + j$ rellena la tabla; el número de celdas es $(m+1)(n+1)$ . Para el alineamiento local la opción adicional $0$ significa “empiécese aquí un alineamiento nuevo”, lo que hace de $F(i,j)$ la mejor puntuación de un alineamiento que termina en $(i,j)$ , y el mejor alineamiento local termina en algún sitio. ∎

Ejemplo 5.3 (Una tabla de cuatro por tres)

Alinéese GAT con GCAT, puntuando $+1$ una coincidencia, $-1$ una discrepancia y $d = 1$ . Los bordes son $0, -1, -2, -3, -4$ a lo largo de la fila superior y $0, -1, -2, -3$ por el lado. Rellenando fila a fila: $F(\text{G},\text{G}) = 1$ , $F(\text{G},\text{C}) = 0$ , $F(\text{G},\text{A}) = -1$ , $F(\text{G},\text{T}) = -2$ ; $F(\text{A},\text{G}) = 0$ , $F(\text{A},\text{C}) = 0$ , $F(\text{A},\text{A}) = 1$ , $F(\text{A},\text{T}) = 0$ ; $F(\text{T},\text{G}) = -1$ , $F(\text{T},\text{C}) = -1$ , $F(\text{T},\text{A}) = 0$ , $F(\text{T},\text{T}) = 2$ . El óptimo es $2$ , y siguiendo hacia atrás — diagonal desde (T,T), diagonal desde (A,A), luego a la izquierda de (G,C) a (G,G) y después diagonal — se obtiene

\begin{array}{c} \texttt{G-AT}\\ \texttt{GCAT} \end{array}

tres coincidencias y un hueco: $3 - 1 = 2$ .

La tabla de Needleman–Wunsch para GAT contra GCAT (coincidencia +1, discrepancia -1, hueco -1). Cada celda es la mejor puntuación de los dos prefijos que terminan ahí; el camino rojo trazado hacia atrás desde la esquina es el alineamiento óptimo. — La tabla de Needleman–Wunsch para GAT contra GCAT (coincidencia $+1$ , discrepancia $-1$ , hueco $-1$ ). Cada celda es la mejor puntuación de los dos prefijos que terminan ahí; el camino rojo trazado hacia atrás desde la esquina es el alineamiento óptimo.

Método 5.4 (Alinear dos secuencias)

(1) Elíjase la puntuación: una matriz de sustitución adecuada a la divergencia esperada (BLOSUM62 para proteínas de distancia desconocida; coincidencia y discrepancia para el ADN) y penalizaciones afines por hueco (típicamente apertura $-11$ y extensión $-1$ con BLOSUM62). (2) Decídase global o local: global para dos secuencias que se creen homólogas en toda su longitud, y local en los demás casos. (3) Rellénese la tabla con la recurrencia, guardando en cada celda un puntero a la decisión que dio su máximo. (4) Trácese hacia atrás desde $(m,n)$ (en el global) o desde la celda máxima hasta un cero (en el local), escribiendo el alineamiento de derecha a izquierda. (5) Júzguese el resultado no por su puntuación bruta sino por su significación estadística (más abajo), y mírese: los huecos largos, las tiradas de baja complejidad y los alineamientos confinados a una repetición son avisos.

5.2 Puntuar: cuánto vale una coincidencia

Definición 5.5 (Matrices de sustitución)

Una matriz de sustitución da $s(a,b)$ para cada par de aminoácidos como una puntuación de log-verosimilitudes:

s(a,b) = \frac{1}{\lambda}\,\log\frac{q_{ab}}{p_{a}\,p_{b}},

donde $q_{ab}$ es la frecuencia con la que $a$ y $b$ aparecen alineados en alineamientos fiables de proteínas emparentadas, $p_{a} p_{b}$ la frecuencia con la que se emparejarían por azar y $\lambda$ una escala elegida para que las entradas sean enteros cómodos. Una puntuación positiva significa que el par aparece más a menudo en homólogos que por azar; las puntuaciones de identidad son mayores para los aminoácidos raros (triptófano $+11$ , cisteína $+9$ en BLOSUM62) y menores para los frecuentes (leucina $+4$ , alanina $+4$ ), y las sustituciones conservadoras (isoleucina–valina $+3$ ) puntúan positivo mientras que las radicales (triptófano–glicina $-2$ ) puntúan negativo. Las matrices PAM (Dayhoff, 1978) se dedujeron de proteínas muy próximas y se extrapolaron a distancias mayores por multiplicación de matrices; las matrices BLOSUM (Henikoff y Henikoff, 1992) se contaron directamente en bloques de secuencias alineadas agrupadas a una identidad dada — BLOSUM62 a partir de bloques al $62\,\%$ — y son las predeterminadas porque se midieron, en lugar de extrapolarse, a la distancia a la que se usan.

Proposición 5.6 (Por qué log-verosimilitudes)

Para que un esquema de puntuación pueda usarse en alineamiento local, la puntuación esperada de una columna emparejada al azar, $\sum_{a,b} p_{a} p_{b}\, s(a,b)$ , ha de ser negativa, y algunas puntuaciones han de ser positivas; de lo contrario los alineamientos aleatorios crecerían sin límite y el segmento de mayor puntuación sería la secuencia entera. Dado eso, todo esquema de esa clase es equivalente a un esquema de log-verosimilitudes para algunas frecuencias objetivo $q_{ab}$ — los alineamientos que encontrará como óptimos son aquellos cuyos pares de residuos se distribuyen como $q_{ab}$ . Elegir la matriz es por tanto elegir la divergencia que se espera detectar: una matriz para parientes próximos (BLOSUM80, PAM30) tiene positivos más marcados y negativos más duros, y una para parientes lejanos (BLOSUM45, PAM250) es más plana.

Demostración. Admitido a este nivel. ∎

Ejemplo 5.7 (Identidad, similitud y la zona crepuscular)

Dos secuencias de proteína aleatorias alineadas de forma óptima con huecos alcanzan por azar entre el $15\text{ a }20\,\%$ de identidad. Por encima del $35\,\%$ de identidad a lo largo de cien residuos, dos proteínas son casi con seguridad homólogas; entre el $20\,\%$ y el $35\,\%$ está la zona crepuscular, donde la identidad por sí sola no decide y debe hacerlo la estadística de más abajo. Los homólogos pueden caer muy por debajo de la zona: las subunidades de la hemoglobina y la mioglobina comparten un $25\,\%$ de identidad, la lisozima y la $\alpha$ -lactalbúmina un $40\,\%$ , y muchos pares de proteínas con el mismo plegamiento comparten menos del $15\,\%$ , detectable solo comparando perfiles o estructuras.

5.3 Buscar en una base de datos

Definición 5.8 (BLAST)

Alinear una consulta de $300$ residuos contra una base de datos de $10^{11}$ por programación dinámica completa costaría $3\times 10^{13}$ actualizaciones de celda por búsqueda. BLAST (Altschul y colaboradores, 1990) cambia un poco de sensibilidad por mil veces más velocidad en tres pasos: (1) se listan las palabras de la consulta (tres residuos para proteínas, once bases para ADN) y sus vecinas de alta puntuación; (2) se rastrea la base de datos en busca de coincidencias exactas de palabra — las semillas; (3) se extiende cada semilla en ambas direcciones sin huecos hasta que la puntuación cae una cantidad fijada por debajo de su máximo, conservando los pares de segmentos de alta puntuación (HSP), y luego se unen los HSP próximos con programación dinámica con huecos en una banda estrecha. Un homólogo verdadero contiene casi siempre al menos una palabra exacta de tres residuos en común; un parecido casual rara vez la tiene, y nunca se extiende.

La heurística de BLAST. Las palabras exactas cortas compartidas por la consulta y la entrada de la base de datos (rojo) son semillas; cada una se extiende a lo largo de su diagonal mientras la puntuación sigue subiendo, y solo las extensiones que se mantienen altas llegan a ser pares de segmentos de alta puntuación.

Teorema 5.9 (La estadística de un acierto por azar)

Para una consulta de longitud $m$ buscada en una base de datos de longitud total $n$ , con un esquema de puntuación de esperanza negativa, el número de alineamientos locales sin huecos que puntúan al menos $S$ y que surgen por azar sigue una Poisson de media

E = K\,m\,n\,e^{-\lambda S},

donde $\lambda$ y $K$ solo dependen del esquema de puntuación y de las frecuencias de residuos ( $\lambda$ es la escala de la matriz de log-verosimilitudes). $E$ es el valor esperado de la puntuación $S$ . Escribiendo la puntuación en bits, $S' = (\lambda S - \ln K)/\ln 2$ , la fórmula queda $E = m n\, 2^{-S'}$ , y la probabilidad de que al menos un alineamiento por azar alcance $S$ es $P = 1 - e^{-E}$ , que vale $E$ cuando $E$ es pequeño.

Demostración parcial. La cola exponencial es el teorema de Karlin–Altschul y se admite: la puntuación máxima de un segmento en un paseo aleatorio con deriva negativa tiene una distribución cuya cola decae como $e^{-\lambda S}$ , con $\lambda$ la raíz positiva de $\sum_{a,b} p_{a} p_{b} e^{\lambda s(a,b)} = 1$ — que es exactamente la ecuación que hace coherente la matriz de log-verosimilitudes. Dada esa cola, el resto es contar. Los segmentos de alta puntuación pueden empezar en cualquiera de los $mn$ pares de posiciones, son raros y son casi independientes; el número de los que superan $S$ sigue por tanto una Poisson de media proporcional a $mn$ y a la probabilidad de la cola, $E = Kmn\,e^{-\lambda S}$ . La probabilidad de que no haya ninguno es $e^{-E}$ . La sustitución por la puntuación en bits es álgebra: $e^{-\lambda S} K = 2^{-(\lambda S - \ln K)/\ln 2}$ . Para los alineamientos con huecos vale la misma forma con $\lambda$ y $K$ estimados por simulación. ∎

Ejemplo 5.10 (Leer un valor E)

Una consulta de $250$ residuos contra una base de datos de $5\times 10^{10}$ residuos tiene $mn = 1.25\times 10^{13} \approx 2^{43.5}$ . Un acierto con una puntuación en bits de $60$ tiene $E = 2^{43.5 - 60} = 2^{-16.5} \approx 10^{-5}$ : es homólogo con certeza prácticamente absoluta. Un acierto con $S' = 40$ tiene $E = 2^{3.5} \approx 11$ : se esperan once puntuaciones así por azar, y el acierto no significa nada. El mismo alineamiento, con la misma puntuación en bits, buscado en una base de datos diez veces mayor, tiene un $E$ diez veces mayor — la significación es una propiedad de la búsqueda, no del par. El umbral de uso común es $E < 10^{-3}$ para un homólogo fiable; $E \approx 0.01$ – $1$ merece una segunda mirada con un método de perfil.

E = mn\,2-S': cada bit adicional reduce a la mitad el número esperado de aciertos por azar, y una base de datos diez veces mayor cuesta 3.3 bits de significación para el mismo alineamiento. — $E = mn\,2^{-S'}$ : cada bit adicional reduce a la mitad el número esperado de aciertos por azar, y una base de datos diez veces mayor cuesta $3.3$ bits de significación para el mismo alineamiento.

5.4 Perfiles, estados ocultos y motivos

Definición 5.11 (Alineamiento múltiple y perfiles)

Un alineamiento múltiple de secuencias dispone una familia de secuencias en columnas de residuos homólogos. La programación dinámica exacta sobre $k$ secuencias cuesta $n^{k}$ y es imposible más allá de tres; los programas prácticos alinean de forma progresiva, primero el par más próximo según un árbol guía y luego secuencias y grupos al alineamiento creciente, con rondas de refinamiento. Un alineamiento terminado se resume en un perfil: para cada columna, la frecuencia de cada residuo y de los huecos. Un modelo oculto de Markov de perfil formaliza esto como una cadena de estados de coincidencia, uno por columna conservada, cada uno emitiendo residuos con sus propias probabilidades, con estados de inserción y de deleción que permiten residuos adicionales o ausentes en cada posición; el modelo de una familia (una entrada de Pfam) puntúa una secuencia nueva por la probabilidad del mejor camino a través de los estados, y encuentra homólogos muy por debajo de la zona crepuscular de la comparación por pares, porque una columna que solo tolera residuos hidrófobos lo dice, mientras que una secuencia aislada no puede.

Un modelo oculto de Markov de perfil de una familia de cuatro columnas. Cada estado de coincidencia M emite un residuo con las frecuencias propias de la columna; los estados de inserción I (con bucles sobre sí mismos) admiten residuos adicionales, y los de deleción D se saltan una columna. Puntuar una secuencia es hallar su camino más probable.

Definición 5.12 (Motivos y contenido de información)

Un motivo es un patrón corto — un sitio de factor de transcripción, una señal de corte, un sitio de fosforilación — representado por una matriz de pesos por posición con la frecuencia $f_{i}(b)$ de cada base o residuo $b$ en cada posición $i$ . El contenido de información de la posición $i$ es $R_{i} = 2 - H_{i}$ bits para el ADN, donde $H_{i} = -\sum_{b} f_{i}(b)\log_{2} f_{i}(b)$ es su entropía: $2$ bits para una base invariante y $0$ para una posición en la que las cuatro son igualmente probables. El total $R = \sum_{i} R_{i}$ se dibuja como un logotipo de secuencia, con cada posición una pila de letras cuya altura total es $R_{i}$ y cuyas letras se dimensionan por frecuencia.

Proposición 5.13 (Cuánta información necesita un sitio)

Un sitio que debe encontrarse $\gamma$ veces en un genoma de $G$ posiciones, y en ningún otro sitio, necesita unos $R_{\text{nec}} = \log_{2}(G/\gamma)$ bits de contenido de información: el motivo ha de reducir las $G$ posiciones candidatas a las $\gamma$ verdaderas, y cada bit reduce a la mitad las candidatas. Los motivos observados de los reguladores bacterianos bien estudiados coinciden con esta predicción — los sitios de E. coli de un represor que se une a unas pocas decenas de lugares en un genoma de $4.6\,\mathrm{Mb}$ llevan $16\text{ a }18$ bits —; los motivos de los factores de transcripción eucariotas, con $8\text{ a }12$ bits en un genoma de $3\times 10^{9}$ , no pueden especificar sus dianas por sí solos, y de ahí que actúen en combinaciones y en la cromatina abierta del Capítulo 1.

Demostración. Una posición al azar coincide con un motivo de contenido de información $R$ con probabilidad cercana a $2^{-R}$ (cada bit de especificidad reduce a la mitad la probabilidad), de modo que el número esperado de coincidencias por azar en $G$ posiciones es $G\,2^{-R}$ . Para que los sitios verdaderos destaquen, esto ha de ser del orden de $\gamma$ o menos: $G\,2^{-R} \le \gamma$ , es decir, $R \ge \log_{2}(G/\gamma)$ . ∎

Ejemplo 5.14 (Coincidencias esperadas por azar)

Un sitio de restricción de seis bases fijas tiene $R = 12$ bits y coincide con una posición al azar con probabilidad $4^{-6} = 2^{-12}$ : unas $1100$ veces en un genoma de E. coli de $4.6\,\mathrm{Mb}$ leído en ambas hebras (el sitio es palindrómico, de modo que una vez por posición), y $7\times 10^{5}$ veces en el genoma humano. Un factor eucariota cuyo motivo lleve $10$ bits coincide con $3\times 10^{9}\times 2^{-10} \approx 3$ millones de posiciones del genoma humano, varios miles de veces más que los genes que regula. Un motivo por sí solo es un predictor débil en un genoma grande; lo que permite una predicción son el estado de la cromatina, los motivos vecinos y la conservación del sitio entre especies.

Un logotipo de secuencia de un motivo de promotor parecido a la caja TATA. La altura de cada pila es el contenido de información de esa posición, 2 - H_i bits; las cuatro primeras posiciones son casi invariantes y llevan la mayor parte de los 12 bits del motivo. — Un logotipo de secuencia de un motivo de promotor parecido a la caja TATA. La altura de cada pila es el contenido de información de esa posición, $2 - H_{i}$ bits; las cuatro primeras posiciones son casi invariantes y llevan la mayor parte de los $12$ bits del motivo.

5.5 De la secuencia a la función

Método 5.15 (Anotar una proteína desconocida)

Dada una secuencia codificante nueva: (1) tradúzcase en el marco correcto y búsquense un péptido señal, segmentos transmembranarios y regiones de baja complejidad; (2) búsquese en las bases de datos de proteínas con BLAST y léanse los aciertos con $E < 10^{-3}$ , anotando si el alineamiento cubre toda la proteína (un ortólogo verdadero) o un segmento (un dominio compartido); (3) búsquese en las bases de datos de dominios con HMM de perfil, que encuentran familias que BLAST pierde y reparten la proteína en dominios; (4) infiérase ortología, y no mera similitud, comprobando que el mejor acierto en el otro genoma tiene la consulta como su mejor acierto (mejores aciertos recíprocos) o situando la proteína en un árbol génico (Capítulo 25); (5) transfiérase la función de los ortólogos con cautela — un residuo catalítico conservado apoya que la química se conserve, y uno ausente lo desmiente — y predígase la estructura; (6) trátese toda predicción como una hipótesis para el laboratorio.

Proposición 5.16 (La estructura a partir de la secuencia)

El plegamiento de una proteína lo determina su secuencia (Capítulo 7), y calcularlo a partir de la secuencia fue durante cincuenta años el problema central sin resolver del campo. Tres enfoques lo consiguieron, uno tras otro. El modelado por homología construye la estructura de una proteína sobre la de un homólogo resuelto, de forma fiable por encima del $30\,\%$ de identidad. El análisis de coevolución aprovecha que dos residuos en contacto en el plegamiento tienden a mutar juntos a lo largo de un alineamiento múltiple profundo, de modo que los pares de columnas acoplados estadísticamente son contactos predichos, y con suficientes contactos queda definido un plegamiento. Los métodos de aprendizaje profundo, entrenados con las cien mil estructuras resueltas y con esos alineamientos, predicen hoy la mayoría de las estructuras de proteínas globulares con una exactitud casi experimental (las evaluaciones CASP de 2020), y las bases de datos guardan una estructura predicha para prácticamente toda secuencia de proteína conocida. Lo que predicen peor es lo que una sola estructura no capta: las regiones desordenadas, las conformaciones alternativas, el efecto de una mutación puntual y los complejos.

Izquierda: una estructura de proteína predicha, coloreada según la confianza del modelo, de alta (azul) a baja (naranja) en un bucle desordenado. Derecha: un despacho de bioinformática — navegadores de genomas y árboles en las pantallas, y ni un solo poyo de laboratorio a la vista.

Observación 5.17 (Los límites de la inferencia)

La mayoría de las anotaciones funcionales de las bases de datos nunca se pusieron a prueba; se transfirieron de un homólogo, que a su vez había sido anotado por transferencia. Los errores se propagan y se multiplican, y una anotación equivocada en una proteína bien conectada puede infectar a toda una familia. Los remedios son los de más arriba: distíngase la ortología de la homología, léase el alineamiento, búsquense los residuos catalíticos y recuérdese que “proteína hipotética” es una etiqueta honrada que un tercio de los genes de la mayoría de los genomas todavía merece.

5.6 Ejercicios

Ejercicio 5.1 ★

Defínanse el alineamiento global y el local y dese una situación biológica que reclame cada uno.

Solución

Solución de Ejercicio 5.1.

Global: las dos secuencias alineadas de extremo a extremo, con todo residuo en una columna — para dos proteínas que se creen homólogas en toda su longitud, como los ortólogos de una enzima constitutiva. Local: el par de subcadenas de mejor puntuación, ignorando el resto — para encontrar un dominio compartido (un dominio SH2 en dos proteínas de señalización por lo demás no emparentadas), o un gen en una secuencia genómica larga.

Ejercicio 5.2 ★

Rellénese la tabla de Needleman–Wunsch de AGC contra AAC con coincidencia $+1$ , discrepancia $-1$ y hueco $-1$ , y dense el alineamiento óptimo y su puntuación.

Solución

Solución de Ejercicio 5.2.

Bordes $0,-1,-2,-3$ en ambos sentidos. Fila A: $1, 0, -1$ . Fila G: $0, 0, -1$ . Fila C: $-1, -1, 1$ . Óptimo $F(3,3) = 1$ : AGC sobre AAC sin huecos (coincidencia, discrepancia, coincidencia: $1 - 1 + 1 = 1$ ).

Ejercicio 5.3 ★

En BLOSUM62, triptófano–triptófano puntúa $+11$ y leucina–leucina $+4$ . Explíquese a partir de la fórmula de log-verosimilitudes por qué la identidad del residuo más raro vale más.

Solución

Solución de Ejercicio 5.3.

$s(a,a) = \lambda^{-1}\log\bigl(q_{aa}/p_{a}^{2}\bigr)$ . El triptófano es raro ( $p_{W} \approx 0.013$ ), de modo que la probabilidad de que dos triptófanos se alineen al azar, $p_{W}^{2}$ , es minúscula, y un par de triptófanos conservado es un indicio de homología mucho más fuerte que un par de leucinas conservado ( $p_{L} \approx 0.1$ ); la razón de log-verosimilitudes es proporcionalmente mayor.

Ejercicio 5.4 ★

¿Qué es un valor E? Una búsqueda devuelve un acierto con $E = 3$ . ¿Qué significa ese número, y es homólogo el acierto?

Solución

Solución de Ejercicio 5.4.

El valor E es el número de alineamientos con una puntuación al menos tan alta que cabría esperar por azar en una búsqueda de esta consulta contra una base de datos de este tamaño. $E = 3$ significa que se esperan tres puntuaciones así por azar: el acierto no es prueba de homología (puede serlo de todos modos, pero la búsqueda no lo puede decir).

Ejercicio 5.5 ★★

Se busca una consulta de $400$ residuos contra $2\times 10^{11}$ residuos. Calcúlense los valores E de los aciertos con puntuaciones en bits $45$ , $55$ y $65$ . ¿Qué puntuación en bits da $E = 10^{-3}$ ? ¿Cómo cambia la respuesta si la consulta mide $40$ residuos?

Solución

Solución de Ejercicio 5.5.

$mn = 400\times 2\times 10^{11} = 8\times 10^{13} = 2^{46.2}$ . $E(45) = 2^{1.2} \approx 2.3$ ; $E(55) = 2^{-8.8} \approx 2\times 10^{-3}$ ; $E(65) = 2^{-18.8} \approx 2\times 10^{-6}$ . $E = 10^{-3}$ exige $S' = 46.2 + 10.0 = 56$ bits. Una consulta de $40$ residuos tiene un $mn$ diez veces menor, $2^{42.9}$ : bastan $53$ bits — pero una consulta corta rara vez llega siquiera a eso.

Ejercicio 5.6 ★★

Calcúlese el contenido de información de un motivo cuyas cuatro posiciones tienen frecuencias de bases (A, C, G, T) de $(1,0,0,0)$ , $(0.5,0,0.5,0)$ , $(0.25,0.25,0.25,0.25)$ y $(0.7,0.1,0.1,0.1)$ . ¿Cuántas coincidencias por azar tiene en un genoma de $4.6\,\mathrm{Mb}$ ?

Solución

Solución de Ejercicio 5.6.

Contenidos de información: $2$ , $1$ , $0$ , y $2 - H$ con $H = -(0.7\log_{2} 0.7 + 3\times 0.1\log_{2} 0.1) = 0.36 + 1.00 = 1.36$ , o sea $0.64$ . Total $R = 3.64$ bits. Coincidencias por azar: $9.2\times 10^{6}$ posiciones en dos hebras $\times 2^{-3.64} \approx 7\times 10^{5}$ — el motivo es casi inútil por sí solo.

Ejercicio 5.7 ★★

Explíquese por qué las penalizaciones afines por hueco son más realistas que las lineales, y por qué una penalización de apertura muy alta y una muy baja dan las dos malos alineamientos.

Solución

Solución de Ejercicio 5.7.

Una inserción de varios residuos es un solo suceso mutacional, de modo que su coste no debería crecer linealmente con su longitud: un coste de apertura más un pequeño coste de extensión lo modela. Una penalización de apertura demasiado alta fuerza discrepancias donde corresponde un hueco y desalinea todo lo que sigue a una inserción real; una demasiado baja esparce huecos por todas partes, empareja residuos por azar e infla la identidad.

Ejercicio 5.8 ★★

Una búsqueda BLAST de una proteína humana contra una base de datos de mosca da un mejor acierto con $E = 10^{-30}$ que cubre los residuos 50–180 de la consulta de $600$ residuos. ¿Es la proteína de la mosca el ortólogo de la humana? ¿Qué prueba adicional se haría?

Solución

Solución de Ejercicio 5.8.

No necesariamente: el alineamiento cubre un segmento de $130$ residuos, que es la firma de un dominio compartido más que la de un ortólogo alineado en toda su longitud. Prueba: búsquese la proteína de la mosca contra el proteoma humano (¿es la consulta su mejor acierto, en toda la longitud?), identifíquese el dominio con un HMM de perfil y constrúyase un árbol génico de la familia en varias especies.

Ejercicio 5.9 ★★

¿Por qué los métodos de perfil detectan homólogos que el alineamiento por pares pierde? Dese un ejemplo de patrón de columna que un perfil capta y una secuencia aislada no.

Solución

Solución de Ejercicio 5.9.

Un perfil registra, columna a columna, lo que la familia tolera: una posición que es siempre hidrófoba pero nunca el mismo residuo, un residuo catalítico invariante, una posición que es siempre un hueco en la mitad de la familia. Un alineamiento por pares puntúa cada residuo contra otro único residuo y no puede saber que una valina en la posición 40 es “tan buena como” la isoleucina que hay ahí en la consulta. El perfil además pondera las columnas conservadas, de modo que una similitud débil concentrada allí donde la familia está conservada se vuelve significativa.

Ejercicio 5.10 ★★★

Muéstrese que, bajo un esquema de puntuación de esperanza positiva, el alineamiento local de Smith–Waterman de dos secuencias aleatorias largas tiene una puntuación que crece linealmente con su longitud, y explíquese por qué esto hace fracasar la teoría del valor E. ¿Qué implica para alinear ADN con coincidencia $+1$ y discrepancia $-1$ a un $60\,\%$ de contenido de GC?

Solución

Solución de Ejercicio 5.10.

Con una puntuación esperada positiva $\mu > 0$ por columna, la puntuación acumulada a lo largo de la diagonal de dos secuencias aleatorias es un paseo aleatorio con deriva positiva: tras $n$ columnas vale unos $\mu n$ , de modo que el mejor alineamiento local es esencialmente el todo y su puntuación crece como $\mu n$ y no como $\log n$ . La teoría de Karlin–Altschul, que exige una deriva negativa para que las puntuaciones altas sean excursiones raras, no se aplica y no existe ningún $\lambda$ . Para ADN al $60\,\%$ de GC la probabilidad de una coincidencia es $2(0.3^{2}) + 2(0.2^{2}) = 0.26$ , de modo que la puntuación esperada es $0.26 - 0.74 = -0.48$ : sigue siendo negativa y la estadística vale; pero un esquema como coincidencia $+1$ y discrepancia $-0.3$ tendría esperanza $+0.04$ y daría todo el genoma como un solo alineamiento.

Ejercicio 5.11 ★★★

La tabla de Needleman–Wunsch necesita $mn$ celdas de memoria; para dos cromosomas de $100\,\mathrm{Mb}$ son $10^{16}$ . Descríbanse dos ideas con las que los alineadores de genomas lo evitan (semillas y encadenamiento; bandas), y a qué renuncia cada una.

Solución

Solución de Ejercicio 5.11.

Semillas y encadenamiento: se buscan coincidencias exactas o casi exactas de $k$ -meros entre las dos secuencias con una tabla de dispersión, se conservan las que se alinean en diagonales coherentes, se encadenan y se aplica programación dinámica solo en los huecos entre semillas encadenadas; se renuncia a los alineamientos en regiones sin semilla (tramos muy divergentes). Bandas: si se sabe que las dos secuencias son casi colineales, se calculan solo las celdas de una banda de anchura $w$ alrededor de la diagonal, a un coste $wn$ en lugar de $mn$ ; se renuncia a cualquier alineamiento con una inserción mayor que la banda.

Ejercicio 5.12 ★★★

Un modelo oculto de Markov de búsqueda de genes para bacterias tiene estados para las tres posiciones del codón y para el ADN no codificante. Explíquese cómo el modelo puede distinguir la secuencia codificante de la no codificante sin información alguna sobre codones de parada (considérese el uso de codones), y por qué el mismo enfoque es mucho más difícil en un genoma humano.

Solución

Solución de Ejercicio 5.12.

La secuencia codificante tiene un período de tres: las tres posiciones del codón tienen composiciones de bases distintas (la tercera es la más sesgada), y el uso de codones es desigual en cada especie. Un modelo con tres estados codificantes en serie, cada uno emitiendo bases con la composición de esa posición del codón, asigna al ADN codificante una probabilidad mayor que el estado no codificante, a lo largo de una ventana de unas pocas decenas de codones, incluso sin paradas. En un genoma humano los exones son cortos ( $150\,\mathrm{bp}$ ) y están separados por intrones de kilobases, de modo que la señal codificante es breve e interrumpida; el modelo debe reconocer además los sitios de corte, que son señales débiles, y la enorme cantidad de secuencia no codificante produce muchos segmentos codificantes falsos.

5.7 Problema: una secuencia del fondo del mar

Problema 5.1

Problema de fin de semana — una proteína desconocida alineada a mano, buscada en las bases de datos con su significación calculada, su motivo regulador pesado en bits y su gen contrastado con la estadística de los marcos de lectura abiertos aleatorios, hasta llegar al valor E del mejor acierto, a los bits que necesita un sitio y a la longitud que ha de tener un marco de lectura para creérselo

Datos: una proteína de $300$ residuos de un anélido abisal. Base de datos de proteínas: $1.2\times 10^{11}$ residuos. Genoma del gusano: $1.6\,\mathrm{Gb}$ , $38\,\%$ de GC. Puntuación de los alineamientos a mano: coincidencia $+1$ , discrepancia $-1$ , hueco $-1$ . Puntuación en bits del mejor acierto BLAST: $92$ ; del décimo acierto: $38$ .

Parte I — A mano.

Alinéense los péptidos KQT y KAQT con la recurrencia de Needleman–Wunsch: escríbase la tabla y dense el alineamiento óptimo y su puntuación.
Repítase con Smith–Waterman (local) para GATCAT contra ACAT: hállense el mejor alineamiento local y su puntuación.
¿Cuántas actualizaciones de celda cuesta un alineamiento global de la proteína de $300$ residuos contra una proteína de $450$ residuos? ¿Y contra toda la base de datos?
Si un ordenador ejecuta $10^{9}$ actualizaciones por segundo, ¿cuánto tarda el alineamiento de la pregunta 3 contra toda la base de datos? ¿Por qué se usa BLAST en su lugar?
Una puntuación de identidad de BLOSUM62 es $+4$ para la alanina ( $p_{A} = 0.074$ ) y $+11$ para el triptófano ( $p_{W} = 0.013$ ). Con $\lambda = 0.347$ (unidades de medio bit), calcúlense la frecuencia objetivo $q_{AA}$ y $q_{WW}$ y la razón $q/p^{2}$ de cada una. Interprétese.
Dos proteínas comparten el $24\,\%$ de identidad a lo largo de $250$ residuos. Dígase por qué la identidad por sí sola no puede zanjar aquí la homología y qué sí podría.

Parte II — La búsqueda.

Calcúlense $mn$ para la consulta contra la base de datos y $\log_{2}(mn)$ .
Calcúlense el valor E del mejor acierto ( $S' = 92$ ) y el del décimo acierto ( $S' = 38$ ).
¿Qué puntuación en bits corresponde a $E = 10^{-3}$ en esta búsqueda? ¿Y a $E = 1$ ?
El mismo mejor acierto se encuentra cuando la base de datos ha crecido hasta $1.2\times 10^{12}$ residuos. ¿Su valor E?
El décimo acierto alinea los residuos 200–260 de la consulta con un $40\,\%$ de identidad a lo largo de $60$ residuos. Usando el valor E, dígase si es prueba de homología y qué podría añadir una búsqueda por perfil.
El mejor acierto es una quinasa humana, alineada a lo largo de los residuos 10–290. Su mejor acierto en el proteoma del gusano es la consulta. ¿Qué establece esta prueba recíproca, y qué no?

Parte III — Un motivo.

Aguas arriba del gen hay un sitio candidato de factor de transcripción de ocho posiciones con contenidos de información $2, 2, 1.6, 2, 0.8, 1.2, 0.4, 0.3$ bits. ¿ $R$ total?
¿Cuántas coincidencias por azar tiene el motivo en el genoma de $1.6\,\mathrm{Gb}$ (ambas hebras, $3.2\times 10^{9}$ posiciones)?
El factor regula unos $200$ genes. ¿Cuántos bits necesitaría un motivo para especificar $200$ sitios por sí solo en este genoma?
¿Cuánto del déficit podría aportar un segundo motivo contiguo de $8$ bits, si los dos han de aparecer juntos con un espaciado fijo?
Una posición con frecuencias $(0.5, 0.5, 0, 0)$ para (A, C, G, T): calcúlense su entropía y su contenido de información.
Explíquese, con el argumento de la información, por qué los factores de transcripción bacterianos suelen tener sitios más largos y más conservados que los eucariotas.

Parte IV — El gen mismo.

En ADN aleatorio de composición de bases uniforme, ¿cuál es la probabilidad de que un codón sea de parada? ¿Cuál es el número esperado de codones antes de que aparezca una parada (una distribución geométrica)?
El genoma del gusano tiene un $38\,\%$ de GC. Recalcúlense la probabilidad de que un codón al azar sea de parada (TAA, TAG, TGA) con las frecuencias de bases reales y la longitud esperada del marco de lectura. ¿En qué sentido empuja un contenido de GC bajo a la búsqueda de genes?
¿Cuál es la probabilidad de que un marco de lectura abierto al azar tenga al menos $100$ codones? ¿Y al menos $300$ ?
En el genoma de $1.6\,\mathrm{Gb}$ , seis marcos en dos hebras dan unos $3.2\times 10^{9}$ inicios de codón. ¿Cuántos marcos de lectura abiertos al azar de al menos $100$ codones se esperan? ¿Y de al menos $300$ ?
Explíquese por qué “marco de lectura abierto de más de $100$ codones” es un buscador de genes utilizable en una bacteria pero no en este genoma, y qué usa en su lugar un buscador de genes eucariota.
El gen del gusano tiene seis exones de $150\,\mathrm{bp}$ de media. Explíquese cómo las lecturas de secuenciación de ARN resuelven la estructura de exones que la secuencia genómica por sí sola deja ambigua.
Resúmase: el valor E del mejor acierto (pregunta 8), los bits necesarios para especificar $200$ sitios (pregunta 15) y el número esperado de marcos de lectura al azar de $300$ codones en el genoma (pregunta 22).

Solución

Solución de Problema 5.1.

1. Filas K, Q, T; columnas K, A, Q, T; bordes $0,-1,-2,-3,-4$ y $0,-1,-2,-3$ . Fila K: $1, 0, -1, -2$ ; fila Q: $0, 0, 1, 0$ ; fila T: $-1, -1, 0, 2$ . Óptimo $2$ : K-QT sobre KAQT. 2. Mejor puntuación local $3$ : CAT contra CAT (residuos 4–6 de GATCAT con 2–4 de ACAT); ATCAT contra A-CAT también puntúa $4 - 1 = 3$ . 3. $300\times 450 = 1.35\times 10^{5}$ actualizaciones; contra la base de datos, $300\times 1.2\times 10^{11} = 3.6\times 10^{13}$ . 4. $3.6\times 10^{4}$ s, diez horas por consulta; las semillas de BLAST se saltan casi toda la tabla y responden en segundos. 5. $q_{ab} = p_{a}p_{b}e^{\lambda s}$ . Alanina: $e^{1.39} = 4.0$ , $q_{AA} = 0.074^{2}\times 4.0 = 0.022$ , razón $4$ . Triptófano: $e^{3.82} = 45$ , $q_{WW} = 0.013^{2}\times 45 = 0.0077$ , razón $45$ . Un par de triptófanos alineado es $45$ veces más frecuente en homólogos que por azar, y un par de alaninas solo cuatro veces; aun así los pares de alaninas son más frecuentes en términos absolutos porque la alanina es frecuente. 6. El $24\,\%$ cae en la zona crepuscular, donde los alineamientos aleatorios alcanzan el $15\text{ a }20\,\%$ ; lo zanjarían el valor E del alineamiento, unos motivos conservados en las posiciones correctas, una coincidencia con un HMM de perfil de una familia conocida o un plegamiento compartido. 7. $mn = 300\times 1.2\times 10^{11} = 3.6\times 10^{13}$ ; $\log_{2}(mn) = 45.0$ . 8. $E(92) = 2^{45 - 92} = 2^{-47} \approx 7\times 10^{-15}$ ; $E(38) = 2^{7} = 128$ . 9. $E = 10^{-3}$ con $S' = 45 + 10 = 55$ bits; $E = 1$ con $45$ bits. 10. Diez veces $mn$ : $E \approx 7\times 10^{-14}$ , aún abrumador. 11. Con $E = 128$ , el décimo acierto es lo que produce el azar; un $40\,\%$ de identidad a lo largo de $60$ residuos no es prueba. Una búsqueda por perfil de los residuos 200–260 contra la base de datos de dominios podría mostrar si ese segmento es un dominio conocido, con una estadística de la que carece la comparación por pares. 12. Los mejores aciertos recíprocos en toda la longitud son compatibles con una ortología uno a uno; no la demuestran — una duplicación en un linaje posterior a la separación da dos coortólogos, y la pérdida del ortólogo verdadero puede dejar un parálogo como mejor acierto. La prueba es un árbol génico con varias especies. 13. $R = 2 + 2 + 1.6 + 2 + 0.8 + 1.2 + 0.4 + 0.3 = 10.3$ bits. 14. $3.2\times 10^{9}\times 2^{-10.3} \approx 2.5\times 10^{6}$ coincidencias por azar. 15. $\log_{2}(3.2\times 10^{9}/200) = \log_{2}(1.6\times 10^{7}) \approx 24$ bits. 16. La coaparición con un espaciado fijo suma los bits: $10.3 + 8 = 18.3$ , lo que aporta $8$ de los $13.7$ que faltan; unos $5.7$ bits (un factor de $50$ en las coincidencias por azar) han de venir de otro sitio — la accesibilidad de la cromatina, más compañeros. 17. $H = -(0.5\log_{2}0.5 + 0.5\log_{2}0.5) = 1$ bit; $R = 2 - 1 = 1$ bit. 18. Un factor bacteriano ha de encontrar sus pocos sitios en un genoma de $4.6\,\mathrm{Mb}$ sin ayuda de la cromatina: necesita unos $19$ bits, y sus sitios son largos y conservados. Un genoma eucariota es mil veces mayor y exige diez bits más, y sin embargo sus factores tienen sitios cortos; consiguen la especificidad por combinación y por la restricción a la cromatina accesible, lo que además hace la regulación más evolucionable, ya que un sitio corto se gana o se pierde con facilidad. 19. $3/64 = 0.047$ ; el número esperado de codones antes de una parada es $64/3 \approx 21$ . 20. $p_{A} = p_{T} = 0.31$ , $p_{G} = p_{C} = 0.19$ : $P(\text{TAA}) = 0.31^{3} = 0.030$ , $P(\text{TAG}) = P(\text{TGA}) = 0.31^{2}\times 0.19 = 0.018$ ; total $0.066$ , longitud esperada del marco, $15$ codones. El ADN rico en AT está lleno de paradas, de modo que los marcos abiertos al azar son más cortos y los largos destacan más. 21. $(61/64)^{100} = e^{-4.80} = 0.008$ ; $(61/64)^{300} = e^{-14.4} = 5.6\times 10^{-7}$ . 22. Cada marco abierto maximal termina en una parada, y $3.2\times 10^{9}$ inicios de codón contienen $3.2\times 10^{9}\times 3/64 = 1.5\times 10^{8}$ paradas: unos $1.5\times 10^{8}\times 0.008 = 1.2\times 10^{6}$ marcos al azar de al menos $100$ codones, y $1.5\times 10^{8}\times 5.6\times 10^{-7} \approx 80$ de al menos $300$ . 23. Una bacteria de $4.6\,\mathrm{Mb}$ tiene unas $4\times 10^{5}$ paradas y por tanto unos $3500$ marcos por azar de $100$ codones pero casi ninguno de $300$ ; sus genes promedian $300$ codones y el $88\,\%$ del ADN es codificante, de modo que un marco abierto largo es casi siempre un gen. En el gusano codifica el $1.5\,\%$ del ADN, los exones promedian $50$ codones — menos que el umbral del azar — y un millón de marcos al azar de $100$ codones los sepultan. Los buscadores de genes eucariotas usan señales de sitios de corte, el sesgo de codones en un modelo oculto de Markov, la homología con proteínas conocidas y, sobre todo, los transcritos secuenciados. 24. Una lectura de un mensajero cortado y empalmado se alinea con el genoma en dos trozos separados por un intrón: la partición marca los dos sitios de corte hasta la base; la cobertura de lecturas delimita los exones y las lecturas emparejadas enlazan exones sucesivos en un solo transcrito, lo que resuelve cuál de varios sitios de corte candidatos se usa. 25. $E \approx 7\times 10^{-15}$ para el mejor acierto; unos $24$ bits para especificar $200$ sitios en el genoma; y unos $80$ marcos de lectura por azar de $300$ codones en todo el genoma.

Términos definidos en este capítulo

Ver los 479 términos del glosario