Assingment 1. NGS

Assignment 1. Práctica NGS

Introducción

El siglo XX se caracterizó por una multitud de hallazgos y avances en el campo de la genética, tales como la definición de la estructura de doble hélice del DNA, el splicing del RNA, los codones formados por tripletes de nucleótidos, la técnica de PCR y los primeros intentos de secuenciación. A finales de dicho siglo, se desarrolló la técnica de secuenciación Sanger, más tarde empleada en el Proyecto Genoma Humano (de sus siglas en inglés, HGP) para secuenciar, por primera vez, la totalidad del genoma humano. Esta fue la primera generación de técnicas de secuenciación.

En estas últimas dos décadas, se han desarrollado numerosas técnicas de secuenciación formando las técnicas de segunda (Next Generation Sequencing, NGS) y tercera generación. Estas nuevas técnicas han permitido ir disminuyendo notablemente el coste monetario, el tiempo y el capital humano invertido a la hora de secuenciar. A día de hoy, la compañía Illumina podría decirse que es el principal referente en desarrollo de plataformas de secuenciación, pues poseen una amplia gama de dispositivos de secuenciación para adaptarse a las distintas necesidades de sus clientes. Aunque Illumina actualmente lidera el sector, la competencia también crece, desarrollando tecnologías y dispositivos más económicos.

En este proyecto, vamos a emplear un archivo de lecturas (reads) genómicas de E. coli obtenidas de un dispositivo Illumina. Este archivo se encuentra en formato FASTQ, un formato de archivo de texto codificado en ASCII, muy utilizado para guardar secuencias genéticas y un puntaje de calidad asociado. Se caracteriza por estar formado por cuatro filas distintas: la primera fila empieza por ‘@’ y contiene un identificador de la secuencia, la segunda fila incluye la secuencia como tal, la tercera es un ‘+’ que actúa como separador y en la cuarta está codificada la calidad de la secuencia.

El objetivo del trabajo será realizar un análisis NGS a partir del archivo fastq con las reads de E. coli, el análisis incluirá etapas de procesamiento, reportes de calidad, alineamiento y visualización de los datos. Todos estos pasos los discutiremos en profundidad en el siguiente apartado.

Materiales y Métodos

El flujo de trabajo seguido en esta práctica puede resumirse en los siguientes pasos:

• Obtención del archivo fastq (raw data)
• Control de calidad del raw data (FastQC)
• De-multiplexing de las muestras utilizando DNA barcodes (Cutadapt)
• Control de calidad de las reads ya procesadas (MultiQC)
• Alineamiento frente a genoma de referencia (Bowtie2)
• Transformar archivos sam en bam, extraer estadísticas relevantes (Samtools)
• Recoger todas las estadísticas obtenidas (MultiQC)
• Visualizar las reads alineadas (IGV)

A continuación, se explican los programas y el sistema operativo empleados durante el proyecto:

UNIX

Utilizamos el terminal del sistema operativo de UNIX para el desarrollo de la totalidad de la práctica para: creación de directorios, descarga de archivos, ejecución de programas, etc.

FastQC

Esta herramienta lee archivos de secuencia y produce un informe de control de calidad con distintos apartados que nos permite discernir si las secuencias analizadas presentan o no problemas (Andrews, S., 2010).

Cutadapt

Nos permite encontrar y eliminar secuencias adaptadoras, primers, colas poli-A y otras secuencias no deseadas en nuestras secuencias de alto rendimiento (Martin, M., 2011).

MultiQC

Se trata de una herramienta que permite crear un reporte único recogiendo los resultados de varias herramientas bioinformáticas a la vez (Ewels et al., 2016).

Bowtie2

Este programa permite realizar alineamientos de reads de secuenciación con secuencias largas de referencia de forma rápida y altamente eficiente (Langmead and Salzberg, 2012)

Samtools

Es un conjunto de programas (Samtools, BCFtools y HTSlib) que permiten realizar distintas acciones sobre datos de alto rendimiento. Nosotros empleamos únicamente los usos del paquete de Samtools, que permite leer, escribir, indexar, editar y visualizar archivos con formato SAM, BAM y CRAM (Danecek et al., 2021).

Integrative Genome Viewer (IGV)

Se trata de una herramienta interactiva de fácil uso para la exploración visual de datos genómicos. Es altamente flexible en los formatos admitidos (Robinson et al., 2011).

