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