Tienes secuencias unidas (joined or extended, MrDNA, BGI) provenientes de reads forward y reverse de una secuenciación masiva de amplicones r16S (r16S metabarcoding), esta guía te servirá para realizar un análisis de taxonomía usando el protocolo Vsearch en Qiime2. Esta guía también te sirve para secuencias ITS pero cambiando la base de datos de referencia.

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

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

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.

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 unidas (reads forward y reverse alineados) en formato fasta o fasta.gz

1.2 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

qiime tools import 
  --input-path All-rename.fasta 
  --output-path sequences.qza 
  --type 'SampleData[Sequences]'

2.2 Dereplicar secuencias (eliminar secuencias idénticas)

qiime vsearch dereplicate-sequences --i-sequences sequences.qza --o-dereplicated-table table.qza --o-dereplicated-sequences rep-seqs.qza

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

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

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

2.3 Identify sequence variants with Vsearch (representative sequences)

Dado que las secuencias ya vienen alineadas, filtradas por calidad y de-replicadas, no podemos usar Dada2 (el cual requiere secuencias "crudas"). La forma de proceder es con vsearch lo cual es un símil a lo que se realiza en Qiime1 con usearch.

Las estrategias a utilizar pueden ser 3 (más info en este link): 

1. closed-reference: Las secuencias son mapeadas a la base de datos directamente. Es rápido pero la desventaja es que no se encuentran OTUs "nuevos" (sólo identificas lo que ya se ha identificado en las bases de datos).
2. open-reference: Las secuencias son mapeadas a la base de datos generando un primer set de OTUs. Las sequencias que no mapearon son utilizadas para encontrar OTUs por la técnica de novo (ver punto 3) y finalmente se clusterizan todos los OTUs.
3. de novo: Las secuencias son analizadas y clusterizadas para obtener OTUs representantes sin usar una base de datos. Requiere mucho poder de cómputo, tiempo… y paciencia.

La recomendación general es usar open-reference.

Vamos a realizar la búsqueda de OTUs con vsearch usando la estrategia open-reference a un 97% de similitud.

El comando necesitará definir el número de hilos de procesador a usar (p-threads). Cambiar el número según la capacidad del computador.

Según el volumen de datos, este paso puede demorar horas.

qiime vsearch cluster-features-open-reference 
  --p-threads 50 
  --i-table table.qza 
  --i-sequences rep-seqs.qza 
  --i-reference-sequences silva-138-99-seqs.qza 
  --p-perc-identity 0.97 
  --p-strand both 
  --o-clustered-table table-or-97.qza 
  --o-clustered-sequences rep-seqs-or-97.qza 
  --o-new-reference-sequences new-refs-seqs-or-97.qza

# resumir table-or-97.qza en qzv
qiime feature-table summarize --i-table table-or-97.qza --o-visualization table-or-97.qzv --m-sample-metadata-file metadata.txt

# Tabla con información de las secuencias representativas (no hacer)
tqiime feature-table --i-data rep-seqs-or-97.qza --o-visualization rep-seqs-or-97.qzv

2.4 Import the latest r16S data into QIIME2

La clasificación taxonómica usando base de datos de referencia (en este caso la base de datos ribosomal Silva), se puede realizar con distintas estrategias/métodos:
El método classify-consensus-vsearch necesitará un archivo con las secuencias de referencia y otro con sus taxonomía separados. Según el número de secuencias a alinear puede demorar muuuucho y consumir muchos recursos (OPCION 1 [no recomendada]).
El método classify-sklearn necesitará un "clasificador" de taxonomía, lo cual es un archivo preparado (entrenado) a partir de un archivo de taxonomía y de secuencias de referencia. Requiere menos recursos que el método anterior (OPCION 2 [recomendada]).

2.4.1 OPCION 1: Archivos para usar método classify-consensus-vsearch [no recomendado]

La base de datos que utilizaremos será Silva v138. En el sitio web oficial de Silva no está disponible la versión 138 para qiime, sólo podremos usar la v132 y tendremos que "manipular" un poco los archivos para importarlos a Qiime2. En cambio, en los recursos proporcionados en el sitio oficial de Qiime2, podemos encontrar los archivos de la v138 listos, llegar y usar

Opción 1a: Descargar archivos de la v138 listos para usar desde Qiime2

# secuencias
wget https://data.qiime2.org/2021.2/common/silva-138-99-seqs.qza
# taxonomías
wget https://data.qiime2.org/2021.2/common/silva-138-99-tax.qza

Opción 1b: Preparando los archivos desde datos crudos (sólo Silva v132)

