ENSAMBLADO DE UN GENOMA BACTERIANO PEQUEÑO (PARTE I)

Prefacio

La presente clase consistirá en explicar los pasos que conlleva realizar un ensamblado de novo de un genoma mediante la utilización de un ensamblador a partir de las lecturas (reads) obtenidas por secuenciación masiva. En este caso utilizaremos el ensamblador VELVET ya que ha sido probado en genomas bacterianos pequeños con muy buen resultado, y CANU específicamente para tecnologías de secuenciación de tercera generación (reads largos).

Objetivo

El objetivo de esta práctica es introducir un workflow clásico de trabajo para lograr el ensamblado de un genoma bacteriano de tamaño pequeño a partir de reads pareados cortos y reads largos.

Luego de esta clase práctica usted debería poder:

  • Computar, investigar y evaluar la calidad de los datos de secuencia provenientes de un experimento de secuenciación masiva.
  • Describir los pasos involucrados en el pre-procesamiento/limpieza de los datos de secuenciación.
  • Distinguir entre una buena y mala corrida/experimento de secuenciación.
  • Computar, interpretar y evaluar el ensamblado de un genoma pequeño completo.

Introducción

El desarrollo de métodos para elucidar la secuencia del ADN dentro de los genomas de los organismos ha revolucionado la investigación biomédica y ha llevado a la actual revolución genómica. La secuencia genómica de decenas de miles de bacterias, virus y humanos y cientos a miles de genomas de otros organismos han sido reconstruidos. Debido a que las tecnologías de secuenciación de ácidos nucleicos de alto rendimiento son cada vez más baratas, rápidas y sensibles, hoy pueden considerarse como parte vital de la caja de herramientas de un investigador. A pesar de los tremendos avances en las tecnologías de secuenciación del ADN, los instrumentos solo pueden generar lecturas de ADN cortas (en comparación al tamaño genómico general) desde 100 pb hasta 10-20 Kb en tecnologías de cuarta generación (o tercera, según autor) como PacBio o Nanopore. Teniendo en cuenta que el genoma humano tiene 3000 millones de bases y algunas bacterias llegan a tener 10 millones de bases (10 Mb) reconstruir un genoma a partir de pequeños fragmentos es un desafío computacional muy grande. El secuenciamiento de genomas bacterianos y de pequeños eucariotas puede realizarse in-house en pequeños laboratorios debido a la existencia de secuenciadores (como: MiSeq (Illumina); Ion GeneStudio S5 System (Ion Torren); Roche 454 (FLX Junior); MinION (ONT)) de pequeño tamaño y costo. En los últimos años han aparecido nuevas tecnologías que producen lecturas más largas, pero a un poseen un mayor costo (PacBio).

Diseño Experimental

Como cualquier buen proyecto, un buen ensamblado de novo empieza con un correcto diseñoexperimental. Aspectos biológicos, técnicos y experimentales deben de ser considerados.

Aspectos biológicos

  • ¿Qué tan grande es el genoma que deseo secuenciar? El tamaño del genoma es importante para elegir el ensamblador (programa que ensamblara las lecturas (reads) en secuencias continuas de mayor tamaño), existen ensambladores que funcionan mejor con genomas pequeños mientras que otros están optimizados para genomas más grandes.
  • ¿Cuál es el contenido de repetidos del genoma que deseo ensamblar? Cuanto mayor sea el contenido de repetidos de un genoma más difícil será el correcto ensamblaje del genoma y seguramente requiera la utilización de read largos o de mate-pairs para obtener un buen ensamblado que abarquen grandes regiones genómicas.
  • ¿Qué tan rico (o pobre) es en contenido AT? Genomas con un gran desbalance AT/GC (en cualquier sentido) tienden a tener bajo contenido de información y complejizan la secuenciación y el ensamblado.
  • ¿El genoma en cuestión, es haploide, diploide o poliploide? Actualmente los ensambladores genómicos tienen mejor performance con muestras haploides, y algunos ensambladores proveen genomas haploides con sitios heterocigotos. Los genomas poliploides (ejemplo: plantas) suelen ser más problemáticos.

