SITUACIÓN DONDE ESTA GUÍA ES ÚTIL: Tienes reads forward y reverse limpios (sin adaptadores) de una secuenciación masiva de amplicones ITS (ITS metabarcoding), esta guía te servirá para realizar un análisis de taxonomía usando el protocolo Dada2 en Qiime2.

Si tienes reads forward y reverse limpios (sin adaptadores) de una secuenciación masiva de amplicones r16S, debes leer esta otra guía.

Si tus reads ya están unidos (típico del postproceso de MrDNA o BGI) y quieres realizar em Qiime2 un análisis similar a Qiime1 puedes adaptar para ITS el protocolo VSearch que detallo en esta guía.

Esta guía se basa en un análisis de secuencias de ITS usando el plugin Q2-ITSxpress. Puedes saltarse este paso para lo cual recomiendo seguir el protocolo de la guía de r16S.

No es una guía de conceptos y procedimientos
Sólo son los comandos con algunas explicaciones útiles para quienes saben del tema o tiene un jomaldon para guiarlos.

Adaptado desde este sitio web

Glosario de extensiones de archivos

• Archivo .gz = Cualquier archivo comprimido con el programa gzip
• Archivo .tar =Paquete de archivos (Un archivo que en su interior tiene varios archivos y carpetas) pero no está comprimido.
• Archivo .tar.gz o .tgz = Es un archivo .tar comprimido con gzip.
• Archivos .fastq o .fq = Archivos de secuencias en un formato que almacena secuencia y calidad.
• Archivos .fasta o .fa = Archivos de secuencia en un formato que no almacena calidad.
• Archivos .qza: Archivos comprimidos (.zip) que contienen toda la info que sale desde un determinado proceso de Qiime2. Se puede descomprimir (winzip, winrar).
• Archivos .qzv: Archivo comprimido (.zip) que contiene información para "visualizar" los resultados de un archivo .qza. Se puede abrir arrastrándolo a la interfáz web view.qiime2.org o, se puede descomprimir (winzip, winrar) y revisar su contenido. De hecho, dentro de la carpeta "data" encontrarás un archivo index.html el cual puedes abrir y se desplegará en tu navegador web toda la info y gráficos correspondientes al tipo de .qza de origen.
• Archivos .tsv: Archivos de texto (.txt) donde las columnas está separadas por tabuladores. Se puede abrir en Excel como un archivo .csv.

1. Archivos necesarios para trabajar en Qiime2

1.1 Archivos de secuencias crudas en formato fastq o fastq.gz

1.2 Archivo manifest.txt (crear)

