Identificación genomas

Si se necesita identificar a que especie pertenece un genoma desconocido o para confirmar su taxonomía, podemos usar la base de datos GTDB online o en el servidor Biobacter. Para esto tenemos que tener ya un ensamble del genoma, que no tiene que ser muy bueno ni completo.

Ensamble rápido

Un ensamble muy rápido, no muy bueno pero suficiente para identificación se puede hacer con clover en Biobacter a partir de secuencias pareadas ya limpias; si las secuencias fastq aún no están limpias, podemos limpiarlas con nextera_cleaner primero para generar dos nuevos archivos llamados sample.R1.fastq y sample.R2.fastq con los que trabajar:

$ /opt/clover-2.0/clover -k 40 -i1 R1.fastq -i2 R2.fastq -o clover

Calidad del ensamble

También es conveniente checar la calidad del ensamble  para ver que tan completo está y si hay contaminación. Para esto usaremos checkm2 que analiza la presencia y cantidad de genes marcadores de bacterias o arqueas y mediante esto calcular si están todos presentes (100% completo) y si hay varios repetidos, lo que podría implicar contaminación (más de una especie en el ensamble).

$ conda activate checkm2

$ checkm2 predict --threads 8 --input DIR --output-directory checkm2 --database_path /dbs/checkm2/CheckM2_database/uniref100.KO.1.dmnd

En donde DIR es el directorio que contenga los contigs en un archivo multifasta con terminación .fna , en caso que el ensamble sea .fasta, habrá que decirle esto al checkm2 agregando al comando: -x .fasta

Genera un archivo quality_report.tsv con la siguiente info:

Name Completeness Contamination Completeness_Model_Used Translation_Table Used

4A2NS.contigs 87.45 7.14 Neural Network (Specific Model) 11

UJAT02_100K.contigs 95.4 4.86 Gradient Boost (General Model) 11

Identificación rápida

Una buena y rápida clasificación taxonómica de nuestro ensamble la podemos realizar con GSearch en el servidor Biobacter.

Clasificación

$ /opt/gsearch_Linux_x86-64_v0.1.5/gsearch --pio 2000 --nbthreads 24 request -r DIR -n 10 -b /dbs/GTDB-Tk/GTDB/nucl/k16_s12000_n128_ef1600_canonical/

en donde:

DIR es un directorio en donde tengamos el o los genomas a clasificar en formato fasta ya ensamblados.

-b es la ruta a la base de datos

-n número de genomas cercanos a reportar por cada genoma analizado, en este caso 10.

Este análisis nos genera un archivo llamado gsearch.neighbors.txt con varias genomas cercanos con base en su similitud de ANI. Este archivo lo podemos reformatear para tener los valores mas claros:
$ /opt/gsearch_Linux_x86-64_v0.1.5/reformat 16 1 gsearch.neighbors.txt clean_AAI.txt

Reformateo

Ahora generamos un nuevo archivo clean_AAI.txt con los siguientes valores:
Query_Name Distance Neighbor_Fasta_name Neighbor_Seq_Len ANI

El mejor hit será el que tengo el valor de ANI mas alto, normalmente el primer genoma que estará en la segunda linea. 

Query_Name Distance Neighbor_Fasta_name Neighbor_Seq_Len ANI

XIM1831.fasta 0.249833 GCF_003544875.1_genomic.fna.gztotal sequence 4910231 99.03735336271534

Identificación

De este primer hit podemos extraer el nombre del ensamble (Neighbor_Fasta_name, tercer columna: GCF_003544875.1_genomic.fna.gztotal sequence) para busca a que especie corresponde en NCBI:

$ awk 'NR==2 {gsub(/_genomic.*/, "", $3); print $3}' clean_AAI.txt

Así obtendremos sólo el nombre del ensamble: GCF_003544875.1 sin todo el resto.

Con este nombre podemos buscarlo en un archivo de NCBI donde están los nombres de todos los ensambles de genomas:
$ grep "GCF_003544875.1" /dbs/NCBI/assembly_summary_refseq.txt | cut -f8

Así obtendremos el nombre de la especie a la que pertenece, para este ejemplo Vibrio alfacsensis

Identificación mucho más lenta pero más exacta. 

La base de datos y el programa (gtdbtk) para clasificar e identificar está instalada en Biobacter en CONDA, para lo cual hay que activar el ambiente CONDA y luego ejecutar el siguiente comando.

Nota revisar antes si no hay una nueva versión instalada de gtdbtk: conda env list

$ conda activate gtdbtk-2.4.0

$ gtdbtk classify_wf --genome_dir DIR --out_dir DIR/GTDB --cpus 8

DIR es el directorio donde está el ensamble o los contigs en formato fasta con terminación .fna

NOTA. En caso que conda no se active (aparezca un error) es necesario configurar conda como se explica aquí. Si se activó, debería aparecer gtdbtk-2.4.0 antes del prompt $.

Los resultados estarán en el directorio GTDB creado. Para saber que especie es la más cercana podemos ejecutar:
$ cut -f2,6 DIR/GTDB/gtdbtk.bac120.summary.tsv | tail -1 | cut -d ";" -f7 | sed 's/s__//'

y para saber el porcentaje de ANI (si es mayor a 97%, es la misma especie):

$ cut -f2,6 DIR/GTDB/gtdbtk.bac120.summary.tsv | tail -1