Blast

Blast proviene del acrónimo “Basic Local Alignment Search Tool”, se trata de un algoritmo de búsqueda utilizado para comparar secuencias de ADN, ARN o proteínas con secuencias biológicas almacenadas en bases de datos públicas (Altschul et al., 1990).

Resultados y Discusión

Cuestión 1

Enunciado: En la primera práctica se le entregó un archivo fastq que contenía lecturas de un experimento simulado de secuenciación del genoma completo. Cuando las lecturas se dividieron utilizando los cuatro códigos de barras (Negative, Positive, Long y BQ), el archivo Negative.fq dio unos resultados de mapeo muy pobres (véase la figura siguiente).

Siguiendo las buenas prácticas de NGS, reasigne Negative.fq para que las estadísticas de mapeo mejoren. Debería poder obtener resultados de mapeo similares a los obtenidos con Positive.fq.

Partimos del archivo original que posee todas las reads de un experimento de secuenciación del genoma completo. Posteriormente, realizamos un “trimming” de los DNA barcodes de las secuencias, pues conocemos la secuencia de los barcodes. De este modo, dividimos las reads, según el DNA barcode, en cuatro: Negative, Positive, Long y BQ.

El archivo Negative.fq posee unos resultados de mapeo pobres y nuestro objetivo es obtener resultados de mapeo similares a los obtenidos con Positive.fq. Esta primera suposición puede comprobarse con un informe FastQC, donde podemos observar en el gráfico de Mean Quality Scores que trimmed_negative (resaltado en morado), posee una calidad baja en las posiciones iniciales y finales.

Entonces, al observar la gráfica de Per Base N Content, nos podemos dar cuenta del problema. El inicio y el final de estas reads tienen N, nucleótidos indeterminados (donde el mapping no ha sabido identificar qué nucleótido corresponde). Por tanto, el siguiente paso será realizar un trimming de estos nucleótidos y hacer un posterior alineamiento de las reads ya procesadas, para ver si el cambio realizado ha sido suficiente.

Vamos a probar distintas formas de resolver el problema, para ello, probaremos distintas opciones de la herramienta Cutadapt. Las opciones trim5 y trim3 nos permiten trimmear selectivamente, los extremos 5’ y 3’ respectivamente, la cantidad de nucleótidos que le especifiquemos. Viendo las secuencias podemos ver que en el extremo 5’ hay 5 N y en el extremo 3’ hay 4 N.

Por último, vamos a probar cambiar los métodos de alineamiento que posee Cutadapt, end-to-end y local.

Tras diversas pruebas, utilizando las distintas opciones por separado, las mejores formas de pulir las reads fue utilizar ambas opciones de trim5 y trim3, junto con un método end-to-end o local. Aquí obtuvimos dos formas:

• Negative stats2 (trim5 + trim3, end-to-end)

• Negative stats3 (trim5 + trim3, local)

Como podemos observar en la imagen de las puntuaciones de alineamiento de Bowtie 2, Negative_bowtie_stats2 es el que mejor se acerca al resultado obtenido en Positive, siendo prácticamente el mismo porcentaje de secuencias sin mapear y de secuencias mapeadas una única vez o multimapeadas. Mientras que Negative_bowtie_stats3, mejora el porcentaje de secuencias sin mapear, pero el multimapeo aumenta, lo que realmente es contraproducente.

No siempre va a ser mejor tener menos reads sin alinear, si para ello sacrificamos secuencias mapeadas una única vez y empieza a aumentar el multimapeo, esto nos refleja que hay ciertas secuencias que comienzar a alinear en distintos puntos del genoma de referencia, cuando se supone que cada read debería ser única.

Del mismo modo, podemos ver que, aparentemente Negative stats3 podría ser mejor que la primera opción, porque tiene un error rate menor (0.10% frente a 0.11%) y ambas tienen un 100% de mapeo. Pero esto no es así, porque consideramos mejor tener menor multimapping.

Cuestión 2

Enunciado: El archivo BQ.fq también dio unas estadísticas de mapeo muy pobres durante la primera práctica. Es posible lograr tasas de alineación mucho más altas, pero en este caso también queremos minimizar la tasa de error de alineación (es decir, los nucleótidos no sólo deben alinearse, sino alinearse correctamente).

Con este objetivo en mente, utilice cutadapt para reprocesar el archivo BQ.fq y bowtie2 para mapear las lecturas con el genoma de referencia.

