Shotgun heces v7

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

Ensamble de contigs

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')

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

$ megahit -1 $R1 -2 $R2 --min-contig-len 1000 -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 y la cantidad de memoria asignada -m 0.8 80%) 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 2,962 contigs y una N50 = 1,614 pb. 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 anvio-7

Nótese que ahora aparece (anvio-7) 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:

(anvio-7)$ bowtie2-build contigs.fa contigs

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

(anvio-7)$ bowtie2 --threads 2 -x contigs -f sample1.fna -S sample1.sam

Esto tenemos que hacer para cada uno de las cinco metagenomas:

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

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

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

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

(anvio-7)$ 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:

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

Ahora debemos inicializar los archivos generado con Anvio:

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

Limpieza de archivo ya no necesarios:

(anvio-7)$ rm *.raw *.sam

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

(anvio-7)$ rename 's/.sam.raw.bam/.bam/' *.bam

(anvio-7)$ 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

(anvio-7)$ 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:

(anvio-7)$ 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:

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

Profile 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.

(anvio-7)$ 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:

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

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

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

(anvio-7)$ 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:

(anvio-7)$ 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:

(anvio-7)$ 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.

En caso que hallamos corrido la clasificación taxonómica aparte (ver aquí) y tengamos el resultado, podemos importarla a la base de datos.

Descargar el archivo: gene-calls_anvio7.zip Este archivo comprimido contiene cuatro archivos:

1. gene-calls.fa Archivo fasta de los genes encontrados en los contigs.

2. gene-calls.tax Taxonomía asociada a los genes anotados.

3. gene-calls.out Anotación de los genes sin taxonomía.

4. kaiju.krona Archivo para generación de gráfica Krona.

El archivo gene-calls_anvio7.tax hay que incorporarlo a la base de datos contigs.db:

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

Podemos generar una gráfica Krona con el archivo kaiju.krona

(anvio-7)$ ktImportText -o kaiju_gene_taxonomy.html kaiju.krona

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:

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

Metadatos

Por último podemos importar los metadatos a anvio:

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

Visualización de resultados

(anvio-7)$ 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:

Al terminar, podemos desconectarnos del ambiente virtual anvio-7:

(anvio-7)$ conda deactivate