Shotgun heces v6.1

Se usarán los datos contenidos en el archivo de heces de camarón (Dietas100K.fastq.tar.gz) que se puede descargar de: 

• Google Drive

• Figshare

Las secuencias ya están limpias y listas para usarse.

Necesitamos primero descomprimir los archivos de las secuencias en todos los subdirectorios:

$ find ./ -type f -name "*.gz" -exec gzip -d "{}" \;

Ensamble de contigs

Tenemos primero que ensamblar contigs con todos los metagenomas, usaremos megahit para este propósito pues está diseñado para metagenomas.

El primer paso es 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')

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

$ megahit -1 $R1 -2 $R2 --min-contig-len 500 -m 0.8 -t 2 -o megahit

Este análisis puede tardar varios minutos dependiendo del número de núcleos (16 min con -t 2) que le hayamos asignado; nos generará una carpeta llamada megahit con archivos, el ensamble estará en final.contigs.fa, podemos inspeccionarlo y ver que tiene 11,734 contigs y una N50 = 929 pb. Alternativamente, para reducir el numero de contigs se puede usar --min-contig-len 1000 pero tener cuidado para que en paso de taxonomía usar el archivo correspondiente.

Podemos ver también que el encabezado (heading) de los contigs es así:

$ grep ">" megahit/final.contigs.fa | head -3

>k99_69 flag=1 multi=3.0000 len=1393

>k99_70 flag=1 multi=2.0000 len=1036

>k99_77 flag=1 multi=3.0000 len=1587

Este tipo de encabezado no le gusta a Anvio, por lo que tenemos que cambiarlo a algo mas sencillo, sin espacios, corto y con números consecutivos:

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

Ahora tendremos los encabezados mas sencillos:

$ grep ">" contigs.fa | head -3

>contig_1

>contig_2

>contig_3

Importante, megahit crea muchos archivos intermedios grandes que solo nos ocupan espacio, y no tenemos mucho en la imagen virtual, borrémoslos, pero antes asegurémonos que tenemos el archivo contigs.fa (generado anteriormente) en la carpeta. Como comprobación de que tanto espacio ocupa, chequemos:

$ du -h megahit/

280M  megahit/intermediate_contigs

285M megahit/

Ahora si borrémoslo:

$ rm -fr megahit/

Mapeo de contigs

Primero debemos activar el ambiente Anvio:

$ conda activate anvio6

Nótese que ahora aparece (anvio6) antes del prompt $.

Teniendo un archivo de contigs formado con todos los metagenomas, debemos ahora mapear cada metagenoma a ese esos contigs, esto nos dirá cuantas y cuales secuencias de cada metagenoma contribuyeron a cada contig, primero tenemos que crear un índice:

(anvio6)$ bowtie2-build contigs.fa contigs

Ahora si mapeamos los cinco metagenomas usando ese índice (contigs):

(anvio6)$ bowtie2 --threads 2 -x contigs -1 C08/C08.R1.fastq -2 C08/C08.R2.fastq -S C08.sam

(anvio6)$ bowtie2 --threads 2 -x contigs -1 P08/P08.R1.fastq -2 P08/P08.R2.fastq -S P08.sam

(anvio6)$ bowtie2 --threads 2 -x contigs -1 P18/P18.R1.fastq -2 P18/P18.R2.fastq -S P18.sam

(anvio6)$ bowtie2 --threads 2 -x contigs -1 S04/S04.R1.fastq -2 S04/S04.R2.fastq -S S04.sam

(anvio6)$ bowtie2 --threads 2 -x contigs -1 S19/S19.R1.fastq -2 S19/S19.R2.fastq -S S19.sam

Convertir los archivos generados a formato sam:

(anvio6)$ for f in *.sam; do samtools view -F 4 -bS $f > $f.raw; done

Ahora debemos inicializar los archivos generado con Anvio:

(anvio6)$ for f in *.raw; do anvi-init-bam $f -o $f.bam; done

Limpieza de archivo ya no necesarios:

(anvio6)$ rm *.raw *.sam

Como los nombres de los archivos van creciendo, es conveniente acortarlos:

(anvio6)$ rename 's/.sam.raw.bam/.bam/' *.bam

(anvio6)$ rename 's/.sam.raw.bam.bai/.bai/' *.bai