Para empezar, vamos a analizar una de las gráficas obtenidas del FastQC, en rojo tenemos resaltado los trimmed_BQ. Esto nos da una idea de los principales problemas que presentan estas reads. El inicio y el final de estas tiene un puntaje de calidad extremadamente bajo, dejándonos ver que en esas bases de cada read hay algún problema subyacente. Las bases intermedias presentan una calidad mediocre, de 21 puntos, las bases no tienen una buena calidad, en general.

Como podemos observar, el primer problema que poseen las reads BQ es que tienen N’s tanto al principio como al final de las reads.

Esta vez vamos a eliminarlas utilizando la función “-- trim-n” de Cutadapt, a la vez que trimeamos los adaptadores de la secuencia.

El segundo problema que presentan estas reads es la calidad media de los nucleótidos, de hecho, si abrimos el archivo .fq podemos observar el Phred score dado a cada nucleótido. La gran mayoría de ellos posee un Phred score expresado en ASCII de 7 e, intercaladamente, podemos encontrar valores inferiores de 3, 1, 0 y /. Estos valores son interpolables a la escala 0-41 que es la que podemos observar en el reporte de FastQC, para ello, buscamos en la tabla de Phred Quality Score por L - Illumina 1.8+ Phred+33 (Illumina, 2011).

A pesar de conocer cuál es el posible problema de las reads de BQ, no logré resolverlo aun intentándolo de muy diversas maneras. Una buena aproximación fue utilizar la función “-q <int>” de Cutadapt que permite filtrar las reads por un quality score que indiques en <int>, aun así esta función no me funcionó.

Por tanto, continué probando con Bowtie2 y la mejor opción fue la siguiente:

• La función “-N 1” establece el número máximo de mismatches en el seed alignment (puede ser 0 ó 1)

De esta forma, como se puede observar en la siguiente imagen, se obtiene un porcentaje de mapeo del 97,8% con un error rate del 5,13%. Esta es la opción que mejoraba notablemente el porcentaje de mapeo sin sacrificar tanto el error rate. Como se trata de un experimento de secuenciado genómico, considero que es importante tener un porcentaje de mapeo, por lo menos, por encima del 90%. El error rate pude llegar a reducirlo a un 3-4%, pero con un porcentaje de mapeo del 89%. Un mapeo cercano al 100% considero que compensa un aumento de una unidad porcentual de error rate.

Cuestión 3

En esta tarea tiene que familiarizarse con el uso de samtools y los SAM flags. Su objetivo es dividir el archivo alineado BQ (SAM) en archivos SAM que contengan subconjuntos del conjunto completo de lecturas alineadas.

Divida el archivo SAM original en:

Pregunta a)

Un archivo que contenga las lecturas de mapeo múltiple y otro archivo que contenga sólo las lecturas de mapeo único únicamente.

Para obtener las secuencias de lectura única (uniquely mapped), podemos excluir usando “-F” con la flag 260 (256 + 4) tanto las secuencias no alineadas (flag 4) y las secuencias provenientes de alineamientos no primarios (flag 256) (Broad Institute of MIT and Harvard)

Para obtener las secuencias que mapean multiples veces, podemos utilizar como filtro “-f” el flag 256 correspondiente a los alineamientos no primarios y excluir con “-F” el flag 4 que corresponde a las secuencias no alineadas (el contrario con -F 4 serían todas las secuencias alineadas) (Broad Institute of MIT and Harvard).

Pregunta b)

Un archivo que contenga lecturas no mapeadas y otro archivo que contenga únicamente lecturas mapeadas.

Para este caso, he empleado el flag 4 del formato SAM/BAM, que hace referencia a lecturas “unaligned” o no alineadas. Con “-f” filtramos aquellas lecturas por el flag que le indiquemos como filtro. Y con la opción “-h” imprimimos el encabezado (header).

Por otro lado, para las lecturas que sí han mapeado, utilizamos la opción “-F” que filtra las lecturas de un archivo SAM/BAM que no contienen la flag indicada como filtro. Entonces, usamos la flag 4 que era para las no alineadas (hacemos el contrario de las no alineadas –> las alineadas).

En esta última imagen, recogemos los resultados obtenidos en los apartados a y b. BQ_stats3 corresponde al archivo usado de partida, por tanto, incluye la totalidad de las secuencias, tanto alineadas como no alineadas.

