Aplicación del modelo de programación CUDA en la simulación de la evolución de secuencias genéticas

Procedimientos de simulación basados en la técnica de Markov

Para la realización de simulaciones de la evolución molecular de secuencias genéticas han sido utilizados diferentes métodos de Markov (Weber, 2012), (Link & Eaton, 2012), dadas las capacidades de este tipo de simulación para representar procesos estocásticos y la divergencia genética a lo largo del tiempo. Todos estos métodos clásicos (Yang, 2006) requieren como entradas un árbol filogenético con una longitud de una rama específica, una secuencia genética de raíz, y un modelo Markoviano de la evolución molecular. La simulación es realizada por cada rama del árbol filogenético, empezando por la raíz y terminando en las puntas de las hojas.

Tiempo de espera evolutivo

En (Counsell, 2005) se plantea que un proceso de Markov de tiempo continuo, determina el tiempo de espera evolutivo entre las sustituciones, o sea, el período que el sitio se mantiene en el estado i, y que este es distribuido exponencialmente con el parámetro |q*ii|. Dado que este último es la tasa de sustitución de un nucleótido determinado i, es fácil construir una variable aleatoria de distribución exponencial a partir de la otra uniformemente distribuida entre 0 y 1 (u) utilizando la ecuación 1:

(1)

Sea T el conjunto representado por los tiempos transcurridos para cada sitio de la secuencia, durante cada paso de la simulación, y t el tiempo de espera generado utilizando dicho número aleatorio. Si se cumple que t > T entonces se mantendrá el sitio en su estado actual, de lo contrario, ocurre una sustitución nucleotídica. Cuando t < T, o sea, que va a ocurrir una mutación, el nucleótido de remplazo es seleccionado de acuerdo con la cadena de Markov de tiempo discreto D^* = {d^* _ij} definida por la ecuación 2:

(2)

El nucleótido sustituto, se obtiene igualmente generando un número aleatorio, distribuido uniformemente entre 0 y 1. Si este número aleatorio es menor que d, el nucleótido 1 sustituye al nucleótido i. De lo contrario, si el número aleatorio es menor que d*i1 + d*i2, el nucleótido 2 sustituye al i. Este proceso se repite hasta que se encuentre el nucleótido de sustitución. Por el hecho que se cumple que: , siempre se encontrará un nucleótido de sustitución, y dado que , el nucleótido sustituido no será el mismo.

Cadena de salto

En este método se ignoran los tiempos de espera (exponenciales con media 1/|q*ii|) entre transiciones, la secuencia de los estados constituirá una cadena de Markov de tiempo discreto, con matriz de transición que se muestra en la ecuación 3:

(3)

Tanto el método de los tiempos de espera evolutivo, como la cadena de salto, pueden ser aplicados a toda la secuencia (todos los sitios en lugar de uno), con tasa de sustitución igual a la sumatoria de las tasas en todos los sitios, y con un tiempo de espera hasta que una sustitución ocurra en cualquier sitio de toda la secuencia, que sigue una distribución exponencial de media igual a la inversa de la tasa total, para cada uno de los sitios de la secuencia.

Metropolis-Hasting (MH)

Este último método es “conocido como el más popular de la familia de métodos Markov Chain Monte Carlo, y comúnmente usado para la evolución de secuencias” (Andrieu & Doucet, 2010; Cornebise & Peters, 2009). MH consiste básicamente en generar un nuevo valor candidato x* dado el valor actual x, de acuerdo con la probabilidad de distribución propuesta q(x*|x). La probabilidad de aceptación que la cadena de Markov se mueva hacia x* está dada por la ecuación 4, en el cual p(x) es la distribución invariante definida en x:

(4)

Existen varios algoritmos prácticos interpretados como casos especiales o extensiones de este método, cuya principal diferencia es la especificación de la probabilidad de distribución propuesta q(x*|x). Es necesario destacar, que las probabilidades de distribución y mutación deben ser determinadas para un nuevo estado, siempre que ocurra una mutación de nucleótido, o que el valor candidato x* sea aceptado acorde con la probabilidad previamente definida para el nuevo nucleótido. Por el contrario, si la cadena de Markov se mantiene en el estado actual con el mismo valor nucleotídico, no existirá la necesidad de cambiar las probabilidades de distribución y mutación de la base en ese sitio.