Ahora tendremos nombre mas cortos para los archivos .bam y .bai

Creación de base de datos

(anvio6)$ anvi-gen-contigs-database -f contigs.fa -o contigs.db

Posteriormente es buena idea buscar hidden Markov models (hmm) en nuestros contigs, éstos son Single Copy Genes para bacterias y nos ayudarán a saber si un genoma bacteriano encontrado está "completo" o no:

(anvio6)$ anvi-run-hmms -c contigs.db -I Bacteria_71 -T 2

También podemos extraer las secuencias de los genes existentes en la base de datos para una posterior clasificación taxonómica:

(anvio6)$ anvi-get-sequences-for-gene-calls -c contigs.db -o gene-calls.fa

Perfiles para cada muestra

Teniendo ya todos los datos, debemos crear un "profile" para cada muestra, tener cuidado por que el subdirectorio va  a ser borrado y crear uno nuevo, pero de todos modos ya no necesitamos los archivos .fastq y solo nos quitan nuestro preciado espacio.

(anvio6)$ anvi-profile -i C08.bam -c contigs.db --output-dir C08 --sample-name C08 -W -T 2

Esto hay que hacerlo para cada una de las cinco muestras por separado:

(anvio6)$ anvi-profile -i P08.bam -c contigs.db --output-dir P08 --sample-name P08 -W -T 2

(anvio6)$ anvi-profile -i P18.bam -c contigs.db --output-dir P18 --sample-name P18 -W -T 2

(anvio6)$ anvi-profile -i S04.bam -c contigs.db --output-dir S04 --sample-name S04 -W -T 2

(anvio6)$ anvi-profile -i S19.bam -c contigs.db --output-dir S19 --sample-name S19 -W -T 2

Luego tenemos que unir todos los cinco profiles:

(anvio6)$ anvi-merge */PROFILE.db -o SAMPLES-MERGED -c contigs.db --sample-name DIETAS

Limpieza

Ya no necesitamos los archivos .bam ni .bai, ni los creados durante la generación de índices, podemos eliminarlos:

(anvio6)$ rm *.bai *.bam contigs.1.bt2 contigs.2.bt2 contigs.3.bt2 contigs.4.bt2 contigs.rev.*

Taxonomía

Como no corrimos clasificación taxonómica ni funcional, por falta de espacio en la imagen virtual, no podremos obtener información al respecto en la visualización. Los archivos de clasificación se pueden bajar de aquí. Si se bajó el set de datos, tendremos ya una carpeta con ellos (anvio_v6).

(anvio6)$ anvi-import-taxonomy-for-genes -c contigs.db -i gene-calls-500.tax -p kaiju --just-do-it

En caso que hayamos corrido la clasificación taxonómica aparte (ver aquí)  y tengamos el resultado, podemos importarla a la base de datos. El archivo gene-calls_Dietas.tax hay que incorporarlo a la base de datos contigs.db.

Los archivos para taxonomía se pueden descargar de las siguientes rutas y descomprimirlos, usar solo el archivo correspondiente al tamaño de contigs que se seleccionó al crearlos con Megahit:

• Archivo para contigs de 1,000 bases: gene-calls-1000.tax.zip

• Archivo para contigs de 500 bases: gene-calls-500.tax.zip

Descripción

Tenemos también un archivo (Dietas_descripcion.md) con una breve descripción de los datos y métodos, podemos incluirlos en la base de datos para que sea desplegado en la visualización:

(anvio6)$ anvi-update-db-description --description Dietas_descripcion.md SAMPLES-MERGED/PROFILE.db

Metadatos

Por último podemos importar los metadatos a anvio:

(anvio6)$ anvi-import-misc-data metadata.txt -p SAMPLES-MERGED/PROFILE.db --target-data-table layers

Visualización de resultados

(anvio6)$ anvi-interactive -p SAMPLES-MERGED/PROFILE.db -c contigs.db --taxonomic-level t_species

Si todo salió bien, se debe abrir una ventana de Chrome con la siguiente imagen. Una explicación detallada de lo que podemos analizar con esta imagen interactiva se encuentra aqui.

Al terminar, podemos desconectarnos del ambiente virtual Anvio6:

(anvio6)$ conda deactivate