Empezando por las secuencias alineadas, vemos que corresponde correctamente con las BQ_stats3 alineadas, pues existen aproximadamente 1050k en ambos casos. Por otro lado, las no alineadas son unos 25k y corresponden con la fracción roja de BQ_stats3.

Las mapeadas únicamente son algo menos de 1050k y se puede corroborar con Bowtie 2 (no enseñado) donde se ve que las mapeadas únicamente corresponden con un 92% aproximadamente, frente al 97% total del mapeado, por tanto, las multimapeadas son el 5% de diferencia. Ese sería el resultado esperado para las multimapeadas, pero que finalmente, no se pudo conseguir.

Pregunta c)

¿Cómo puede obtener la salida de (b) utilizando Bowtie2 en su lugar?

Para esta cuestión, utilicé la opción de antes de cutadapt para trimear el DNA barcode y adicionalmente eliminar las N’s de BQ.

Posteriormente, utilicé las mismas opciones de Bowtie 2 para conseguir un 97,8% de mapeo, pero para conseguir la salida de la cuestión (b), empleé la opción “--un” para generar un output en formato .fastq que llamé unmapped3c. Por último, usé la opción “--al” para crear un archivo con formato .sam con las secuencias alineadas, bajo el nombre de aligned3c. Los formatos escogidos en los outputs son los que me permitían las opciones de Bowtie 2.

Pregunta d)

Inspeccione las unmapped reads y utilice BLAST para identificar su origen.

Partiendo del archivo generado en la cuestión 3(b), el primer paso para realizar un BLAST es pasar nuestro archivo en formato .bam a formato .fasta. Posteriormente, usando el comando “sed” de UNIX, pude eliminar los headers (cabeceros), para quedarme únicamente con las secuencias.

Para realizar el BLASTn (Nucleotide BLAST), escogí las diez primeras reads a modo de ejemplo.

El resultado obtenido nos indica que las secuencias que han obtenido un Max score mayor con un E value muy significativo corresponden con secuencias provenientes de Homo sapiens.

Esto explica la falta de alineamiento de estas secuencias, pues estamos trabajando con E. coli y, parece ser, que estas secuencias que no alinean provienen de Homo sapiens, pudiéndonos indicar que nuestros reads se encuentran contaminados con secuencias de humano, posiblemente de la extracción y purificación de las reads originales.

Bibliografía

Altschul, S.F., Gish, W., Miller, W., Myers, E.W. & Lipman, D.J. (1990) “Basic local alignment search tool.” J. Mol. Biol. 215:403-410. PubMed

Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data [Online]. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

Broad Institute of MIT and Harvard (no date) Picard, Explain SAM Flags. Available at: https://broadinstitute.github.io/picard/explain-flags.html (Accessed: March 9, 2023).

Ewels, P., Magnusson, M., Lundin, S., & Käller, M. (2016). MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics, 32(19), 3047-3048. https://doi.org/10.1093/bioinformatics/btw354

Illumina (2011). Quality Scores for Next-Generation Sequencing. Available at: https://www.illumina.com/documents/products/technotes/technote_Q-Scores.pdf. (Accessed: March 7, 2023)

James T. Robinson, Helga Thorvaldsdóttir, Wendy Winckler, Mitchell Guttman, Eric S. Lander, Gad Getz, Jill P. Mesirov. Integrative Genomics Viewer. Nature Biotechnology 29, 24–26 (2011).

Langmead B, Salzberg S. Fast gapped-read alignment with Bowtie 2. Nature Methods. 2012, 9:357-359

Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal, 17(1), pp. 10-12. doi:https://doi.org/10.14806/ej.17.1.200

Philip Ewels, Måns Magnusson, Sverker Lundin, Max Käller, MultiQC: summarize analysis results for multiple tools and samples in a single report, Bioinformatics, Volume 32, Issue 19, October 2016, Pages 3047–3048, https://doi.org/10.1093/bioinformatics/btw354

Petr Danecek, James K Bonfield, Jennifer Liddle, John Marshall, Valeriu Ohan, Martin O Pollard, Andrew Whitwham, Thomas Keane, Shane A McCarthy, Robert M Davies, Heng Li. Twelve years of SAMtools and BCFtools. GigaScience, Volume 10, Issue 2, February 2021, giab008, https://doi.org/10.1093/gigascience/giab008