# Descargamos y descomprimimos las secuencias
wget https://www.arb-silva.de/fileadmin/silva_databases/qiime/Silva_132_release.zip
unzip Silva_132_release.zip

# Transformamos minúsculas en mayúsculas
sed -e 's/_//g' -e 's/a/A/g' -e 's/c/C/g' -e 's/t/T/g' -e 's/g/G/g' -e 's/n/N/g' SILVA_132_QIIME_release/rep_set/rep_set_16S_only/99/silva_132_99_16S.fna > silva_132_99_16S-edited.fna

# Importamos las secuencias a Qiime2
qiime tools import 
  --type 'FeatureData[Sequence]' 
  --input-path silva_132_99_16S-edited.fna 
  --output-path silva_132_99_16S.qza
  
# Importamos la taxonomía a Qiime2
qiime tools import --type 'FeatureData[Taxonomy]' 
  --input-format HeaderlessTSVTaxonomyFormat 
  --input-path SILVA_132_QIIME_release/taxonomy/16S_only/99/taxonomy_7_levels.txt
  --output-path 16S99.taxonomy_7_levels.qza

# Con estos archivos podemos usar el método classify-consensus-vsearch

# También podemos usarlos para generar el archivo necesario para el método classify-sklearn donde necesitamos un archivo classify (ver 2.4.2). Este paso puede demorar varias horas.

qiime feature-classifier fit-classifier-naive-bayes --i-reference-reads silva_132_99_16S.qza --i-reference-taxonomy 16S99.taxonomy_7_levels.qza --o-classifier silva-138-99-nb-classifier.qza

2.4.2 OPCION 2: Archivos para usar método classify-sklearn [recomendado]

El método classify-sklearn necesitará un "clasificador" de taxonomía, lo cual es un archivo preparado (entrenado) a partir de un archivo de taxonomía y de secuencias de referencia. Existen dos formas de proceder en este punto: Generar nuestro propio archivo "clasificador" a partir de los archivos crudos descargados desde el sitio web oficial o descargar el "clasificador" listo, el cual ha sido preparado por el equipo de Qiime2.

Claramente usaremos el archivo listo para usar. Más info en los links que dejé en cada opción.

wget https://data.qiime2.org/2021.2/common/silva-138-99-nb-classifier.qza

Como alternativa, en caso que quieras entrenar tu propio clasificador (por ejemplo si usas otra base de datos como ITS-UNITE), necesitarás tener el qza de las secuencias de referencia y el de la taxonomía (ver Opción 1). Con esos archivos 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.

# Este paso puede demorar varias horas
qiime feature-classifier fit-classifier-naive-bayes --i-reference-reads unite82s.dyn-seqs.qza --i-reference-taxonomy silva-138-99-tax.qza --o-classifier silva-138-99-nb-classifier.qza

2.5 Classify the sequence variants (representative sequences)

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

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 silva-138-99-nb-classifier.qza 
  --i-reads rep-seqs-or-97.qza 
  --o-classification table-or-97.tax.sklearn.qza 
  --p-n-jobs -1

Summarize the results

qiime metadata tabulate --m-input-file table-or-97.tax.sklearn.qza --o-visualization table-or-97.tax.sklearn.qzv

2.6 Create an interactive bar plot figure

El siguiente comando permite crear un gráfico interactivo para visualizar las abundancias de cada taxa, en cada muestra mediante barras apiladas.

qiime taxa barplot 
  --i-table table-or-97.qza 
  --i-taxonomy table-or-97.tax.sklearn.qza 
  --m-metadata-file metadata.txt 
  --o-visualization taxa-bar-plots-97-sklearn.qzv

A partir de la tabla de abundancias (table.qza) y la tabla de taxonomías (table-or-97.tax.sklearn.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 diversidad alfa y 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 rep-seqs-or-97.qza 
        --o-alignment aligned-rep-seqs-or-97.qza 
        --o-masked-alignment masked-aligned-rep-seqs-or-97.qza 
        --o-tree unrooted-tree.qza 
        --o-rooted-tree rooted-tree.qza

# En varios pasos (alternativa)
qiime alignment mafft 
  --i-sequences rep-seqs-or-97.qza 
  --p-n-threads 64 
  --o-alignment aligned-rep-seqs-or-97.qza

qiime alignment mask 
  --i-alignment aligned-rep-seqs-or-97.qza 
  --o-masked-alignment masked-aligned-rep-seqs-or-97.qza

qiime phylogeny fasttree 
  --i-alignment masked-aligned-rep-seqs-or-97.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 table-or-97.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 table-or-97.tax.sklearn.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 de QIIME 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 [no recomendada]

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' otu_table-filt.w_tax.txt

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' otu_table-filt.w_tax.txt

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