Una guía paso a paso para el análisis de datos de secuenciación de ADN

El Dr. Javier Quilez Oliete, un experimentado consultor de bioinformática en Kolabtree, proporciona una guía completa para el análisis de datos de secuenciación de ADN, incluyendo las herramientas y el software utilizado para leer los datos.

Introducción

El ácido desoxirribonucleico (ADN) es la molécula que transporta la mayor parte de la información genética de un organismo. (En algunos tipos de virus, la información genética es transportada por el ácido ribonucleico (ARN)). Los nucleótidos (representados convencionalmente por las letras A, C, G o T) son las unidades básicas de las moléculas de ADN. Conceptualmente, la secuenciación del ADN es el proceso de lectura de los nucleótidos que componen una molécula de ADN (por ejemplo, «GCAAACCAAT» es una cadena de 10 nucleótidos de ADN). Las tecnologías actuales de secuenciación producen millones de estas lecturas de ADN en un tiempo razonable y a un coste relativamente bajo. Como referencia, el coste de secuenciar un genoma humano -un genoma es el conjunto completo de moléculas de ADN de un organismo- ha bajado la barrera de los 100 dólares y puede hacerse en cuestión de días. Esto contrasta con la primera iniciativa para secuenciar el genoma humano, que se completó en una década y tuvo un coste de unos 2.700 millones de dólares.

Esta capacidad para secuenciar el ADN con un alto rendimiento y bajo coste ha permitido el desarrollo de un número creciente de métodos y aplicaciones basados en la secuenciación. Por ejemplo, la secuenciación de genomas enteros o de sus regiones codificadoras de proteínas (dos enfoques conocidos, respectivamente, como secuenciación del genoma completo y del exoma) en personas enfermas y sanas puede indicar las alteraciones del ADN causantes de enfermedades. Asimismo, la secuenciación del ARN que se transcribe a partir del ADN -técnica conocida como secuenciación del ARN- se utiliza para cuantificar la actividad de los genes y cómo ésta cambia en diferentes condiciones (por ejemplo, sin tratamiento frente a con tratamiento). Por otro lado, los métodos de secuenciación de captura de la conformación cromosómica detectan las interacciones entre moléculas de ADN cercanas y, por tanto, ayudan a determinar la distribución espacial de los cromosomas dentro de la célula.

Común a estas y otras aplicaciones de la secuenciación del ADN es la generación de conjuntos de datos del orden de los gigabytes y que comprenden millones de secuencias de lectura. Por lo tanto, para dar sentido a los experimentos de secuenciación de alto rendimiento (HTS) se necesitan importantes capacidades de análisis de datos. Afortunadamente, existen herramientas computacionales y estadísticas específicas y flujos de trabajo de análisis relativamente estándar para la mayoría de los tipos de datos de HTS. Aunque algunos de los pasos (iniciales) del análisis son comunes a la mayoría de los tipos de datos de secuenciación, los análisis posteriores dependerán del tipo de datos y/o del objetivo final del análisis. A continuación, proporciono una cartilla sobre los pasos fundamentales en el análisis de datos HTS y hago referencia a herramientas populares.

Algunas de las secciones siguientes se centran en el análisis de datos generados a partir de tecnologías de secuenciación de lectura corta (principalmente Illumina), ya que éstas han dominado históricamente el mercado de HTS. Sin embargo, las nuevas tecnologías que generan lecturas más largas (por ejemplo, Oxford Nanopore Technologies, PacBio) están ganando terreno rápidamente. Como la secuenciación de lecturas largas tiene algunas particularidades (por ejemplo, mayores tasas de error), se están desarrollando herramientas específicas para el análisis de este tipo de datos.

Control de calidad (QC) de las lecturas en bruto