Aspectos experimentales

¿Es posible extraer mucho ADN? Si solo se dispone de poco material es posible que se deba amplificar la muestra, pero hay que ser consiente que este procedimiento traerá sesgos en los datos. - ¿El ADN proviene de una única célula, de una población clonal, o de una colección heterogénea de células? La diversidad presente en la muestra podrá traer menor o mayor cantidad de ruido con la cual los diferentes ensambladores tratan de diferente manera.

Aspectos técnicos

  • ¿Cuánto cuesta la secuenciación?
  • ¿Cuál es la profundidad de secuenciación adecuada? Cuanto más ruido en la muestra mayor será la profundidad necesaria.
  • ¿Cuán largo son los reads? ¿ Single-end reads o Paired-end reads? Cuantos más largos sean las lecturas (reads) mayor facilidad habrá para desambiguar secuencias repetidas.
  • ¿Conviene utilizar más de un tipo de read, que estrategia debo utilizar? Por ejemplo, reads cortos en gran profundidad (de menor costo) junto a reads largos en menor profundidad (mayor costo).

Diseño Computacional

Al proceso de descifrar la secuencia genómica a partir de pequeños fragmentos de ADN, en conjunto con alguna información adicional disponible, se le denomina ensamblado de genomas.Las estrategias para el ensamblado de genomas se pueden dividir en dos categorías: ensamblaje por comparación, en el que se utiliza un genoma como referencia (ensamblaje con referencia); y ensamblaje de novo, en el cual se utiliza solo la información obtenida de la secuenciación para reconstruir el genoma de interés, sin conocimiento a priori de la organización del mismo. Sin embargo, en esta última estrategia algunas informaciones previas son útiles, como el tamaño esperado del genoma, el contenido en GC y el contenido de regiones repetitivas, ya que ayudan a elegir la mejor estrategia a seguir. Estos datos pueden ser inferidos a partir de secuencias de organismos relacionados.

De una manera general podemos decir que el flujo de trabajo para realizar un ensamblaje de novo involucran:

  1. control de calidad y corrección de errores en los datos
  2. ensamblaje
  3. evaluación de la calidad del ensamblaje
  4. terminación del genoma.

Actividades

Pre-procesamiento de los datos

El propósito de esta sección es demostrar cómo se puede analizar los datos de secuenciación crudos (raw data) de manera de evaluar su calidad. El control de calidad de las lecturas (reads) sirve como un chequeo rápido para identificar y excluir datos con problemas serios de calidad, lo cual permite ahorrar gran cantidad de tiempo en los análisis posteriores. Las herramientas empleadas chequean la calidad de las bases (probabilidad de que la base asignada sea la correcta), la distribución de los nucleótidos, la distribución del contenido de GC, secuencias repetidas entre otros parámetros.

El primer paso es obtener los reads. Ya sean datos propios (un experimento in-house) o recuperados de base de datos públicas. Para esta demostración utilizaremos datos de Illumina MiSeq (pair-ends) de E. coli K-12 cepa MG1655 que se encuentran disponibles en la base de datos Illumina. Es importante que bajemos ambos archivos de lectura (correspondientes a los reads forward y reverse). En el caso que quisiéramos bajar los datos directamente del NCBI la herramienta que podríamos utilizar (parte del paquete SRA-TOOLKIT) sería:

fastq-dump --split-3 ERR*
wget https://www.dropbox.com/s/4686pqo0oiwjbju/MiSeq_Ecoli_MG1655_50x_R1.fastq
wget https://www.dropbox.com/s/d10h4qulfzr4c0k/MiSeq_Ecoli_MG1655_50x_R2.fastq

Ustedes NO tienen que descargarlos, ya que cuentan con los archivos FASTQ en la carpeta DATA/Ensamblado_DATA/E_coli_K12_illumina/. Lo primero que tienen que hacer es copiarse estos archivos a sus carpetas personales.