Los métodos de tiempo de espera evolutivo y de cadena de salto poseen la ventaja de no requerir el cálculo de la matriz de probabilidades P(t) sobre el tiempo t requerida por MH, puesto que, tanto la matriz de transición para la cadena de salto como los tiempos de espera, se especifican en su totalidad por la matriz de tasas instantáneas Q* definida por el modelo evolutivo Markoviano seleccionado para la simulación de la evolución molecular de las secuencias. Los tres métodos descritos con anterioridad no son los únicos posibles, de hecho, estos tres procedimientos en la práctica son comúnmente combinados dando lugar a nuevos.

Simulación basada en el método de MH y el tiempo de espera

Mediante la estimación del tiempo de espera evolutivo se evalúa la ocurrencia o no de una mutación en un sitio determinado. Dado el tiempo evolutivo T transcurrido en cada paso de la simulación, si t < T, o sea, que hay un cambio de nucleótido en el sitio dado, se aplica el método de MH de forma repetitiva hasta que un nuevo nucleótido de sustitución sea encontrado.

Simulación basada en el tiempo de espera y la matriz de transición

Mediante el proceso de simulación del tiempo de espera se determinan los nucleótidos que cambiarán, de la misma forma que en el procedimiento anterior. En caso que haya cambio de nucleótido, se buscará el nucleótido de remplazo mediante la matriz de transición de la cadena de salto, de la siguiente manera: al suponer que i es el nucleótido de partida ocupado en el sitio dado (el estado actual del proceso es i), a partir de la matriz de tasas de sustitución Q*, se obtiene la fila correspondiente a ese nucleótido. Los valores en esta fila constituyen una cadena de Markov de tiempo discreto D^* = {d^* _ij}, cuya definición fue descrita con anterioridad.

Simulación basada en el método de MH con repetición

Mediante el uso de MH se determinará el cambio de nucleótido en cada sitio de la secuencia. Posteriormente, mediante la matriz de probabilidades de sustitución P(t), se obtendrá la fila correspondiente al nucleótido ocupado en el sitio dado, esta fila es el vector de distribución de probabilidades de ese nucleótido y se denota por p. Para saber si ocurre la mutación del nucleótido actual i, se genera un número aleatorio con distribución uniforme entre 0 y 1 llamado u. Si se cumple que u < p(i), no ocurrirá cambio en el sitio dado, en caso contrario, se buscará un nucleótido de sustitución diferente al nucleótido actual mediante la aplicación repetitiva del método MH, garantizando el cambio del proceso a un nuevo estado.

Simulación basada en el método de MH sin repetición

Este procedimiento no examina si una mutación ocurre o no en un sitio determinado. El método de MH se aplica una sola vez para cada sitio de la secuencia, y el nucleótido puede cambiar su valor a otro diferente o mantenerse en su estado actual, según el valor devuelto por MH. Si MH retorna una base igual que la actual en el sitio dado, se dice que no hay mutación. De otro modo, el nucleótido ancestral se remplaza por el nuevo.

Modelos Markovianos de la evolución molecular

Para estimar el número real de sustituciones entre secuencias relativamente alejadas en el proceso evolutivo, se necesita un modelo probabilístico que describa los cambios entre nucleótidos. Las cadenas continuas de Markov son la herramienta más usada con este propósito (Yang, 2006), donde normalmente se asume que cada sitio de la secuencia evoluciona de forma independiente. Las sustituciones en cada sitio en particular son descritas mediante una cadena de Markov, que tiene como estados las bases nucleotídicas. La imposición de algunas restricciones en las tasas de sustitución entre nucleótidos lleva a diferentes modelos. En la literatura (Yang & Rodríguez, 2013), han sido diversos los modelos de Markov propuestos para la evolución de las secuencias, y se han desarrollado también, varias técnicas para estimar experimentalmente los parámetros de dichos modelos de Markov mediante el análisis de secuencias reales observadas.

