Limpieza de fastq

Primero debemos partir de archivos pair-end limpios, podemos usar nextera_cleaner como se explica aquí, pero no es necesario al final unir las secuencias.

Si tenemos varias carpetas con los archivos fastq crudos y comprimidos (.fastq.gz), podemos limpiarlos todos con el script:

$ multiple_fastq_cleaner

Nos genera una nueva carpeta (clean_fastq) ya con los archivos limpios en subcarpetas, ver los detalles aquí. Básicamente, el script realiza el siguiente comando en todos los subdirectorios:

$ cutadapt -a CAAGCAGAAGACGGCATACGAGAT -a CTGTCTCTTATACACATCT -A AATGATACGGCGACCACCGAGATCTACAC  --times 2 -q 30  -m 20 --trim-n -j 8 -o SAMPLE.R1.fastq -p SAMPLE.R2.fastq -u 15 SAMPLE_R1_001.fastq.gz SAMPLE_R2_001.fastq.gz

Ensamble de contigs

El primer paso es ensamblar las secuencias de todas las muestras para obtener contigs generales, esto lo podemos hacer con megahit y los archivos fastq pareados (pair-end) ya sea de solo una muestra, o de varias. Megahit puede utilizar archivos fasta, fastq y fastq.gz (archivos comprimidos).

Una muestra

Si estamos en el directorio con los dos archivos fastq pareados de una sola muestra, el comando es muy sencillo:

$ megahit -1 Sample.R1.fastq -2 Sample.R2.fastq --min-contig-len 1000 -m 0.4 -t 8 -o megahit

Si solo tenemos un archivo fasta o fastq, no pareados:

$ megahit -r Sample.R1.fastq --min-contig-len 1000 -m 0.4 -t 8 -o megahit

Megahit nos produce un archivo con todos los contigs que pudo formar con las secuencias en el subdirectorio: megahit/final.contigs.fa Si éste archivo esta vacío, es probable que no haya encontrado contigs superiores a 1000 bases, por lo que el parámetro --min-contig-len 1000 habría que bajarlo (quizá a 200) o de plano quitarlo.

Varias muestras

Tenemos primero que crear un listado de las muestras que usaremos para cada uno de las pareadas; para esto tenemos que tener todas los archivos en una carpeta o dentro de subdirectorios pero que los archivos pareados tengan algo que los distinga y que sea del mismo tipo siempre. Por ejemplo: Sample01.R1.fastq y Sample01.R2.fastq. 

Con los siguientes comandos podemos crear una variable para que automáticamente genere un listado para cada archivo pareado separado por comas que se encuentren en subcarpetas:

$ R1=$(find . -name '*.R1.*' | tr '\n' ',' | sed s'/.$//; s/\.\///g')

$ R2=$(find . -name '*.R2.*' | tr '\n' ',' | sed s'/.$//; s/\.\///g')

si las muestras pareadas están todas en la misma carpeta y NO en subcarpetas, entonces debemos cambiar un poco la sintaxis:

$ R1=$(ls -1 *R1* | tr '\n' ',' | sed s'/,$//')

$ R2=$(ls -1 *R2* | tr '\n' ',' | sed s'/,$//')

Y entonces estas variables ya las podemos usar en el comando de Megahit:

$ megahit -1 $R1 -2 $R2 --min-contig-len 1000 -m 0.4 -t 8 -o megahit

Dependiendo de la cantidad de metagenomas, este análisis podría durar horas, por lo que es recomendable usar nohup para prevenir que se corte la conexión y se aborte el proceso:

$ nohup megahit -1 $R1 -2 $R2 --min-contig-len 1000 -m 0.4 -t 8 -o megahit &

Mejor aún es correrlo con screen.

Megahit nos produce un archivo con todos los contigs que pudo formar con las secuencias: megahit/final.contigs.fa

Secuencias sin ambigüidades

En los contigs finales, es necesario que no haya ambigüedades, excepto N, en las secuencias pues al parecer Anvio puede presentar errores en el análisis.

Si sabemos que hay ambigüedades en nuestras secuencias fasta, podemos cambiarlas por N, que si es válido:

perl -pe 's/[^AGTC\n]/N/gi unless m/>/;' old.fa > new.fa

Limpieza de headers

Debido a que los contigs (final.contigs.fa) provenientes de megahit tiene un encabezado largo, es necesario limpiarlo para que Anvi'o los pueda utilizar sin problemas, se puede limpiar con un simple comando:

$ awk '/^>/{print ">'contig'_"++i; next}{print}' megahit/final.contigs.fa > contigs.fa

Este comando cambia el header largo a uno con solo "contig" más un número consecutivo deparado por guión bajo (contig_1, conti_2, etc.)

En general para anvi'o, es importante que los archivos fasta:

• No empiecen con un número ni tengan espacios, solo se permiten en los nombres de los archivos:

  • Letras

  • Números

  • Guion bajo _

• Los nombre de las secuencias (headers) sean cortos y sin espacios

• De preferencia que las secuencias estén en una sola línea

• De preferencia que tengan la extensión .fa

Una vez teniendo los contigs de todas las muestras (contigs.fa), empieza el proceso de Anvi'o que puede ser bastante largo (ver el proceso entero aquí). Para simplificar, tenemos un script en bash (metagenomic_workflow) the realiza todos los pasos necesarios esenciales para tener una gráfica interactiva de Anvi'o lista para análisis.

Debemos tener en un directorio solo las secuencias a analizar de los metagenomas y los contigs (contigs.fa), la estructura del directorio tiene que ser similar a la siguiente:

|-contigs.fa

|_Sample01

|  |_Sample01.R1.fastq

|  |_Sample01.R2.fastq

|_Sample02

   |_Sample02.R1.fastq

   |_Sample02.R2.fastq

Debido a la gran carga de trabajo que implica este proceso, puede tardar horas en realizarlo, por lo que se recomienda ejecutarlo por las tardes o en fin de semana.

$ metagenomic_workflow contigs.fa PROYECTO

contigs.fa es el archivo con los contigs creado anteriormente y que le cambiamos la extensión.

PROYECTO es cualquier nombre (corto y sin espacios) que haga referencia a lo estamos analizando.

Este script genera un archivo llamado PROYECTO.yyyy-mm-dd.hh:mm con todos los archivos de salida.

Importante. Como este script puede tardar horas, es mejor correrlo con screen.

Visualización remota

Como normalmente estos análisis se realizan en un servidor, pero la visualización y manipulación se tiene que hacer en un navegador (se recomienda usar Chrome) en nuestra computadora personal, para lograr esto tenemos dos opciones:

1. Bajar los archivos esenciales a nuestra computadora, que además tiene que tener también instalado Anvio, o

2. Lograr que el servidor mande la señal, que es un archivo html, a nuestra computadora. Para activar esto seguir los pasos explicados aquí.