El analista entusiasta comenzará el análisis a partir de archivos FASTQ; el formato FASTQ ha sido durante mucho tiempo el estándar para almacenar datos de secuenciación de lectura corta. En esencia, los archivos FASTQ contienen la secuencia de nucleótidos y la calidad de llamada por base para millones de lecturas. Aunque el tamaño del archivo dependerá del número real de lecturas, los archivos FASTQ suelen ser grandes (del orden de megabytes y gigabytes) y estar comprimidos. Cabe destacar que la mayoría de las herramientas que utilizan archivos FASTQ como entrada pueden manejarlos en formato comprimido, por lo que, para ahorrar espacio en el disco, se recomienda no descomprimirlos. Como convención, aquí equipararé un archivo FASTQ a una muestra de secuenciación.

FastQC es probablemente la herramienta más popular para llevar a cabo el QC de las lecturas crudas. Puede ejecutarse a través de una interfaz visual o mediante programación. Mientras que la primera opción puede ser más conveniente para los usuarios que no se sienten cómodos con el entorno de línea de comandos, la segunda ofrece una escalabilidad y reproducibilidad incomparables (piense en lo tedioso y propenso a errores que puede ser ejecutar manualmente la herramienta para decenas de archivos). En cualquier caso, el resultado principal de FastQC es un archivo HTML que informa de las estadísticas resumidas clave sobre la calidad general de las lecturas de secuenciación en bruto de una muestra determinada. Inspeccionar decenas de informes FastQC uno por uno es tedioso y complica la comparación entre muestras. Por lo tanto, es posible que desee utilizar MultiQC, que agrega los informes HTML de FastQC (así como de otras herramientas utilizadas posteriormente, por ejemplo, el recorte de adaptadores, la alineación) en un único informe.

La información sobre el QC está pensada para que el usuario pueda juzgar si las muestras tienen una buena calidad y, por lo tanto, pueden utilizarse para los pasos posteriores o deben descartarse. Lamentablemente, no existe un umbral consensuado basado en la métrica FastQC para clasificar las muestras como de buena o mala calidad. El enfoque que utilizo es el siguiente. Espero que todas las muestras que han pasado por el mismo procedimiento (por ejemplo, la extracción de ADN, la preparación de la biblioteca) tengan estadísticas de calidad similares y una mayoría de indicadores de «aprobado». Si algunas muestras tienen una calidad inferior a la media, las seguiré utilizando en el análisis posterior teniendo esto en cuenta. Por otro lado, si todas las muestras del experimento obtienen sistemáticamente banderas de «advertencia» o «suspenso» en múltiples métricas (véase este ejemplo), sospecho que algo ha ido mal en el experimento (por ejemplo, mala calidad del ADN, preparación de la biblioteca, etc.) y recomiendo repetirlo.

Recorte de lecturas

El QC de las lecturas en bruto ayuda a identificar las muestras problemáticas pero no mejora la calidad real de las lecturas. Para ello, necesitamos recortar las lecturas para eliminar las secuencias técnicas y los extremos de baja calidad.

Las secuencias técnicas son restos del procedimiento experimental (por ejemplo, adaptadores de secuenciación). Si estas secuencias son adyacentes a la verdadera secuencia de la lectura, el alineamiento (ver más abajo) puede asignar las lecturas a una posición incorrecta en el genoma o disminuir la confianza en un alineamiento dado. Además de las secuencias técnicas, también podemos querer eliminar las secuencias de origen biológico si éstas están muy presentes entre las lecturas. Por ejemplo, los procedimientos de preparación del ADN subóptimos pueden dejar una alta proporción de ARN ribosómico (ARNr) convertido en ADN en la muestra. A menos que este tipo de ácido nucleico sea el objetivo del experimento de secuenciación, mantener las lecturas derivadas del ARNr sólo aumentará la carga computacional de los pasos posteriores y puede confundir los resultados. Cabe destacar que si los niveles de secuencias técnicas, ARNr u otro contaminante son muy altos, lo que probablemente ya habrá sido puesto de manifiesto por el control de calidad, es posible que desee descartar toda la muestra de secuenciación.

