Las siguientes estadísticas son utilizadas para determinar la calidad de ensamblaje, que tan fraccionado o continuo quedó.

Note que otros estadísticos, como N75, N90, etc., se calculan de manera similar.

Ejemplo:

Tenemos 10 contigs de distintos tamaño: 10, 20, 30, 40, 50, 60, 70, 80, 90, 100. Tamaño total: 550

La mitad del tamaño es 275, que es menor que la suma de los contigs de tamaño (100 + 90 + 80 + 70). El N50 por lo tanto es 70, y el L50 es 4.


Evaluación de Ensamblajes Illumina

Lo primero será copiar los ensamblajes realizados con abyss-pe en scoville3e6.fcien.edu.uy a su maquina virtual.

  1. Crear un directorio para cada ensamblaje. Ejemplo:
prompt@maquina:$ mkdir abyss1
  1. Ingresar al directorio creado
prompt@maquina:$ cd abyss1
  1. Luego copiar el ensamblaje desde la maquina scoville3e6.fcien.edu.uy a su maquina virtual con scp.
prompt@maquina:$ scp cursogN@scoville3e6.fcien.edu.uy:/media2/cursogN/abyss_k94/* .

(Nota: cambie lo que está en rojo por su grupo y el nombre que le haya asignado a cada ensamblaje)

Repita para los puntos a, b y c para los dos ensamblajes que realizó con abyss-pe


1- listar los archivos de cada ensamblaje con ls -l. En particular, los archivos coverage.hist, prefijo-3.hist y prefijo-6.hist, contienen información de interés. Explore su contenido.

2- Determine la profundidad de los ensamblajes que realizó. Calcule en base a la cantidad y largo de los reads usados y el largo estimado del genoma de S. cerevisiae la profundidad esperada para cada ensamblaje.

3- Determine la cantidad de contigs que obtuvieron sus ensamblajes. Nota: puede usar grep y recuerde que el nombre del resultado final del ensamblaje es del tipo ABYSS-8.fa (es igual a ABYSS-scaffolds.fa), o más en general, prefijo-8.fa o prefijo-scaffolds.fa

4- El programa infoseq del paquete EMBOSS nos brinda información sobre un archivo de secuencias fasta

Corra el programa infoseq para alguno de sus resultados de ensamblaje con abyss-pe

prompt@maquina:$ infoseq ABYSS-scaffolds.fa

Identifique las columnas (puede valerse del comando head) y vea las opciones del programa con infoseq -help

En especial nos interesa solamente el nombre de los contigs, el largo y el contenido de GC, podemos obtener un versión resumida de la tabla con la opcion -only. Ejemplo:

prompt@maquina:$ infoseq -only -name -length -pgc ABYSS-scaffolds.fa > tabla.out

Usando el programa infoseq, determinar para cada ensamblaje los siguientes parámetros:

  1. El contig más largo. Puede usar sort con la opción -k2 o -k6 (dependiendo de si usan una tabla resumida o no) para ordenar por dicha fila, y la opción -n para ordenar numéricamente. Ejemplo:
prompt@maquina:$ sort -n -k2 tabla.out

  1. El número de contigs mayores a 1 kb (puede usar awk para filtrar por la columna 2 la tabla generada)

  2. El largo total del ensamblaje.

Hay varias formas de obtener el largo total del ensamblaje. Una forma es cargar tabla.out en una hoja de cálculo como excel y sumar la columna de los largos. Otra forma es usar awk, filtrar por la columna de los largos y sumar (con una sola línea de comando se obtiene el largo total!).

prompt@maquina:$ awk '{largo+= $2 } END {print largo}' tabla.out

Si quisiera obtener el largo de la suma de los contigs mayores a 5Kb, solo debería agregar esa condición en el awk. ¿Cuánto es la suma de los contigs mayores a 5Kb?

5- Calcule el N50 y el L50 para alguno de sus ensamblajes abyss-pe. Para esto puede abrir el archivo tabla.out, generado a partir de correr infoseq, en excel u otra hoja de cálculo. Recuerde que deberá ordenar los contigs de mayor a menor y realizar la suma acumulada. Aquellos que prefieran lo pueden hacer usando awk.

6- Obtenga las estadísticas de los tres ensamblajes con el programa abyss-fac

prompt@maquina:$ /usr/lib/abyss/abyss-fac ABYSS-scaffolds.fa # use el nombre que corresponda

Cuál de los ensamblajes obtenidos usando ABYSS cree que es de mayor calidad? ¿A qué atribuye la diferencia de calidad? ¿Está conforme con el resultado? ¿Como se podría mejorar este resultado?

A partir de ahora se va a trabajar SOLO con el mejor ensamblaje obtenido.

7- Con la ayuda del programa infoseq, grep y awk, obtenga un archivo multifasta que contenga solamente los contigs superiores a 1 kb del ensamblaje seleccionado anteriormente (es decir, un sub-ensamblaje)

prompt@maquina:$ awk '$2>1000 {print ">"$1}' tabla.out > nombres.txt # Note como con awk se puede agregar texto!

Visualizar el archivo nombres.txt y corregir si identifica líneas no deseadas.

prompt@maquina:$ grep -w -A1 -f nombres.txt ABYSS-MEJORENSAMBLAJE.fa > contigs1Kb_abyss.fas

La opción -A1 permite seleccionar la siguiente línea, y la opción -f nos permite dar como entrada un archivo con todos los nombres a buscar.

Observar y corroborar el archivo multifasta generado. ¿Recuerda cómo eliminar las líneas con guiones no deseadas?

8- Calcule el largo total de este genoma filtrado, la cantidad de contigs, el N50 y el L50. Discuta los resultados comparando con el ensamblaje del cual se originó.

9- ¿Hay algo a destacar en relación a los largos o contenidos GC? ¿Alguno de los contigs presenta comportamiento anómalo?

10- Compararemos ahora nuestro ensamblaje (o sub-ensamblaje) con el genoma de referencia de la cepa S288c, para lo cual debemos mapear los contigs de este sub-ensamblaje (obtenido en el punto 7), usando BLASTN contra los cromosomas de levadura cepa S288c. Este genoma está localizado en /media2/YEAST/genoma_ncbi. El archivo se llama YEAST_chromosome.fas.

prompt@maquina:$ blastn -query contigs1Kb_abyss.fas -subject /media2/YEAST/genoma_ncbi/YEAST_chromosome.fas -outfmt "6 std qlen slen" -evalue 1e-200 > tabla_blast

  1. Determinar en cuántos contigs quedó fragmentado cada cromosoma (recuerde que vimos cómo hacerlo en el primer práctico: tomar las columnas 1 y 2 con awk, y hacer un pipeline con los comandos grep, sort y uniq -c).

  2. ¿Son todos los hits de blast igualmente informativos?

  3. ¿Qué criterio usaría para filtrar los hits de blast?

11- Ahora vamos a centrarnos en la comparación de un solo cromosoma. Elija uno de los cromosomas del genoma con longitud mayor a 600 KB, seleccionar y aislar los contigs de su ensamblaje que mapean sobre dicho cromosoma generando un archivo multifasta que contenga solamente estos contigs.

Primero generamos el archivo con los nombres de los contigs que pegan con el cromosoma elegido:

prompt@maquina:$ grep CROMOSOMA tabla_blast | awk '{print ">"$1}' > contigs.lista

Donde CROMOSOMA es el identificador del cromosoma que usamos para buscar con grep. Por ejemplo para el cromosoma XV podríamos usar “C_XV|”, o “C_II|” para el cromosoma II. En ese caso el comando sería:

prompt@maquina:$ grep "C_XV|" tabla_blast | awk '{print ">"$1}' > contig.lista # visualice el resultado contigs.lista con more o less

Nota: es necesario incluir la barra vertical para no “encontrar” otros cromosomas. Por ejemplo, si queremos buscar las coincidencias con el cromosoma XV, debemos cuidar de que nuestro patrón sólo coincida con “C_XV” y no con “C_XVI”, “C_XVII”, etc.

Luego extraemos los contigs (del archivo de salida de ABYSS) cuyos nombres están en contig.lista con el comando grep, nuevamente usando las opciones -A1 y -f.

prompt@maquina:$ grep -A1 -w -f contig.lista ABYSS-scaffolds.fa > contigs_cromosoma.fas

También extraemos la secuencia del cromosoma elegido del genoma de referencia:

prompt@maquina:$ grep -A1 CROMOSOMA /media2/YEAST/genoma_ncbi/YEAST_chromosomeU.fas > chromosome.fasta

Por ejemplo, para el cromosoma XV, el comando sería:

prompt@maquina:$ grep -A1 "C_XV|" /media2/YEAST/genoma_ncbi/YEAST_chromosomeU.fas  > chromosome.fasta

Nota: usamos el archivo YEAST_chromosomeU.fas en lugar del YEAST_chromosome.fas porque el primero contiene cada secuencia en una única línea, lo que permite usar grep -A1.

12- Una forma práctica de visualizar la salida de blast en forma gráfica es utilizar NCBI-blastn. Puede buscar NCBI blast en su navegador (ATENCIÓN: usar el navegador de la máquina virtual) o utilizar el siguiente link: https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome

Clickear en la opción “Align two or more sequences”, y subir sus archivos multifasta pinchando en Browse. Use como query la secuencia completa del cromosoma y como subject los contigs de su ensamblado que pegaron contra ese cromosoma.


Identifique las zonas de “quiebre” del ensamblaje, y determine qué características tienen estas zonas (repetición, contenido G+C).

13- Realizar un dot-plot de estos contigs versus el cromosoma elegido.

Los dot-plot los realizaremos en el sitio web YASS. Puede buscar YASS genome en su navegador o utilizar el link: https://bioinfo.lifl.fr/yass/yass.php

YASS es de uso fácil e intuitivo, se suben los archivos a comparar o directamente copy-and-paste en la ventana de entrada de secuencias. Recuerde usar “select” para que el archivo o la secuencia queden seleccionados. Use “run YASS” para correr el programa.

  1. Compare primero el cromosoma elegido consigo mismo.

  2. En una nueva ventana compare el multifasta generado con los contigs de ABYSS vs el mencionado cromosoma. Observe en el dot-plot las zonas ruptura del ensamblaje.