Práctica de metagenómica para el tópico Hackeando las comunidades microbianas impartido en el Instituto de Ciencias del Mar y Limnología (ICMyL) de la UNAM. El material de esta práctica inicialmente se preparó para los Talleres Internacionales de Bioinformática (TIB2022) organizados por la Sociedad Mexicana de Bioinformática (SMB) y la Comunidad de Desarrollo de Software Bioinformático (CDSB). El material se ha ido modificando con el tiempo y sigue en desarrollo.
Autoras: Mirna Rosas Vázquez Landa @MirnaVazquez y Diana Hernández Oaxaca @DianaOaxaca
Sinopsis: En este trabajo se secuenciaron librerias de Illumina MiSeq (2 X 149 pb) de muestras a diferentes tiempos de fermentación: Aguamiel, 0, 3, 6 y 12 horas. En general observaron que la abundancia de los géneros cambió durante la fermentación y se asoció con una disminución de sacarosa, aumento de etanol y ácido láctico (Medidos por HPLC). Predijeron, entre otros, la biosíntesis de folato y vitamina B. Uno de los MAG obtenidos correspondió a Saccharomyces cereviseae relacionado filogenéticamente con S. cerevisiae aislado de sake y bioetanol y otro correspondió a Zymomonas mobilis que propusieron como un nuevo linaje.
Trabajaremos con un subset de los datos originales.
Los datos crudos: PRJNA603591 https://www.ebi.ac.uk/ena/browser/view/PRJNA603591
El articulo: Genomic profiling of bacterial and fungal communities and their predictive functionality during pulque fermentation by whole-genome shotgun sequencing https://www.nature.com/articles/s41598-020-71864-4
Referencia: Chacón-Vargas, K., Torres, J., Giles-Gómez, M. et al. Genomic profiling of bacterial and fungal communities and their predictive functionality during pulque fermentation by whole-genome shotgun sequencing. Sci Rep 10, 15115 (2020). https://doi.org/10.1038/s41598-020-71864-4
# Creamos los directorios en donde empezaremos a trabajar
mkdir -p 03.Análisis_de_metagenomas/{data/{raw,clean},results}
#Accedemos al directorio principal
cd 03.Análisis_de_metagenomas
Una herramienta común para verificar la calidad de las lecturas es FastQC, genera reportes que nos permiten darnos una idea de como están las lecturas, los posibles problemas y que futuros análisis son necesarios.
Para obtener los reportes y el filtrado de calidad ejecutamos la siguiente linea:
mkdir -p results/{01.fastqc,02.trimgalore}
fastqc data/raw/*.fastq -o results/01.fastqc/
conda activate /botete/diana/.conda/envs/qiime2-2022.11
multiqc results/01.fastqc/*.zip -o results/01.fastqc/multiqc
La herramienta TrimGalore nos permite eliminar lecturas de baja calidad, adaptadores, etc. Y con MultiQC, podemos ver las calidades del conjunto de lecturas. Existen otros programas para limpiar las lecturas como Trimmomatic.
#Ejecutamos trim_galore para filtrar
trim_galore --fastqc -j 45 --paired data/raw/pulquet0_1_10M.fastq data/raw/pulquet0_2_10M.fastq -o results/02.trimgalore/pulquet0_trimgalore
Para evaluar los resultados de calidad tras el filtrado, podemos ejecutar multiqc.
Creamos el directorio donde alojaremos las entradas de multiqc y movemos estos archivos al nuevo directorio.
mkdir -p results/02.trimgalore/zips
mv results/02.trimgalore/*/*.zip results/02.trimgalore/zips/
Ahora si ejecutamos multiqc que esta en en el ambiente de qiime2, por eso lo vamos a activar. OJO, multiqc puede estar alojado en un ambiente independiente o instalarse fuera de ambientes, en el servidor esta dentro de qiime pero no es una regla.
#Activamos el ambiente
conda activate /botete/diana/.conda/envs/qiime2-2022.11
#Ejecutamos multiqc
multiqc results/02.trimgalore/zips/*.zip -o results/02.trimgalore/zips/multiqc
Movemos nuestros datos limpios al subdirectorio data/clean
mv results/02.trimgalore/*/*.fq data/clean/
Y descativamos el ambiente
conda deactivate
En este emplo ensamblaremos todas las muestras que dan lugar a la fermentación del pulque, es decir, haremos un coensamble. Sin embargo, también podríamos correr un esamble por muestra de manera independiente.
Ensamblaremos con MetaSPAdes y Megahit, después compararemos ambos coensambles basándonos en longitudes de los contigs con QUAST y el porcentaje de lecturas que dan lugar al coensamble, esto con bbmap de bbtools.
Para acceder a metaspades y megahit activamos el ambiente metagenomics
conda activate metagenomics
Ejecutamos megahit ...
megahit -t 34 --k-list 21,33,55,77,99,127 --min-contig-len 1000 -1 data/clean/fermentation_1.fastq -2 data/clean/fermentation_2.fastq -o results/03.megahit
Para SPAdes necesitamos crear el directorio en el que se alojarán los resultados
mkdir -p results/03.metaspades
Y lo ejecutamos ....
spades.py --meta -k 21,33,55,77,99,127 -t 30 -1 data/clean/fermentation_1.fastq -2 data/clean/fermentation_2.fastq -o results/03.metaspades
Desactivamos el ambiente
conda deactivate
En spades no pudimos delimitar el tamaño mínimo de contig, por lo tanto tenemos que eliminar los de menor tamaño.
Veamos cuanto miden los contigs de nuestro ensamble.
infoseq -auto -nocolumns -delimiter '\t' -only -noheading -name -length results/03.metaspades/contigs.fasta | sort -k2n | cut -f2 | sort | uniq -c | sort -k2n | less
Vamos a filtrarlo ... Para el filtrado necesitamos entrar al directorio de resultados de spades y ejecutar el script desde ahí.
python filter_contig_length_sp.py 1000 1 contigs.fasta
Comparemos el número de contigs que había antes y después de filtrar.
grep -c '>' contigs.fasta
173273
grep -c '>' contigs.fltr.fasta
18673
Creamos el directorio para los mapeos
mkdir -p results/04.depth/
Mapeamos las lecturas al ensamble de metaspades con bbmap
#Hacer ligas simbolicas y correr dentro del directorio 04.depth
nohup bbmap.sh ref=results/03.metaspades/contigsftr.fasta in=data/clean/fermentation_1.fastq in2=data/clean/fermentation_2.fastq out=results/04.depth/fermentation.metaspades.sam scafstats=fermentation_metaspades.scafstats &
Bien, veamos cuantas lecturas mapearon al ensamble
#Número de lecturas iniciales
95588026
#Obtenemos las que mapearon
cut -f8 results/04.depth/fermentation_metaspades.scafstats | grep -v assignedReads | awk '{sum += $1;}END{print sum;}'
91148282
#Sacamos el porcentaje
95.35 %
Hagamos lo mismo para el de megahit
nohup bbmap.sh ref=results/03.megahit/final.contigs.fa in=data/clean/fermentation_1.fastq in2=data/clean/fermentation_2.fastq out=results/04.depth/fermentation.megahit.sam maxindel=80 scafstats=results/04.depth/fermentation_megahit.scafstats &
Cuántas mapearon?
95588026
cut -f8 results/04.depth/fermentation_megahit.scafstats | grep -v assignedReads | awk '{sum += $1;}END{print sum;}'
91218466
95.428
Ok, falta observar las longitudes de los contigs y el tamaño del ensamble
Activemos el ambiente metagenomics para ejecutar QUAST
conda activate metagenomics
Ahora ejecutemos QUAST ...
quast.py results/03.metaspades/contigsftr.fasta results/03.megahit/final.contigs.fa -o results/05.quast
Desactivemos el ambiente
conda deactivate
Y veamos la salida de QUAST para determinar la mejor opción.
Ahora que sabemos que ensamble usar para el binning, empezemos.
Primero vamos a borrar el archivo .sam de megahit porque es muy muy pesado y no lo necesitamos.
rm results/04.depth/fermentation.megahit.sam
Ahora si, para poder agrupar en bins el coensamble, los binneadores requieren de un archivo con las coverturas de cada contig. Afortunadamente ya habíamos mapeado las lecturas, asi que ya tenemos un archivo de mapeo, nos falta convertirlo a formato bam y ordenarlo. Para ello ejecutamos samtools
samtools view -bShu results/04.depth/fermentation.metaspades.sam | samtools sort -@ 80 -o results/04.depth/fermentation_metaspades_sorted.bam
samtools index results/04.depth/fermentation_metaspades_sorted.bam
Y borramos el .sam
rm results/04.depth/fermentation.metaspades.sam
Ahora vamos a obtener la información en el formato que necesitan los programas que harán el binning. Para ello usamos jgi_summarize_bam_contig_depths
.
Activamos el ambiente de metabat
conda activate metabat
Ahora si, obtengamos la información del mapeo
jgi_summarize_bam_contig_depths --outputDepth results/04.depth/fermentation_metaspades_depth.txt results/04.depth/fermentation_metaspades_sorted.bam
Ejecutemos metabat2
metabat2 -i results/03.metaspades/contigsftr.fasta -a results/04.depth/fermentation_metaspades_depth.txt -o results/06.metabat/metabat_bin -t 50 -m 1500
# Y desactivemos el ambiente
conda deactivate
Ahora ejecutemos Maxbin
Activa el ambiente
conda activate maxbin
Creamos el directorio de estos resultados
mkdir -p results/07.maxbin/
Ahora si, corremos Maxbin
run_MaxBin.pl -thread 48 -min_contig_length 1500 -contig results/03.metaspades/contigsftr.fasta -out results/07.maxbin/maxbin -abund results/04.depth/fermentation_metaspades_depth.txt
Y desactivamos el ambiente
conda deactivate
Ya casi....
Corramos Vamb
/botete/diana/.local/bin/vamb
vamb --fasta results/03.metaspades/contigsftr.fasta --jgi results/04.depth/fermentation_metaspades_depth.txt --minfasta 500000 --outdir results/08.vamb
Tenemos los resultados de tres binneadores, muchos de estos bins estarán repetidos, ejecutaremos DAS_Tool para desreplicar estos bins, el flujo para DAS_Tool es el siguiente:
Primero activamos el ambiente.
conda activate das_tool
Creamos un directorio de trabajo para los resultados
mkdir -p results/09.dastool
Y creamos archivos tabulares que son legibles para DAS_Tool
#Para metabat
Fasta_to_Contig2Bin.sh -i results/06.metabat/ -e fa > results/09.dastool/fermentation_metabat.dastool.tsv
#Para maxbin
Fasta_to_Contig2Bin.sh -i results/07.maxbin/ -e fasta > results/09.dastool/fermentation_maxbin.dastool.tsv
#Para vamb
Fasta_to_Contig2Bin.sh -i results/08.vamb/bins/ -e fna > results/09.dastool/fermentation_vamb.dastool.tsv
Ahora si, ejecutemos DAS_Tool:
DAS_Tool -i results/09.dastool/fermentation_maxbin.dastool.tsv,results/09.dastool/fermentation_metabat.dastool.tsv,results/09.dastool/fermentation_vamb.dastool.tsv -l maxbin,metabat,vamb -c results/03.metaspades/contigsftr.fasta -o results/09.dastool/fermentation_bins -t 12 --write_bins
Ya desreplicamos los bins que obtuvimos, pero nos falta evaluar la calidad de estos, para ello ejecutaremos CheckM
checkm lineage_wf results/09.dastool/fermentation_bins_DASTool_bins/ results/10.checkm/ -x fa -t 40 -f results/10.checkm/checkm_dastool_bins.txt
Recordemos que hay más herramientas para desreplicar como dRep y también para refinar los bins que obtengamos, como RefineM. Si te es posible, detente a observar tus datos, analiza los resultados, prueba lo que puedas y toma las mejores decisiones.
Bien, ya hemos visto el flujo de trabajo, revisamos algunas salidas y las tomas de decisiones. Ahora te toca a tí.
¿Cómo serían los resultados si en lugar de un coensamble se analizan las muestras individuales? Averiguémoslo entre todXs.
En equipos, ejecuten el flujo de análisis a partir de los ensamble y su comparación hasta el binning. Tomen sólo las lecturas de una muestra, anotaremos los resultados en una tabla compartida y al final los discutiremos.
-
Elige un equipo y que muestra analizarán
-
Crea los directorios de trabajo, el directorio principal será el del nombre d ela muestra que analizarán
-
Crea una liga simbólica de los datos limpios correspondientes a la muestra que eligieron y ponlos en tu directorio
data
los datos se encuentran en:/botete/diana/Hackeando_las_comunidades_microbianas_v1/03.Analisis_de_metagenomas/data/clean/
-
Para remover los contigs cortos puedes copiar el script
filter_contig_length_sp.py
al directorio de resultados de metaspades, su ubicación es:/botete/diana/Hackeando_las_comunidades_microbianas_v1/03.Analisis_de_metagenomas/src/filter_contig_length_sp.py