En la secuenciación de lectura corta, la secuencia de ADN se determina un nucleótido a la vez (técnicamente, un nucleótido cada ciclo de secuenciación). En otras palabras, el número de ciclos de secuenciación determina la longitud de la lectura. Un problema conocido de los métodos de secuenciación HTS es la disminución de la precisión con la que se determinan los nucleótidos a medida que se acumulan los ciclos de secuenciación. Esto se refleja en una disminución general de la calidad de la llamada por base, especialmente hacia el final de la lectura. Al igual que ocurre con las secuencias técnicas, intentar alinear lecturas que contienen extremos de baja calidad puede llevar a una mala colocación o a una mala calidad de mapeo.

Para eliminar las secuencias técnicas/contaminantes y los extremos de baja calidad, existen herramientas de recorte de lecturas como Trimmomatic y Cutadapt, que se utilizan ampliamente. En esencia, estas herramientas eliminan las secuencias técnicas (disponibles internamente y/o proporcionadas por el usuario) y recortan las lecturas en función de la calidad, maximizando la longitud de la lectura. Las lecturas que quedan demasiado cortas después del recorte se descartan (las lecturas excesivamente cortas, por ejemplo, <36 nucleótidos, complican el paso de alineación, ya que es probable que se asignen a múltiples sitios en el genoma). Es posible que desee observar el porcentaje de lecturas que sobreviven al recorte, ya que una alta tasa de lecturas descartadas es probablemente un signo de datos de mala calidad.

Por último, suelo volver a ejecutar FastQC en las lecturas recortadas para comprobar que este paso ha sido eficaz y ha mejorado sistemáticamente las métricas de control de calidad.

Alineación

Salvo excepciones (por ejemplo, el ensamblaje de novo), la alineación (también denominada mapeo) suele ser el siguiente paso para la mayoría de los tipos de datos y aplicaciones de HTS. La alineación de la lectura consiste en determinar la posición en el genoma de la que se deriva la secuencia de la lectura (típicamente expresada como cromosoma:extremo inicial). Por lo tanto, en este paso se requiere el uso de una secuencia de referencia para alinear/mapear las lecturas.

La elección de la secuencia de referencia estará determinada por múltiples factores. Por un lado, la especie de la que procede el ADN secuenciado. Aunque el número de especies con una secuencia de referencia de alta calidad disponible es cada vez mayor, puede que todavía no sea el caso de algunos organismos menos estudiados. En esos casos, es posible que desee alinear las lecturas con una especie evolutivamente cercana para la que exista un genoma de referencia. Por ejemplo, como no existe una secuencia de referencia para el genoma del coyote, podemos utilizar la del perro, estrechamente relacionado, para la alineación de las lecturas. Del mismo modo, es posible que queramos alinear nuestras lecturas con una especie estrechamente relacionada para la que existe una secuencia de referencia de mayor calidad. Por ejemplo, aunque se ha publicado el genoma del gibón, éste está dividido en miles de fragmentos que no recapitulan completamente la organización de ese genoma en decenas de cromosomas; en ese caso, realizar el alineamiento utilizando la secuencia de referencia humana puede ser beneficioso.

Otro factor a considerar es la versión del ensamblaje de la secuencia de referencia, ya que se publican nuevas versiones a medida que la secuencia se actualiza y mejora. Es importante destacar que las coordenadas de un determinado alineamiento pueden variar entre versiones. Por ejemplo, se pueden encontrar múltiples versiones del genoma humano en el UCSC Genome Browser. En cualquier especie, soy muy partidario de migrar a la versión más reciente del ensamblaje una vez que se haya liberado por completo. Esto puede causar algunas molestias durante la transición, ya que los resultados ya existentes serán relativos a las versiones anteriores, pero vale la pena en el largo plazo.