Recuerden SIEMPRE trabajar en la carpeta que les corresponde. En caso de no recordar dónde se encuentran en el servidor, el comando que deben utilizar es pwd que les devuelve su ubicación.

Para copiar los archivos:

Usando ruta absoluta:

cp /media/Disco10/curso_genomica/estudiantes/2023/DATA/Ensamblado_DATA/E_coli_K12_illumina/* ./

Usando ruta relativa:

cp ../DATA/Ensamblado_DATA/E_coli_K12_illumina/* ./

El * representa todos los archivos que existan en la carpeta Illumina_reads

Con respecto a las lecturas originales que usted descargo en su carpeta:

  1. ¿Cuántas lecturas existen? (Recuerde comando grep)
  2. ¿Cuál es la segunda secuencia en R1? ¿Y en R2?
  3. Obtener los primeros 10 IDs de R1. Reportar el comando utilizado.

Los archivos provenientes del secuenciador están en el formato FASTQ. Uno de los software más utilizados para analizar la calidad de los datos provenientes de tecnologías Illumina es FASTQC. Una de las ventajas de este software es que posee interfaz gráfica, lo que lo hace de muy fácil uso. FASTQC utiliza diversos reportes para determinar la calidad de las secuencias. Dentro de los más destacados están el apartado de “calidad por base” que aporta información de la calidad de cada uno de los nucleótidos que forman la secuencia, es decir, la calidad con la que se puede precisar que cada base esté, realmente, en cada una de esas posiciones. En la mayoría de los casos, las bases nucleotídicas van perdiendo calidad a medida que nos acercamos al extremo 3`debido a la propia naturaleza de los cebadores y química de la síntesis, los cuales se van despegando y la secuenciación (base calling) va perdiendo calidad. Por otro lado, también se puede analizar el “contenido de base en la secuencias” lo que me permite ver la proporción de las bases en cada posición de la lectura. En una muestra fragmentada al azar uno no esperaría que haya diferencias en el contenido de nucleótidos en cada posición. La cantidad relativa de cada base debería reflejar la cantidad total de estas bases en el genoma, pero en cualquier caso no debería ser muy desequilibrada unos de otros. La observación de fuertes sesgos suele indicar la presencia de una secuencia sobrerrepresentada que está contaminando la librería.

Para analizar la calidad de nuestros datos ejecutaremos el siguiente comando:

fastqc MiSeq_Ecoli_MG1655_50x_R1.fastq MiSeq_Ecoli_MG1655_50x_R2.fastq

Este comando crea dos archivos .html y dos archivos .zip. Descomprima las carpetas, entre a cada una de ellas y visualice el reporte (recuerde el comando unzip).

unzip MiSeq_Ecoli_MG1655_50x_R1_fastqc.zip
unzip MiSeq_Ecoli_MG1655_50x_R2_fastqc.zip

Para visualizar los archivos .html:

firefox MiSeq_Ecoli_MG1655_50x_R1_fastqc.html

También pueden visualizar 1 a 1 las imágenes dentro de Images/:

display per_base_quality.png

Una vez que se ha analizado el reporte de FASTQC y se posee conocimiento acerca de la calidad de los datos de secuenciación, es importante utilizar esta información para recortar y filtrar las lecturas para mejorar su calidad global antes del ensamblado. Hay una serie de herramientas disponibles por línea de comandos que pueden realizar este paso, pero tendrá que ser una que pueda manejar las lecturas pareadas (ya que el ejemplo de este práctico es con lecturas pareadas).

  1. ¿Es estrictamente necesario recortas las secuencias? Justifique
  2. ¿Si tuviera que recortar bases de un extremo, de cuál sería?
  3. ¿Qué tipo de adaptadores esperaría encontrar?

La herramienta elegida para recortar las lecturas de acuerdo a su calidad es TRIMMOMATIC. El filtrado ocurre en el orden en que los pasos están especificados en la línea de comandos. En la mayoría de los casos, se recomienda que el recorte de adaptadores, si es necesario, se realice lo antes posible. Por defecto reconoce los adaptadores de Illumina. - Illumina: AGATCGGAAGAGC

El comando para recortar las secuencias utilizando TRIMMOMATIC es:

trimmomatic PE -phred33 -threads 2 -trimlog logfile MiSeq_Ecoli_MG1655_50x_R1.fastq MiSeq_Ecoli_MG1655_50x_R2.fastq Ecoli_MG1655_50x_R1_P Ecoli_MG1655_50x_R1_U Ecoli_MG1655_50x_R2_P Ecoli_MG1655_50x_R2_U ILLUMINACLIP:/home/estudiantes/miniconda3/pkgs/trimmomatic-0.39-hdfd78af_2/share/trimmomatic-0.39-2/adapters/TruSeq3-PE.fa:2:30:10:2:keepBothReads LEADING:25 TRAILING:25 MINLEN:50

Tambien les dejamos otra opción de filtrado utilizando TRIM GALORE (NO es para correrlo en clase)

trim_galore --paired --three_prime_clip_R1 10 --three_prime_clip_R2 10 MiSeq_Ecoli_MG1655_50x_R1.fastq MiSeq_Ecoli_MG1655_50x_R2.fastq

Ejecute el comando FASTQC para visualizar las lecturas trimeadas:

fastqc Ecoli_MG1655_50x_R1_P
fastqc Ecoli_MG1655_50x_R2_P

Visualice nuevamente los archivos .html

¿Cree que el recorte mejoró las chances de lograr un buen ensamblado? Justifique

Una vez filtradas las secuencias estamos listos para el siguiente paso que es el ensamblado de ellas.

Ensamblaje

El propósito de esta sección es delinear el proceso de ensamblaje de los read recortados de alta calidad en contigs/scaffolds(draft contigs). La mayoría de los ensambladores necesitan una serie de parámetros que deben ser establecidos antes de ejecutar el software. Estos parámetros pueden tener y tienen un gran efecto sobre el resultado de cualquier ensamblado. Los ensamblados producidos pueden tener más o menos gaps. Por lo tanto, el conocimiento de los parámetros y sus efectos es esencial para obtener buenos ensamblados. En la mayoría de los casos, un conjunto óptimo de parámetros de los datos se puede encontrar mediante un método iterativo. Todos los métodos de ensamblaje se basan en la simple suposición de que fragmentos de ADN altamente similares se originan de la misma posición dentro del genoma. De esta manera, la similitud entre secuencias de ADN se usa para conectar fragmentos individuales en secuencias contiguas más largas, denominadas contigs (secuencias consenso obtenidas a partir del ensamblaje de los reads). Los reads de corta longitud constituyen un problema cuando hay repeticiones en el genoma. Las regiones repetitivas son segmentos de ADN que aparecen más de una vez a lo largo del genoma. Cuando un read proviene de una región repetitiva, y es más corto que esta, no se sabe con certeza de cuál copia de la repetición se obtuvo. Es por ello que, durante el ensamblaje, se pueden crear falsas uniones en el genoma en las regiones de repeticiones. Además, debido a las regiones repetitivas, muchas veces resulta difícil ensamblar todos los fragmentos para que se logre, en un solo evento, reconstruir la secuencia del genoma completo, incluso en genomas microbianos pequeños.

Las estrategias empleadas por los programas ensambladores de secuencias pueden agruparse en tres paradigmas principales: Greedy, Overlap-Layout-Consensus y grafos de Bruijn. Debido a que el software que utilizaremos utiliza los grafos de Brujin, describiremos muy brevemente su algoritmia. El método de grafos de Brujin para ensamblado de secuencias tiene sus raíces en los trabajos teóricos a finales de la década del 80 realizados por Pevzner. Pevzner estudio el problema de reconstruir un genoma cuando solo un set de k-meros es conocido. Este trabajo teórico y subsiguientes desarrollos establecieron las bases del ensamblado genómico mediante grafos de Brujin. En este tipo de ensamblaje, cada read es fragmentado en secuencias solapantes de k letras (k-meros). Los distintos k-meros son agregados como vértices al grafo y los k-meros que se originan de una posición adyacente de un mismo reads son unidos con conectores. Entonces, el problema de ensamblaje puede formularse como encontrar un camino por el grafo en el que se visite cada nodo una única vez. En la práctica, los errores de secuencia y sesgos de muestreo oscurecen el grafo, por lo que normalmente no se busca un recorrido euleriano completo a través de todo el grafo. Incluso cuando un ciclo euleriano a través de todo el grafo es posible, es poco probable que refleje una secuencia precisa del genoma debido a la presencia de repeticiones. En la mayoría de casos, el ensamblador intenta construir contigs que consisten en las regiones inequívocas, no ramificadas del grafo.

Figura 1: Representación esquemática del grafo de Brujin. Cada nodo, representado aquí por un rectángulo, representa una serie de k-meros solapantes (k = 5), listados directamente abajo (en rojo el último nucleótido de cada k-mero). Los nucleótidos representados en rojo, figuran en letras grandes en el centro del rectángulo y son los que definen el nodo (se representa la secuencia y su complementaria). Las relaciones entre nodos se representan por conectores direccionados. El último k-mero de un nodo de origen solapa con el primero k-mero del nodo de destino.

Ensamblado de novo con lecturas cortas

El software elegido en esta demostración es VELVET. VELVET realiza ensamblados de ADN a partir de lecturas cortas mediante la manipulación de grafos de Brujin. Es capaz de formar largos contigs (N50 mayores a 150 Kb) a partir de lecturas pareadas. El ensamblaje mediante VELVET sucede en cuatro etapas: hashing (hashing the reads), construcción de grafos (Constructing the de Brujin graph), simplificación del grafo (simplification of the graph) y construcción del contig (read off the contigs). Estos pasos son ejecutados por las funciones velveth y velvetg. Velveth fragmenta los reads y construye un diccionario (hash-table) de todas las palabras de largo k (y su complementaria reversa), donde k es un parámetro definido por el usuario. Cada k-mero es almacenado una vez, pero se indica el número de veces que aparece. Velvetg entonces utiliza este diccionario e inicialmente agrega dé a uno en uno los distintos k-meros. Si el k-mero a agregar solapa (se considera solapamiento como cobertura de k-1) con un k-mero existente, se lo coloca adyacente a este y es conectado a través de un conector, por otra parte, si no posee solapamiento con ningún nodo del grafo, este k-mero constituirá un nuevo nodo aislado. Cada nodo guardará información de cuantas veces aparece en el diccionario. Dependiendo del camino que se siga diferentes secuencias pueden ser leídas. El siguiente paso, también realizado por velvetg, es la simplificación del grafo. El primer paso es la unión de cadenas, que implica que cuando existen dos nodos conectados sin ramificaciones (lo que indica un único camino posible) estos sean unidos. El siguiente paso es la eliminación de puntas (tip clipping) que consiste en remover cadenas de nodo que están ramificadas y poseen baja cobertura. Posteriormente se remueven las burbujas, que son caminos redundantes que empiezan y terminan con el mismo nodo y son producidos por errores de secuenciación, variantes biológicas u otros. Finalmente, velvetg remueve nodos mediante el cut-off establecido por el usuario (o establecido automáticamente) y realiza el read-off que se basa en encontrar el mejor camino euleriano para conectar los nodos y de esa manera establecer la secuencia de los contigs que luego formaran los scaffolds.

Existen dos estrategias principales para determinar los parámetros que brindaran el mejor ensamblado. El primero es ir ejecutando velveth y velvetg iterativamente cambiando los parámetros hasta lograr encontrar el mejor ensamblaje y el segundo se basa en utilizar un wrapper de VELVET, VELVET OPTIMISER que realiza esto por nosotros y elige el mejor ensamblaje de acuerdo al N50. Si bien la mejor opción es la segunda alternativa se mostrará como ejecutar velvetg y velveth para introducir los parámetros más importantes y luego se procederá a utilizar el wrapper para obtener el ensamblado final.

Uno de los parámetros más importantes en la ejecución de VELVET es hash length, que usualmente se lo denomina largo del k-mero, y corresponde al largo de las palabras que son almacenadas y comparadas para construir el grafo de Brujin. El largo del k-mero debe ser un número impar (para evitar palíndromos) más corto que la mayoría de los reads. Todos los reads más cortos que el hash serán ignorados en el ensamblaje. Encontrar el hash óptimo es esencialmente un compromiso entre sensibilidad y especificidad. Por un lado, cuanto más largo es el hash, menos probable es que las palabras se encuentren repetidas en el genoma. Una longitud de hash más grande conduce a un grafo de Brujin más simple con menos ramificaciones y errores. Sin embargo, cuanto mayor es la longitud de hash menor es el número de veces que se observa cada palabra en el conjunto de datos, despreciando errores y repeticiones. El número esperado de veces que un k-mero es observado en un dataset, o cobertura de k-mero Ck, puede ser derivado de la cobertura de reads (profundidad de secuenciación) C, el largo del read (L) y la longitud del hash (K) utilizando la siguiente formula:

\[ C_k = C^{\frac{L - K}{L}} \]

A medida que la longitud de hash se acerca a L, Ck puede caer rápidamente, creando así más vacíos de cobertura en el ensamblado. Por tanto, una longitud corta de hash da al ensamblador más sensibilidad para rellenar los gaps, debido a la mayor cobertura. La función que asocia a un N50 para cada longitud de hash tiene en general una forma de campana muy pronunciada. Por tanto, es en general bastante fácil de encontrar el valor óptimo con unas pocas iteraciones. Las diferentes plataformas de secuenciación producen fragmentos de distinta longitud y calidad, por lo tanto, diferentes k-meros permitirán un mejor ensamblado, dependiendo del dataset. Debe buscarse un equilibrio entre la sensibilidad ofrecida por un k-mero de menor tamaño y la especificidad ofrecida por uno de mayor tamaño.

  • Para ejecutar el comando de velveth, los vamos a dividir en 4 grupos, cada uno de ellos debe cambiar el largo del k-mero y elegir uno de los siguientes números: 19, 31, 51, 81.
velveth Assembly_velvet 31 -fastq -shortPaired -separate Ecoli_MG1655_50x_R1_P Ecoli_MG1655_50x_R2_P

¿Cuáles fueron los archivos creados? ¿Qué representan?

Una vez ejecutado velveth se ejecutará velvetg. Velvetg implica el uso de dos parámetros críticos. Expected coverage y coverage cut-off. La cobertura esperada es crucial, ya que permite a VELVET determinar que contigs corresponden a regiones únicas del genoma y cuales corresponden a repeticiones. Como se explicó anteriormente, la cobertura esperada se puede determinar utilizando una formula probabilística. Sin embargo, no toma en cuenta la tasa de error lo que disminuye la cobertura efectiva. Generalmente se opta por observar la distribución de la cobertura en un primer ensamblaje y luego se vuelve a ejecutar velveth/g con los parámetros más optimizados. El cut-off de cobertura es una manera simple pero eficaz de eliminar muchos errores básicos que pasaron los filtros anteriores. Es simplemente eliminar los contigs con una cobertura promedio menor al cut-off. El cut-off debe ser ajustado para asegurar que elimine más errores que contigs genuinos. Si el umbral es demasiado alto, entonces contigs genuinos son eliminados provocando la generación de malos ensamblados. Mientras que si el umbral esdemasiado bajo muchos errores no serán detectados. Si se decide establecer el modo auto el cut-off se ajusta simplemente a la mitad de la cobertura esperada.

  • Seleccione los parámetros de cobertura y tamaño de contig que quiera:
velvetg Assembly_velvet -clean yes -exp_cov 21 -cov_cutoff 2.81 -min_contig_lgth 200

Como resultado de esta acción se generan scaffolds y no contigs (debido a que se han utilizado reads pareados), sin embargo el programa genera un archivo denominado contigs.fa, puede renombrarlo si desea.

Como se ha mencionado anteriormente no sabemos si los parámetros elegidos son los óptimos para nuestros datos de secuencia a menos que iterativamente vayamos ajustando los parámetros en búsqueda de las mejores métricas. Para evitar el gasto humano y computacional de ejecutar iterativamente los diferentes valores utilizaremos el wrapper VelvetOptimiser.

VelvetOptimiser.pl -s 19 -e 81 -d AssemVelvetOptimiser -f '-shortPaired -fastq -separate Ecoli_MG1655_50x_R1_P Ecoli_MG1655_50x_R2_P' -t 8 --

Los parámetros s y e establecen el rango inferior y superior respectivamente de longitud de kmero. Para establecer el rango del k-mero recomendamos mirar los análisis iniciales de FASTQC y establecer el corte donde la calidad de las secuencias comienza a decaer.

Ensamblado de novo con lecturas largas

En esta oportunidad vamos a utilizar Canu, un sucesor del ensamblador Celera (fue desarrollado en Celera Genomics para el ensamblaje del genoma humano), diseñado específicamente para secuencias ruidosas de una sola molécula (PacBio RSII o Oxford Nanopore MinION). Canu utiliza nuevos algoritmos de superposición (OLC) y ensamblaje. El programa permite elegir el tipo de reads usando la opción -pacbio-raw para reads de PacBio, -nanopore-raw para reads de Nanopore. Debemos indicar el tamaño esperado del genoma, que en nuestro caso es 4,8 M. Puesto que Canu no usa información de calidad, también se podrían usar archivos FASTA en lugar de FASTQ.

  • Para la instalación de Canu (ustedes NO van a tener que instalarlo, el programa ya se encuentra instalado en el servidor)
conda install -c conda-forge -c bioconda -c defaults canu

Aplicando lo Aprendido (Tarea Domiciliaria):

La idea es que corran este comando en sus casas para el ensamblado de E. coli con lecturas largas. Canu ya se encuentra instalado en los servidores que van a utilizar. Recuerden siempre correrlo desde sus carpetas personales. - Para descargar los datos (ustedes NO tienen que descargarlos): (Ya están descargados en el servidor en DATA/)

wget http://gembox.cbcb.umd.edu/mhap/raw/ecoli_p6_25x.filtered.fastq
wget https://sra-pub-src1.s3.amazonaws.com/SRR10971019/m54316_180808_005743.fastq.1
wget https://figshare.com/articles/dataset/Ecoli_K12_MG1655_R10_3_HAC/11823087
  • Para copiarse los archivos de lecturas largas a sus carpetas personales:

!!!(parados dentro de su carpeta personal)

  • Usando ruta absoluta:
cp /media/Disco10/curso_genomica/estudiantes/2023/DATA/Ensamblado_DATA/E_coli_K12_PacBio/ecoli_HIFI.fastq ./
  • Usando ruta relativa:
cp ../DATA/Ensamblado_DATA/E_coli_K12_PacBio/ecoli_HIFI.fastq ./

Hagan los mismo para las lecturas de Oxford Nanopore.

  • Para correr el ensamblado de lecturas largas de PacBio/Nanopore con Canu:

    • Pacbio raw
    nohup canu -p ecoli-pacbio-raw -d Assembly_canu_PacbioRaw genomeSize=4.8m -pacbio-raw ecoli_p6_25x.filtered.fastq &
    • Pacbio HI-FI
    nohup canu -p ecoli-pacbio-hifi -d Assembly_canu_Pacbio_HIFI genomeSize=4.8m -pacbio-hifi ecoli_HIFI.fastq &
    • Nanopore
    nohup canu -p ecoli-oxford -d Assembly_canu_Nanopore genomeSize=4.8m -nanopore ecoli_oxford.fasta &

[Dar doble ENTER para que se libere la terminal] El proceso demora aproximadamente 20-30 minutos.

Con las opciones:

  • -d : Crea un directorio de salida con el nombre que especifiquemos
  • -p : Especifica el prefijo que van a tener nuestros archivos de salida

Una vez completado estos pasos, poseemos dos ensamblados de novo de un genoma pequeño, uno a partir de lecturas cortas y otro a partir de lecturas largas. Los pasos siguientes son para evaluar la calidad los ensamblados realizados.

Evaluación de la calidad del ensamblaje

Aunque utópico, el objetivo final del ensamblaje de novo es obtener un número de fragmentos igual al número total de cromosomas. En el caso de las bacterias es un cromosoma, aunque puede haber plásmidos que se obtienen como fragmentos separados. De acuerdo a esta idea, el investigador siempre debe buscar obtener el menor número de contigs posible. Existen algunos indicadores métricos que permiten evaluar la calidad del ensamblaje cuantitativamente. Se calcula generalmente el tamaño mínimo, medio y máximo de los contigs (scaffolds), así como el tamaño total del ensamblaje, el cual debe coincidir con el tamaño esperado del genoma. El principal valor estadístico es el valor N50, el cual corresponde con el menor de los mayores contigs (scaffolds) que cubren la mitad del genoma. Aunque el valor N50 constituye un indicador acerca de la contigüidad del genoma (i.e., cuán acertado fue el programa ensamblador en unir las secuencias contiguas en una única secuencia más larga), no es una señal absoluta de precisión y calidad del genoma ensamblado. Tampoco brinda una correcta estimación de si el tamaño del ensamblaje difiere o no del tamaño esperado. Sin embargo, hay que tener en cuenta que en general cuanto mayor sea el N50 mejor será el ensamblado. Es importante poder evaluar métricamente el ensamblado generado para poder determinar si fue o no un buen ensamblaje. Una manera de evaluar su resultado y obtener métricas es utilizar el script assemblathon_stats.pl, donde el parámetro -genome_size es el tamaño estimado del genoma (visualice la carpeta Referencia para obtener el tamaño del genoma (pista, use el comando: ll -h)).

  • Para obtener el genoma de referencia:
 wget https://www.dropbox.com/s/01qmwi1b5z8ciua/e.coli_reference.fa

Ustedes NO tienen que descargarlo, ya cuentan con el archivo FASTA en: DATA/Ensamblado_DATA/E_coli_K12_Reference_Genome/e.coli_reference.fa

Nuevamente, la idea es que se copien el archivo para sus carpetas personales.

(parados dentro de su carpeta personal)

  • Usando ruta relativa:
cp ../DATA/Ensamblado_DATA/E_coli_K12_PacBio/ecoli_p6_25x.filtered.fastq ./
  • Usando ruta absoluta:
cp /media/Disco10/curso_genomica/estudiantes/DATA/Ensamblado_DATA/E_coli_K12_PacBio/ecoli_p6_25x.filtered.fastq ./

Vamos a realizar el análisis para los ensamblados realizados:

assembly-scan contigs.fa > reporte_Velvet_assembly-scan.txt

assembly-scan ecoli-pacbio.contigs.fasta > reporte_Pacbio_assembly-scan.txt

assembly-scan ecoli-oxford.contigs.fasta > reporte_Nanopore_assemblyscan.txt
  1. ¿Cuál fue la longitud del k-mero utilizada?
  2. ¿Cuál es el tamaño del genoma de referencia utilizado?
  3. ¿Cuál es el N50 del ensamblado? ¿Le parece bueno?
  4. ¿Cuál es el tamaño del contig más largo?
  5. ¿Cuántos scaffold hay en contigs.fa?