Según Yang (Yang, 2006), de forma general, los modelos evolutivos describen el modo y la probabilidad que una secuencia de nucleótidos cambie a otra secuencia de nucleótidos homóloga a lo largo del tiempo. Estos modelos describen para cada uno de los sitios de la matriz, la probabilidad que se produzca el cambio de un nucleótido a otro a lo largo de las ramas de un árbol filogenético dado.

Los modelos de evolución de nucleótidos se definen matemáticamente mediante dos clases de parámetros que determinan el cambio:

Frecuencia de cada nucleótido: parámetro que mide la frecuencia de cada nucleótido en la matriz de datos. En los modelos evolutivos más sencillos, toma una misma frecuencia para los cuatro nucleótidos (πA = πC = πG = πT = 0.25), sin tener en cuenta la frecuencia de aparición de los mismos en la matriz de datos; En el caso de los modelos más complejos asume que las frecuencias asociadas a cada uno de los nucleótidos son diferentes y son calculadas a partir de los datos (πA ≠ πC ≠ πG ≠πT).

Tipos de sustituciones y sus correspondientes tasas de sustitución (rate parameters): las tasas de sustitución se representan con las tasas relativas de cambio de un nucleótido a otro para una posición de un tiempo t0 a un tiempo t1. Así, cada posición de la matriz de tasas de sustitución tendrá una probabilidad asociada de cambio para cada unidad de tiempo (unidad de distancia evolutiva). Los modelos más sencillos asumen una misma tasa relativa para todas las sustituciones posibles, mientras que los más complicados asumen una tasa relativa diferente para cada tipo de sustitución. A partir de estas tasas relativas se calcula la tasa media de sustitución (μ). La matriz de tasas más general de un modelo de Markov tiene solo el requisito de que los elementos no diagonales son positivos y la suma de cada fila es cero.

Modelo de programación en CUDA

La programación en CUDA está basada en tres abstracciones básicas: una jerarquía de grupos de hilos, otra de tipos de memoria y barreras de sincronización. La estructura que conforma la jerarquía de hilos que se ejecutan en el dispositivo NVIDIA (device) se agrupan, en tres elementos: mallas, bloques e hilos. Una malla está conformada por múltiples; un bloque está formado por un grupo de hilos, y el hilo es la unidad elemental, a continuación, en la Figura 1 se muestra esta jerarquía.

Como el dispositivo NVIDIA (tarjeta gráfica) está colocado en una computadora se establece un modelo “híbrido” para la programación en el que se combina el código del programa principal que se ejecuta en la CPU de la computadora anfitrión (host) y las llamadas a funciones que se ejecutan en el dispositivo NVIDIA (device). El proceso comienza en un programa que se ejecuta en la CPU anfitriona y este programa a su vez invoca funciones que se ejecutarán paralelamente por varios hilos en la GPU, conocidas como funciones kernel, especificando para cada una de ellas la configuración de ejecución que define la dimensión de la malla de bloques (dimGrid) y la cantidad de hilos por bloques (dimBlock) con que será ejecutada la función kernel en la GPU.

Fig. 1:

Jerarquía del lote de hilos (Tomado de (Nvidia, 2011))

Los hilos que se ejecutan en un dispositivo CUDA tienen acceso a múltiples espacios de memoria: global, local, compartida, textura y registros. Cada hilo tiene un espacio de memoria local privado; a su vez, cada bloque de hilos posee memoria compartida visible solo a todos los hilos del bloque y con la misma duración de vida que el bloque. Todos los hilos ejecutados por una función kernel tienen acceso a la misma memoria global. Los espacios de memoria global, constante y de textura son persistentes a diferentes activaciones de funciones kernels en la misma aplicación. El trabajo de estos hilos puede ser sincronizado mediante directivas de bloqueo, lo cual posibilita la coordinación entre estos; sin embargo, hilos agrupados en diferentes bloques no pueden comunicarse entre sí.

La estructura general de una aplicación CUDA incluye los pasos siguientes:

El computador (CPU) anfitrión ejecuta el cuerpo principal del programa (main())

Se reserva memoria dentro del dispositivo GPU

Se copian los datos del CPU anfitrión al dispositivo GPU

El CPU anfitrión llama a la función kernel del GPU

El dispositivo GPU ejecuta el código de manera paralela

Se copian los resultados de vuelta a la memoria del CPU anfitrión