Además, el tipo de datos de secuenciación también importa. Las lecturas generadas por los protocolos DNA-seq, ChIP-seq o Hi-C se alinearán con la secuencia de referencia del genoma. Por otro lado, como el ARN transcrito del ADN se procesa posteriormente en ARNm (es decir, se eliminan los intrones), muchas lecturas de ARN-seq no se alinearán con una secuencia de referencia del genoma. En su lugar, tenemos que alinearlas con las secuencias de referencia del transcriptoma o utilizar alineadores que tengan en cuenta la división (véase más adelante) cuando se utiliza la secuencia del genoma como referencia. En relación con esto está la elección de la fuente para la anotación de la secuencia de referencia, es decir, la base de datos con las coordenadas de los genes, transcritos, centrómeros, etc. Yo suelo utilizar la anotación de GENCODE, ya que combina la anotación completa de genes y secuencias de transcripción.

Se ha desarrollado una larga lista de herramientas de alineación de secuencias de lectura corta (véase la sección de alineación de secuencias de lectura corta aquí). Revisarlas va más allá del alcance de este artículo (los detalles sobre los algoritmos detrás de estas herramientas se pueden encontrar aquí). En mi experiencia, entre las más populares están Bowtie2, BWA, HISAT2, Minimap2, STAR y TopHat. Mi recomendación es que elija su alineador teniendo en cuenta factores clave como el tipo de datos HTS y la aplicación, así como la aceptación por parte de la comunidad, la calidad de la documentación y el número de usuarios. Por ejemplo, uno necesita alineadores como STAR o Bowtie2 que sean conscientes de las uniones exón-exón al mapear RNA-seq al genoma.

Común a la mayoría de los mapeadores es la necesidad de indexar la secuencia utilizada como referencia antes de que el alineamiento real tenga lugar. Este paso puede llevar mucho tiempo, pero sólo es necesario hacerlo una vez para cada secuencia de referencia. La mayoría de los mapeadores almacenan los alineamientos en archivos SAM/BAM, que siguen el formato SAM/BAM (los archivos BAM son versiones binarias de los archivos SAM). El alineamiento es uno de los pasos que más computaciones y tiempo consume en el análisis de los datos de secuenciación y los archivos SAM/BAM son pesados (del orden de los gigabytes). Por lo tanto, es importante asegurarse de que se dispone de los recursos necesarios (véase la sección final más abajo) para ejecutar el alineamiento en un tiempo razonable y almacenar los resultados. Del mismo modo, debido al tamaño y al formato binario de los archivos BAM, evite abrirlos con editores de texto; en su lugar, utilice comandos Unix o herramientas dedicadas como SAMtools.

De los alineamientos

Diría que no hay un paso común claro después del alineamiento, ya que en este punto es donde cada tipo de datos HTS y aplicación puede diferir.

Un análisis posterior común para los datos de DNA-seq es la llamada de variantes, es decir, la identificación de posiciones en el genoma que varían en relación con la referencia del genoma y entre individuos. Un marco de análisis popular para esta aplicación es GATK para polimorfismos de un solo nucleótido (SNP) o pequeñas inserciones/deleciones (indels) (Figura 2). Las variantes que comprenden trozos más grandes de ADN (también denominadas variantes estructurales) requieren métodos de llamada específicos (véase este artículo para una comparación exhaustiva). Al igual que con los alineadores, aconsejo seleccionar la herramienta adecuada teniendo en cuenta factores clave como el tipo de variantes (SNP, indels o variantes estructurales), la aceptación por parte de la comunidad, la calidad de la documentación y el número de usuarios.

