Clasificación taxonómica
Para clasificar taxonómicamente las secuencias presentes en un metagenoma, podemos hacerlo de sod formas, 1) comparando las secuencias contra genomas conocidos y 2) extrayendo las secuencias de genes ribosomales (16S y 18S) y comparándolas contra una base de datos.
Una pregunta frecuente es ¿Qué tanto tengo que secuenciar para tener una buena representatividad de la diversidad en mi muestra? Aquí muestro un análisis de que se obtiene secuenciando a diversas profundidades.
Clasificación comparando contra genomas conocidos
Kaiju
Kaiju es quizá el mejor clasificador taxonómico para secuencias de metagenomas. Si bien requiere varios pasos para obtener un análisis completo, estos pasos no son difíciles. Es necesario tener una base de datos de los genomas microbianos para que compare, lo cual sólo se puede hacer en el servidor por lo altos requerimientos computacionales.
Lo primero es comparar las secuencias (fasta o fastq) con la base de datos, tenemos instaladas las siguientes bases de datos:
ProGenomes. 141 millones de secuencias de genomas bacterianos y de arqueas.
refseq. 127 millones de secuencias de genomas bacterianos, arqueas y virus.
nr+euk. 321 millones de secuencias de genomas bacterianos, arqueas, virus y protistas.
Las rutas en Biobacter donde se encuentran las bases de datos son: /db2/kaiju_dbs/
Análisis con Kaiju
$ kaiju -t /dbs/kaiju_dbs/progenomes_2023-05-10/nodes.dmp -f /dbs/kaiju_dbs/progenomes_2023-05-10/kaiju_db_progenomes.fmi -i File.R1.fastq.gz -j File.R2.fastq.gz -o kaiju.out -v -z 4
-f y -t son los archivos de la base de datos que están en /db2/kaiju_dbs/ en Biobacter, ver ruta completa arriba
-i es el archivo de entrada, puede ser la R1.fastq de archivos pair-end o solo uno en fasta o fastq
-j el segundo archivo pair-end, R2.fastq (opcional)
-v verbose, para que también adjunte a la salida información adicional
-z el número de núcleos a usar
-o el archivo de salida
A continuación es conveniente generar un archivo tsv pero con la los taxones
$ /opt/kaiju-v1.9.0-linux-x86_64-static/kaiju2table -t /dbs/kaiju_dbs/progenomes_2023-05-10/nodes.dmp -n /dbs/kaiju_dbs/progenomes_2023-05-10/names.dmp -r species -l superkingdom,phylum,class,order,family,genus,species -o kaiju_summary.tsv kaiju.out
Este script toma la salida del primero (kaiju.out)y genera un reporte con porcentaje, número de secuencias y taxonomía asignada.
También es posible generar un gráfico Krona:
$ /opt/kaiju-v1.9.0-linux-x86_64-static/kaiju2krona -t /dbs/kaiju_dbs/progenomes_2023-05-10/nodes.dmp -n /dbs/kaiju_dbs/progenomes_2023-05-10/names.dmp -i kaiju.out -v -o kaiju2krona.out
que hay que importar a krona para generar el html:
ktImportText -o kaiju.html kaiju.krona
Si se tienen varias muestras, hay que hacer el kaiju para cada una y luego juntarlas, ver aqui o preferentemente correr el script descrito abajo para múltiples muestras.
En Biobacter tenemos ya un script para simplificar todos estos pasos y generar la clasificación por muestra:
$ kaiju_file_classifier
Puede usar un archivo fasta o fastq, pero solo uno y solo de una muestra.
Para analizar varias muestras usando ambos archivos fastq/fastq.gz tenemos el script:
$ multiple_kaiju_classifier
Ambos archivos de los metagenomas (.R1.fastq y .R2.fastq) deben estar dentro de subdiretorios en la carpeta donde lo vayamos a correr y deben tener un nombre seguido de R1 o R2 y la terminación fastq o fast.gz (p. ej. muestra1.R1.fastq muestra1.R2.fastq.gz).
Este script genera una tabla de OTUs que puede ser analizada o graficada como se describe aquí.
O bien un script para analizar sin seleccionar nada un par de archivos (R1.fastq y R2.fastq) de una muestra, que se clasificarán con la base de datos más completa (nr+euk)
$ simple_kaiju_classifier sample.R1.fastq sample.R2.fastq
KRAKEN2
Kraken nos permite clasificar las secuencias dependiendo de la frecuencia de kmers que tengan y que tanto son semejantes a una bases de datos.
$ kraken2 --db /dbs/kraken/standard_20220926/ -report SAMPLE.report --use-names --use-mpa-style -o SAMPLE.kraken2 --threads 8 FILE.fasta/q
--db ruta para la base de datos a usar; la ruta del ejemplo es la que tenemos en biobacter. Tenemos dos bases de datos, una con datos de archaea, bacteria, viral, plasmid, human1, UniVec_Core (/dbs/kraken/standard_20220926/) y la otra, además de lo anterior, también protozoa y fungi (/dbs/kraken/pluspf_20230314/).
--threads Número de procesadores a usar para el análisis; el número depende del servidor o computadora.
FILE.fasta/q Archivo con las secuencias a analizar en formato fasta o fastq. Puede utilizar fastq comprimidos (.fastq.gz) pero habría que poner el flag --gzip-compressed.
NOTA. Se recomienda no correr mas de proceso con Kraken2 a la vez pues consume mucha memoria y podría colapsar el servidor; antes de iniciar checar con htop si no se están corriendo otros procesos que estén consumiendo mas del 40% de la memoria.
Script para clasificar metagenomas con Kraken2 en BioBacter:
$ kraken2_classifier *.fasta
Los archivos de los metagenomas deben estar en formato fasta en la misma carpeta y permite seleccionar entre dos bases de datos: Standard (archaea, bacteria, viral, plasmid, human1, UniVec_Core) y refseq (Standard + protozoa, fungi).
FOCUS
FOCUS: An Alignment-Free Model To Identify Organisms In Metagenomes Using Non-Negative Least Squares
USAGE:
python focus.py -q query_sequence.fna [-k] [-m] [-s]
-q Specify input file or directory with multiples samples. Required: Input files should be in FASTA format
-k Specify k-mer frequency used on the profile (default: 7) 6,7 and 8 frequencies are available.
All the output files have this project name as prefix.
-m minimum relative abundance to show in the results (default: 1%) Cut-off for the showing results.
-d Insert your own data into the database more information README .
-s output in the STAMP format for multiple files [0=False (default) and 1=True]
-l Split STAMP output in different levels (default: all; options: kingdom, phylum, class, order, family, genus, or species)
example 1 > (single file) python focus.py -q query.fasta
example 2 > (multiple files) python focus.py -q samples directory -s 1
Extracción de genes ribosomales (16S/18S)
Metaxa2
Metaxa2 es otro clasificador que extrae secuencias de algunas de las subunidades de los genes ribosomales, SSU (16S/18S) y LSU (23S/28S), y las clasifica de acuerdo a bases de datos propias de arqueas, bacterias, cloroplastos y eucariotas.
$ metaxa2 -i file.fasta -o metaxa --cpu 8
Si tenemos archivos pair-end fastq, podemos ejecutarlo así:
$ metaxa2 -1 R1.fastq -2 R2.fastq -o metaxa --cpu 8
Ambas opciones pueden durar horas, por lo que es conveniente correrlo en screen.
Otras opciones son:
-i {file} : DNA FASTA input file to investigate.
-o {file} : Base for the names of output file(s).
-f {a, auto, f, fasta, q, fastq, p, paired-end} : Specifies the format of the input file, default = auto
-g {ssu, lsu} : Specifies if Metaxa should look for SSU or LSU rRNA genes, default is ssu.
--cpu {value} : The number of CPU threads to use, default is 1.
--help : displays the full options for Metaxa.
-p {directory} : A path to a directory of HMM-profile collections representing rRNA conserved regions, default is in the same directory as metaxa itself.
-d {database} : The BLAST databased used for classification, default is in the same directory as metaxa itself.
-t {b, bacteria, a, archaea, e, eukaryota, m, mitochondrial, c, chloroplast, A, all, o, other} : Profile set to use for the search (comma-separated), default is all.
Ribotagger
$ perl /opt/ribotagger/ribotagger.pl -region v4 -out outfile -in file1.fasta file2.fasta ...
Options
-region Can be any of the variable regions V4 to V7 (v4 v5 v6 v7)
-out Name of an output file
-in fasta file(s)