Se libera la memoria reservada dentro del dispositivo GPU

Repetir los pasos del 2 al 7 tantas veces como sean necesarios para la solución del problema

Debido a los tiempos relativamente grandes de retardo (latencia) y al bajo ancho de banda en las transferencias de memoria entre la computadora anfitrión y el dispositivo GPU, es recomendable dividir la aplicación, de tal manera que cada componente de hardware haga únicamente el trabajo que mejor realiza. El uso del GPU es solamente recomendado si (Nvidia, 2011), (Sánchez, Carbajal, Cortés, & Fernández, 2012):

La complejidad de las operaciones a realizar en la GPU justifica el costo de mover datos desde y hacia el dispositivo GPU. Se deben minimizar las transferencias y mantener en GPU los datos tanto tiempo como sea posible.

La aplicación realizará una misma operación sobre numerosos datos al mismo tiempo en paralelo, mediante operaciones simples que son asignadas a múltiples hilos.

Los tipos de variables y arreglos utilizados, y su tamaño, son congruentes con los patrones de memoria y las instrucciones aritméticas de la GPU, garantizando su buen desempeño.

Tomando en cuenta estas consideraciones, resulta sumamente importante determinar que partes del proceso de simulación que se desea implementar usando CUDA deben ser paralelizadas y cuáles no.

Proceso general para simular la evolución de secuencias

A modo general, el proceso de simulación de la evolución molecular de secuencias genéticas se realiza mediante los pasos algorítmicos descritos a continuación:

Paso 1: Especificar la cantidad de descendientes que pudiera tener cada ancestro.

Paso 2: Generar una secuencia mutante aplicando la técnica de Markov seleccionada, a partir de la secuencia ancestral.

Paso 3: Eliminar los codones de parada existentes en las secuencias mutantes obtenidas.

Paso 4: Clasificar la secuencia recién creada a través de la secuencia de aminoácidos sintetizada. Si coincide con uno de los alelos de la población, se aumentará el contador de este en 1; de lo contrario, se agregará la mutante a la población como una nueva clase. En ambos casos, la cantidad total de secuencias de la población se incrementa en 1.

Paso 5: Calcular la distribución de probabilidad de la población actual.

Paso 6: Tomar la secuencia mutante generada como el nuevo ancestro y volver el algoritmo al paso 2 si la cantidad de descendientes no se ha acabado; de lo contrario, moverse al paso siguiente.

Paso 7: Seleccionar el nuevo ancestro entre las secuencias de la población actual, si el número de ancestros no se ha terminado, el algoritmo vuelve al paso 1; de otro modo, ir al paso siguiente.

Paso 8: Terminar el algoritmo.

Al realizar un análisis de este proceso y de las consideraciones abordadas con anterioridad para el uso de CUDA, se identificaron como procedimientos a ser paralelizados la mutación (Paso 2), y la clasificación de las secuencias (Paso 4). Acorde a esta valoración el esquema de paralelización propuesto para la simulación de la evolución de secuencias es el que se muestra en la Figura 2.

Paralelización de las mutaciones

Los sitios de una secuencia genética mutan de forma independiente, es decir, que las mutaciones que ocurran en uno de ellos no influyen en las que puedan ocurrir en otros sitios de la secuencia.

Fig. 2:

Esquema general de paralelización de la simulación.

Esta independencia facilita notablemente el proceso de paralelización; y dadas sus características se aplicó la técnica SIMD (Single Instructions Multiple Data), que permite conseguir paralelismo a nivel de datos y consiste en aplicar una misma operación sobre un conjunto de datos, en nuestro caso secuencias genómicas, y con una organización en donde una única unidad común de control despacha las instrucciones a diferentes unidades de procesamiento. Todas estas reciben la misma instrucción, pero operan sobre diferentes conjuntos de datos.