Probablemente la aplicación más frecuente de RNA-seq es la cuantificación de la expresión génica. Históricamente, las lecturas debían alinearse con la secuencia de referencia y, a continuación, el número de lecturas alineadas con un determinado gen o transcripción se utilizaba como indicador para cuantificar sus niveles de expresión. Este enfoque de alineación+cuantificación se realiza mediante herramientas como Cufflinks, RSEM o featureCounts. Sin embargo, este enfoque ha sido superado cada vez más por nuevos métodos implementados en software como Kallisto y Salmon. Conceptualmente, con estas herramientas no es necesario alinear la secuencia completa de una lectura con la secuencia de referencia. En su lugar, sólo necesitamos alinear suficientes nucleótidos para estar seguros de que una lectura se originó a partir de un determinado transcrito. En pocas palabras, el enfoque de alineación+cuantificación se reduce a un solo paso. Este enfoque se conoce como pseudo-mapeo y aumenta en gran medida la velocidad de la cuantificación de la expresión génica. Por otro lado, hay que tener en cuenta que el pseudo-mapeo no será adecuado para aplicaciones en las que se necesita la alineación completa (por ejemplo, la llamada de variantes a partir de datos de ARN-seq).

Otro ejemplo de las diferencias en los pasos de análisis posteriores y las herramientas necesarias en las aplicaciones basadas en la secuenciación es ChIP-seq. Las lecturas generadas con esta técnica se utilizarán para la llamada de picos, que consiste en detectar las regiones del genoma con un exceso significativo de lecturas que indica dónde se une la proteína objetivo. Existen varios «peak callers» y esta publicación los estudia. Como último ejemplo mencionaré los datos Hi-C, en los que los alineamientos se utilizan como entrada para herramientas que determinan las matrices de interacción y, a partir de ellas, las características 3D del genoma. Comentar todos los ensayos basados en la secuenciación va más allá del alcance de este artículo (para una lista relativamente completa ver este artículo).

Antes de empezar…

La parte restante de este artículo toca aspectos que pueden no ser considerados estrictamente como pasos en el análisis de datos HTS y que son ampliamente ignorados. Por el contrario, sostengo que es capital que piense en las preguntas planteadas en la Tabla 1 antes de empezar a analizar los datos de HTS (o cualquier tipo de datos, de hecho), y he escrito sobre estos temas aquí y aquí.

Tabla 1

Piensa en ello	Acción propuesta
¿Tienes toda la información de tu muestra necesaria para el análisis?	Recoger sistemáticamente los metadatos de los experimentos
¿Podrá identificar inequívocamente su muestra?	Establecer un sistema para asignar a cada muestra un identificador único
¿Dónde estarán los datos y los resultados?	Organización estructurada y jerárquica de los datos
¿Será capaz de procesar múltiples muestras sin problemas?	Escalabilidad, paralelización, configuración automática y modularidad del código
¿Podrá usted o cualquier otra persona reproducir los resultados?	¡Documente su código y sus procedimientos!

Como se ha mencionado anteriormente, los datos brutos de HTS y algunos de los archivos generados durante su análisis son del orden de gigabytes, por lo que no es excepcional que un proyecto que incluya decenas de muestras requiera terabytes de almacenamiento. Además, algunos pasos del análisis de los datos HTS son intensivos desde el punto de vista informático (por ejemplo, la alineación). Sin embargo, la infraestructura de almacenamiento y computación necesaria para analizar los datos de HTS es una consideración importante y a menudo se pasa por alto o no se discute. Como ejemplo, como parte de un análisis reciente, revisamos decenas de artículos publicados que realizaban análisis de asociación de todo el fenotipo (PheWAS). Los PheWAS modernos analizan entre 100 y 1.000 variantes genéticas y fenotipos, lo que supone una importante capacidad de almacenamiento de datos y de computación. Y, sin embargo, prácticamente ninguno de los artículos que revisamos comentaba la infraestructura necesaria para el análisis PheWAS. No es de extrañar que mi recomendación sea que planifique por adelantado los requisitos de almacenamiento y computación a los que se enfrentará y los comparta con la comunidad.

¿Necesita ayuda para analizar los datos de secuenciación de ADN? Ponte en contacto con bioinformáticos autónomos y expertos en genómica en Kolabtree.