Debe tener un encabezado específico (ver ejemplo) y tantas filas como archivos de secuencias tengas
En la primera columna va el nombre de la muestra, luego el nombre del archivo y luego la dirección de secuenciación (cada dato separado por una coma). El archivo manifest.txt debes colocarlo en el mismo directorio donde están las secuencias y debes ejecutar los comandos de qiime2 en el mismo directorio (a menos que entiendas la forma de localizar y escribir el path hacia tus archivos.
Ejemplo:

sample-id,absolute-filepath,direction
C-400-r1,$PWD/Sample1_1.fq.gz,forward
C-400-r1,$PWD/Sample1_2.fq.gz,reverse

Si sólo posees secuencias forward (o secuencias forward y reverse que ya han sido unidas), entonces no necesitas usar la línea "reverse", todas tus muestras deben quedar con dirección "forward".

1.3 Archivo metadata.txt (crear) [antes llamado mapping.txt]

Debe tener un encabezado específico (ver ejemplo) y tantas filas como muestras
En la primera columna va el nombre de la muestra y en las siguientes, los datos "extra", por ejemplo, el agrupamiento (cada dato separado por un tabulador)
Ejemplo:

ID ID2 Species Altitude Replicate
C-400-r1 C-400 C 400 1
C-1000-r1 C-1000 C 1000 2

2. Manos a la obra

2.1 Import sequence data

# Si los datos provienen en dos archivos por muestra (uno para forward y otro para reverse):
qiime tools import --type SampleData[PairedEndSequencesWithQuality] --input-format PairedEndFastqManifestPhred33 --input-path manifest.txt --output-path sequences.qza

# Si los datos provienen de un archivo single o de un archivo forward y reverse ya fusionados (joined):
qiime tools import --type SampleData[SequencesWithQuality] --input-format SingleEndFastqManifestPhred33 --input-path manifest.txt --output-path sequences.qza

# Resumen del resultado de importar las secuencias
qiime demux summarize --i-data sequences.qza --o-visualization sequences.qzv

2.2 Trimming ITS samples with Q2-ITSxpress for Dada2 and QIIME2

Ojo con el número de hilos del procesador... yo he usado 64 pero debes ajustarlo a la capacidad de tu máquina o eliminar el comando para que sólo se use un hilo.

qiime itsxpress trim-pair-output-unmerged --i-per-sample-sequences sequences.qza --p-region ITS2 --p-taxa F --o-trimmed trimmed.qza --p-threads 64

#No realizar
#qiime itsxpress trim-pair-output-unmerged --i-per-sample-sequences sequences.qza --p-region ITS2 --p-taxa F --o-trimmed trimmed_exact.qza --p-threads 64 --p-cluster-id 1.0

2.3 Use Dada2 to identify sequence variants

Los puntos de corte de largo de secuencias en este ejemplo fueron deteminados mirando la calidad de las secuencias (revisando el archivo sequences.qzv).

Ojo con el número de hilos del procesador... yo he usado 64 pero debes ajustarlo a la capacidad de tu máquina o eliminar el comando para que sólo se use un hilo.

El parámetro p-trun-len debe estar definidos según la revisión de calidad de las secuencias en el archivo sequences.qzv. En caso de no desear establecer este parámetro basta con reemplazarlo por cero [0].

El parámetro p-trunc-q es el equivalente al punto de corte de qvalue. "Reads are truncated at the first instance of a quality score less than or equal to this value. If the resulting read is then shorter than trunc-len, it is discarded. default=2"

El parámetro p-max-ee: "Reads with number of expected errors higher than this value will be discarded. default=2"

# Método paired (alinea los pares para usar la info de secuencia más larga)
qiime dada2 denoise-paired --i-demultiplexed-seqs trimmed.qza --p-trunc-len-f 270 --p-trunc-len-r 220 --p-max-ee-f 2 --p-max-ee-r 2 --p-trunc-q 2 --output-dir dada2out --p-n-threads 64

# Método single (usa las secuencias sin aparearlas, es decir, cortas)
qiime dada2 denoise-single --i-demultiplexed-seqs trimmed.qza --p-trunc-len 220 --p-max-ee-r 2  --p-trunc-q 2 --output-dir dada2out --p-n-threads 64

# Visualización del resultado
qiime feature-table summarize --i-table dada2out/table.qza --o-visualization tableviz.qzv

# Visualización del resultado añadiendo info de grupos
qiime feature-table summarize --i-table dada2out/table.qza --o-visualization tableviz2.qzv --m-sample-metadata-file metadata.txt

# Tabla con información de las secuencias representativas (no hacer)
table tabulate-seqs --i-data dada2out/rep-seqs.qza --o-visualization rep-seqs.qza.qzv

2.4 Import the latest UNITE data into QIIME2

Para realizar la clasificación taxonómica es necesario un artefacto del tipo "classifier" que se prepara a partir con las secuencias de la base de datos deseada y su respectivo archivo de texto con la "clasificación" taxonómica para cada ID.
Acá usaré UNITE82s, versión sólo hongos.

# Importamos el archivo fasta
qiime tools import --type 'FeatureData[Sequence]' --input-path sh_refs_qiime_ver8_dynamic_04.02.2020.fasta --output-path unite.qza

# Importamos el archivo de taxonomías que hace match con el archivo fasta anterior
qiime tools import --type 'FeatureData[Taxonomy]' --input-format HeaderlessTSVTaxonomyFormat --input-path sh_taxonomy_qiime_ver8_dynamic_04.02.2020.txt --output-path unite-taxonomy.qza

2.5 Train the QIIME classifier (este paso puede demorar varias horas)

Vamos a "entrenar" un clasificador de taxonomía, lo cual significa que Qiime2 usará un algoritmo (a alección) para "aprender" a clasificar según el marcador molecular que estemos usando. En esta caso usaremos el método Naive-Bayes. Mas info en este link.

qiime feature-classifier fit-classifier-naive-bayes --i-reference-reads unite.qza --i-reference-taxonomy unite-taxonomy.qza --o-classifier classifier.qza

2.6 Classify the sequence variants

Ahora estamos listos para "clasificar" taxonómicamente nuestos ASVs

Ojo con el número de trabajos a lanzar... yo he usado la opción -1 que significa que se usarán todas las CPUs, pero debes ajustarlo a la capacidad de tu máquina o eliminar el comando para que sólo se use un hilo. Se puede colocar un número negativo (ej -64) ante lo cual amplicará lo siguiente (n_cpus_total + 1 - 64).

qiime feature-classifier classify-sklearn --i-classifier classifier.qza --i-reads dada2out/representative_sequences.qza --o-classification taxonomy.qza --p-n-jobs -1

Summarize the results

qiime metadata tabulate --m-input-file taxonomy.qza --o-visualization taxonomy.qzv

2.7 Create an interactive bar plot figure (each taxonomic level exploration)

El siguiente comando permite crear un gráfico interactivo para visualizar las abundancias de cada taxa, en cada muestra mediante barras apiladas. Además, podrás exportar un archivo csv de las abundancias a "cada nivel taxonómico", muy útil para exploración manual de datos en Excel.

qiime taxa barplot --i-table dada2out/table.qza --i-taxonomy taxonomy.qza --m-metadata-file metadata.txt --o-visualization taxa-bar-plots.qzv

A partir de la tabla de abundancias (dada2out/table.qza) y la tabla de taxonomías (taxonomy.qza) se pueden realizar muchos análisis en qiime2, por ahora llego hasta acá pero mas adelante iré completando esta lista de comandos, como la estimación de la diversidad alfa y la diversidad beta con rarefacción. De todas formas, mas abajo te explico como obtener un archivo BIOM que te permitirá seguir con los análisis en el programa MEGAN.

3. Análisis de Diversidad

3.1 Alineamiento múltiple de secuencias representativas

# En un paso (recomendado)
qiime phylogeny align-to-tree-mafft-fasttree \
        --p-n-threads 64 \
        --i-sequences dada2out/representative_sequences.qza \
        --o-alignment aligned-representative_sequences.qza \
        --o-masked-alignment masked-aligned-representative_sequences.qza \
        --o-tree unrooted-tree.qza \
        --o-rooted-tree rooted-tree.qza

# En varios pasos (alternativa)
qiime alignment mafft \
  --i-sequences dada2out/representative_sequences.qza \
  --p-n-threads 64 \
  --o-alignment aligned-representative_sequences.qza

qiime alignment mask \
  --i-alignment aligned-representative_sequences.qza \
  --o-masked-alignment masked-aligned-representative_sequences.qza

qiime phylogeny fasttree \
  --i-alignment masked-aligned-representative_sequences.qza \
  --p-n-threads 64 \
  --o-tree unrooted-tree.qza

qiime phylogeny midpoint-root \
  --i-tree unrooted-tree.qza \
  --o-rooted-tree rooted-tree.qza

3.2 Alpha diversity with rarefaction

Para determinar el punto de corte de la rarefacción tendrás que revisar cuál es la muestra que tiene el menor número de reads (ver archivo tableviz.qzv) y fijar el punto de corte un poco mas abajo de ese número. Por ejemplo, en mi caso, la muestra con menos reads es 16S-T1 con 7252 reads, entoces mi punto de corte será 7000.

qiime diversity alpha-rarefaction \
  --i-phylogeny rooted-tree.qza \
  --i-table dada2out/table.qza \
  --p-max-depth 7000 \
  --p-metrics observed_features chao1 shannon faith_pd goods_coverage \
  --m-metadata-file metadata.txt \
  --o-visualization rarefaction_r7k

# Las opciones de índices son:
'berger_parker_d', 'brillouin_d', 'enspie', 'ace', 'fisher_alpha', 'simpson_e', 'shannon', 'heip_e', 'simpson', 'dominance', 'gini_index', 'mcintosh_d', 'doubles', 'singles', 'margalef', 'goods_coverage', 'robbins', 'lladser_pe', 'menhinick', 'chao1', 'mcintosh_e', 'faith_pd', 'pielou_e', 'michaelis_menten_fit', 'observed_features'

3.3 Core metrics (alpha and beta diversity) with rarefaction

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table dada2out/table.qza \
  --p-sampling-depth 7000 \
  --p-n-jobs-or-threads 64 \
  --m-metadata-file metadata.txt \
  --output-dir core-metrics-phylogenetic-metrics_r7k

3.4 Obtain the rarefaction table (tsv and biom format)

Para este paso primero se debe hacer el paso 4.1 para obtener el archivo taxonomy2.tsv

cd core-metrics-phylogenetic-metrics_r7k

qiime tools export --input-path rarefied_table.qza --output-path ../exported_r7k

cd ../exported_r7k

biom add-metadata -i feature-table.biom -o feature-table-with-taxonomy.biom --observation-metadata-fp ../exported/taxonomy2.tsv --sc-separated taxonomy --observation-header OTUID,taxonomy

# Ahora exportamos la tabla que podremos revisar en Excel
biom convert -i feature-table-with-taxonomy.biom -o feature-table-with-taxonomy.tsv --header-key taxonomy --output-metadata-id "taxonomy" --to-tsv

# Ahora exportamos la tabla que podremos revisar en MEGAN (continuar en sección 4.3)
biom convert -i feature-table-with-taxonomy.tsv -o feature-table-with-taxonomy.biom --table-type="OTU table" --process-obs-metadata='taxonomy' --to-hdf5

3.5 Create an interactive bar plot figure (each taxonomic level exploration)

El siguiente comando permite crear un gráfico interactivo para visualizar las abundancias de cada taxa, en cada muestra mediante barras apiladas. Además, podrás exportar un archivo csv de las abundancias a "cada nivel taxonómico", muy útil para exploración manual de datos en Excel.

qiime taxa barplot --i-table core-metrics-phylogenetic-metrics/rarefied_table.qza --i-taxonomy ../taxonomy.qza --m-metadata-file metadata.txt --o-visualization taxa-bar-plots_r7k.qzv

4. Obtener archivos BIOM (MEGAN) y tsv (Excel) que incluyen abundancia y taxonomía

Qiime2 no tiene un "artefacto" o tipo de archivo que mantenga la información de abundancias y taxonomías juntas. Como ya hemos visto, esta información quedó en dada2out/table.qza y en taxonomy.qza y no hay forma de juntarlos en un mismo qza.

Los siguientes pasos tiene como objetivo obtener un archivo .biom con taxonomías. BIOM es un formato "estándar" de archivos para almacenar datos de microbioma, tanto abundancias como taxonomías.

4.1 Convertir tabla (.qza) a biom y exportar taxonomía a .tsv

Exportar table.qza a al archivo feature-table.biom (ojo que el BIOM sólo tendrá la info de la abundancia) y luego lo transformamos el BIOM en un archivo de texto separado por tabulaciones (.tsv)

qiime tools export --input-path dada2out/table.qza --output-path exported

biom convert -i exported/feature-table.biom -o exported/feature-table.tsv --to-tsv

Exportar el archivo taxonomy.qza a un archivo de texto separado por tabulaciones (.tsv)

qiime tools export --input-path taxonomy.qza --output-path exported

Ambos archivos feature-table.tsv y taxonomy.tsv pueden ser revisados en la línea de comandos o directamente en Excel y ahí mezclarlos usando el ID como clave.

more exported/feature-table.tsv
more exported/taxonomy.qza

4.2 Añadir taxonomia a biom y extraer datos a tabla de texto (.tsv), OPCION1 (simple)

Como te comenté anteriormente, el archivo .BIOM puede almacenar taxonomía y eso es lo que haremos, agregar la taxonomía al archivo que exportamos desde table.qza. No es directo pero acá te muestro los pasos.

Modificar archivo taxonomy para que tenga el encabezado correcto

cd exported

sed -e 's/Taxon/taxonomy/' -e 's/^Feature ID/#OTUID/' taxonomy.tsv > taxonomy2.tsv

# La primera parte es para cambiar "Taxon" por "taxonomy"
# La segunda parte es para cambiar "Feature ID" por "#OTUID"

Luego, añadir la taxonomía al BIOM con este comando

biom add-metadata -i feature-table.biom -o feature-table-with-taxonomy.biom --observation-metadata-fp taxonomy2.tsv --sc-separated taxonomy --observation-header OTUID,taxonomy

De esta forma, el archivo feature-table-with-taxonomy.biom estará completo pero no se podrá importar en MEGAN pues, por alguna razón que desconozco, el formato biom que usa Qiime2 no es compatible con MEGAN.

4.3 Modificar BIOM de QIIME a BIOM compatible con MEGAN

Por alguna razón que desconozco, el formato biom que usa Qiime2 no es compatible con MEGAN

Lo que haremos será exportar el BIOM a texto y luego volver a importarlo BIOM. Es extraño pero sólo así funciona.

# Convertimos el BIOM en texto
biom convert -i feature-table-with-taxonomy.biom -o feature-table-with-taxonomy.tsv --header-key taxonomy --output-metadata-id "taxonomy" --to-tsv

El archivo feature-table-with-taxonomy.tsv es todo lo que necesitas para ver los datos en Excel.

Si la taxonomía que estás usando proviene de la base de datos Silva original tendrá esta forma D_0__*;D_1__*; D_2__*, etc
y tendrás que aplicar el siguiente comando para que el archivo BIOM final sea compatible con MEGAN. Si descargaste la base de datos Silva desde el sitio de Qiime2 o, si usas otra base de datos cuyo formato sea del tipo k__.*;p__*;c__.*, etc puedes saltar este paso.

## Comando extra sólo para el caso de taxonomía del tipo D_0__*;D_1__*; D_2__*
sed -i -e 's/ D_.__//g' -e 's/D_.__//g' feature-table-with-taxonomy.tsv

# En este comando busco cualquier aparición de D_0__ y la elimino (//). Lo mismo para D_1__, D_2__ etc.
# La gracia es que, como en la expresión que deseo eliminar sólo cambia el número del nivel taxonómico, puedo usar en reemplazo del número el comodín "." que significa "cualquier caracter" y así sed buscará cualquiera de los niveles y los eliminará.
# El reemplazo se realiza para dos casos, cuando comienza con un espacio ( D_.__) y cuando comienza sin espacio (D_.__).

Nos aseguramos que no exitan comillas en el archivo

sed -i 's/"//g' feature-table-with-taxonomy.tsv

Ahora transformamos el archivo de texto a BIOM

biom convert -i feature-table-with-taxonomy.tsv -o feature-table-with-taxonomy.MEGAN.biom --table-type="OTU table" --process-obs-metadata='taxonomy' --to-hdf5

El archivo feature-table-with-taxonomy.MEGAN.biom podrá ser importado en MEGAN

4.4 Extraer estadística básica desde el archivo BIOM.

# Estadística del conteo de reads por muestra
biom summarize-table -i feature-table-with-taxonomy.biom -o feature-table-with-taxonomy.summary.txt

# Estadística del conteo de OTUs (ASVs) por muestra
biom summarize-table -i feature-table-with-taxonomy.biom --qualitative -o feature-table-with-taxonomy.OTUs.summary.txt

4.5 Añadir taxonomía a BIOM, OPCION2

Estos pasos realizan lo mismo que la sección 3.2, pero de una forma alternativa

cd exported

Cambio el nombre de la primera columna del archivo taxonomy.tsv

sed -i "s/^Feature ID/OTU ID/" taxonomy.tsv

Cambio el nombre de la primera columna del archivo feature-table.tsv

sed -i "s/^#OTU ID/OTU ID/" feature-table.tsv

Elimino la primera línea del archivo feature-table.tsv

sed -i '1d' feature-table.tsv

Ordeno con awk los archivos pues el comando join (mas adelante) necesita que los IDs estén en orden.

awk 'NR == 1; NR > 1 {print $0 | "sort"}' taxonomy.tsv > taxonomy.sort.tsv
awk 'NR == 1; NR > 1 {print $0 | "sort"}' feature-table.tsv > feature-table.sort.tsv

Genero un archivo que mezcla feature-table.tsv con taxonomy.tsv usando el comando join

join feature-table.sort.tsv taxonomy.sort.tsv -t $'\t' > otu_table-filt.w_tax.txt

Elimino la última columna de otu_table-filt.w_tax.txt (confidence)

sed -i -r 's/\s+\S+$//' otu_table-filt.w_tax.txt # -r o -E (expresión regular extendida)
# comando alternativo: gawk -i inplace 'BEGIN{OFS="\t"} NF{NF--};1' otu_table-filt.w_tax.txt

sed -i 's/Taxon/taxonomy/' otu_table-filt.w_tax.txt

Nos aseguramos que no exitan comillas en el archivo

sed -i 's/"//g' feature-table-with-taxonomy.tsv

El archivo otu_table-filt.w_tax.txt es todo lo que necesitas para ver los datos en Excel.

En caso que quieras el BIOM, lee lo siguiente.

Finalmente, si la taxonomía proviene de la base de datos Silva (D_0__*;D_1__*; D_2__*, etc) tendrás que aplicar el siguiente comando para que el archivo BIOM final sea compatible con MEGAN, de lo econtrario te lo puedes saltar.

sed -i -e 's/ D_.__//g' -e 's/D_.__//g' otu_table-filt.w_tax.txt

# En este comando busco cualquier aparición de D_0__ y la elimino (//). Lo mismo para D_1__, D_2__ etc.
# La gracia es que, como en la expresión que deseo eliminar sólo cambia el número del nivel taxonómico, puedo usar en reemplazo del número el comodín "." que significa "cualquier caracter" y así sed buscará cualquiera de los niveles.
# El reemplazo se realiza para dos casos, cuando comienza con un espacio ( D_.__) y cuando comienza sin espacio (D_.__).

Nos aseguramos que no exitan comillas en el archivo

sed -i 's/"//g' feature-table-with-taxonomy.tsv

Con el siguiente comando podrás generar el archivo BOM con taxonomía para abrir en MEGAN.

biom convert -i otu_table-filt.w_tax.txt -o otu_table-filt.w_tax.biom --table-type="OTU table" --process-obs-metadata='taxonomy' --to-hdf5

Y bueno, con el archivo biom puedes aplicar los siguientes comandos para obtener la estadística de mapeo.

# Estadística del conteo de reads por muestra
biom summarize-table -i otu_table-filt.w_tax.biom -o otu_table-filt.w_tax.summary.txt

# Estadística del conteo de OTUs (ASVs) por muestra
biom summarize-table -i otu_table-filt.w_tax.biom --qualitative -o otu_table-filt.w_tax.OTUs.summary.txt