Como el proceso de mutar un sitio es invariable y solo interesa la base nitrogenada del ancestro, entonces se programa dicho procedimiento como una operación y al aplicarlo a más de un sitio durante la ejecución se obtiene la paralelización. En la implementación en CUDA de este esquema se realiza un kernel para el proceso de mutar una generación de secuencias. Las secuencias se modelan como arreglos de tipo char, pues solo se cuenta con 4 bases. Esta variante posee un menor consumo en memoria, pero el procesamiento de un sitio necesita de más operaciones debido a la necesaria conversión de la base a una codificación numérica. En una implementación secuencial, si necesitamos obtener 10 descendientes el proceso debe de ejecutarse 10 veces, pues se obtiene solamente una mutación en cada iteración del algoritmo. Por su parte la implementación paralela del algoritmo logra esta tarea en un solo paso, generando de una vez todas las mutaciones deseadas para cada uno de los ancestros de una generación, en la Figura 3 se muestra gráficamente el proceso de mutación paralelizado.

La sincronización de los hilos se logra utilizando _syncthreads(), una función propia de CUDA; esta sincronización resulta necesaria para el cálculo de los elementos del procedimiento evolutivo seleccionado (matriz Q* de tasas de sustitución, matriz P(t) de probabilidades de transición, o matriz

Fig. 3:

Proceso de mutación paralelizado por sitio.

de transición M) que en el esquema de paralelización propuesto se realiza en el hilo con índice 0 de cada bloque, y por tanto los demás hilos del bloque deben esperar a que este hilo coloque su resultado en memoria compartida para poder comenzar el proceso mutacional. Aunque esta forma de obtener los elementos del procedimiento seleccionado contiene divergencia entre los hilos de un bloque, la rapidez con la que se producen los accesos a la memoria compartida por parte de los hilos conduce a un menor consumo temporal.

Algoritmo paralelo para realizar las mutaciones en cada hilo

Entrada: arreglo de ancestros A, arreglo para almacenar las mutaciones M, número de ancestros nag, número de descendiente por cada ancestro ndv, número de sitios N, y otros elementos necesarios según procedimiento implementado (ej. frecuencias de aparición F, alpha1, alpha2, beta, etc.).

Salida: Mutaciones almacenadas en R.

Paralelización de la clasificación de las mutaciones mediante sus proteínas

El proceso de clasificación, se llevaba a cabo con el objetivo de detectar cuando dos secuencias distintas sintetizan una misma proteína. De forma secuencial se hace comparando cada secuencia obtenida con cada una de las secuencias almacenadas, tomando en cuenta tanto las mutaciones generadas como las secuencias iniciales. En la versión paralela del algoritmo se han utilizado algunas estrategias en función de lograr mayor nivel de paralelización durante la ejecución paralela del algoritmo siguiente:

Algoritmo paralelo para clasificar las mutaciones

Entrada: aminoácidos de las mutaciones R, aminoácidos de los ancestros A, índices I, num_sitios, nag, ndv

Salida: Arreglo de índices (I) con valores de las clasificaciones.

Para lograr una adecuada paralelización de este procedimiento de clasificación se ha utilizado un arreglo de índices (I), el cual tendrá una posición por cada una de las secuencias obtenidas en una generación, por lo cual su longitud será igual a (nag * ndv), estará inicializado con −100 en todas sus posiciones y se almacenará en la memoria global. En los índices de este arreglo se tendrá −100 en caso de que la correspondiente secuencia no sintetice una proteína obtenida anteriormente, un valor del intervalo [0, nag − 1] si la misma produce una proteína igual que algunas de las correspondientes a los ancestros o el mismo estará en el intervalo [nag, nag + nag * ndv −1] si la proteína se corresponde con algunas de las mutaciones que se generaron junto con ella en la actual generación.

En el algoritmo paralelo se asigna un hilo a cada una de las secuencias, diferente a como se paralelizaron las mutaciones, en las cuales un hilo era responsable de la mutación de un sitio dentro de la secuencia. Luego que el hilo conoce que mutación le corresponde entonces realiza la comparación entre su respectiva secuencia y todos los ancestros, si la misma no es clasificada por ningún ancestro se procede a compararla con cada una de las mutaciones que se obtuvieron de forma posterior a ella, buscando las que ella clasifica. En el proceso se actualiza la información contenida en el arreglo, como este se encuentra en la memoria compartida y es escrito y leído por los hilos de forma simultánea, entonces constituye una sección crítica y se hace uso de un semáforo para realizar tanto las lecturas como las escrituras de forma atómica.