SITUACIÓN DONDE ESTA GUÍA ES ÚTIL: Tienes reads forward y reverse limpios (sin adaptadores) 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 Dada2 en Qiime2.

Si tienes reads forward y reverse limpios (sin adaptadores) de una secuenciación masiva de amplicones de ITS este protocolo también te sirve, sólo debes cambiar las bases de datos. Por otro lado, si quieres usar el plugin Q2-ITSxpress, revisa esta otra guía. En lo personal prefiero no usar el plugin Q2-ITSxpress pues no he encontrado mayores diferencias. En resumen, la siguiente guía te sirve tanto para r16S como para ITS u otro marcador.

Si tus reads ya están unidos (típico del post-proceso de MrDNA o BGI) y quieres realizar en Qiime2 un análisis similar a Qiime1 puedes adaptar para ITS el protocolo VSearch que detallo en esta 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 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 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
C-400-r1,$PWD/Sample2_1.fq.gz,forward
C-400-r1,$PWD/Sample2_2.fq.gz,reverse
...

Si sólo posees secuencias forward o secuencias forward y reverse que ya han sido unidas [ej. MrDNA o BGI], 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 Identify sequence variants (representative sequences)

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 indica cuál será el tamaño máximo que permitiremos para los reads. Como la calidad de los reads decae hacia el extremo 3', lo recomendado es revisar la calidad promedio de los reads y eliminar las sección a partir de la cual la calidad decae. Esto se hace revisando el archivo sequences.qzv, sección "Interactive Quality Plot", seleccionando desde el gráfico el número de base desde donde comienza la mala calidad. Ese valor seleccionado será el p-trun-len. El valor seleccionado no debe ser mayor al tamaño de los reads. El tamaño promedio de los reads se puede revisar más abajo del gráfico en "Demultiplexed sequence length summary". 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 de reads para usar la info de secuencia más larga)
qiime dada2 denoise-paired --i-demultiplexed-seqs sequences.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-paired --p-n-threads 64

ln -s dada2out-paired dada2out

# Método single (usa las secuencias sin aparearlas, es decir, cortas o el caso que ya estén "unidas" [joined])
qiime dada2 denoise-single --i-demultiplexed-seqs sequences.qza --p-trunc-len 220 --p-max-ee 2  --p-trunc-q 2 --output-dir dada2out.single --p-n-threads 64

ln -s dada2out-single dada2out

# 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)
qiime feature-table tabulate-seqs --i-data dada2out/representative_sequences.qza --o-visualization rep-seqs.qza.qzv

2.3 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]).

Como referencia, este es el sitio web donde podrán encontrar la información mas actualizada de los archivos que se usan en los siguientes pasos: https://docs.qiime2.org/2024.2/data-resources/

2.3.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 pero ojo que el archivo incluye 16S y 18S.

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

# ejecutamos la clasificación
qiime feature-classifier classify-consensus-vsearch \
  --i-query dada2out/representative_sequences.qza \
  --i-reference-reads silva-138-99-seqs.qza \
  --i-reference-taxonomy silva-138-99-tax.qza \
  --p-strand both \
  --p-threads 96 \
  --p-perc-identity 0.97 \
  --o-classification taxonomy.vsearch97.qza \
  --output-dir vsearch97

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
awk '/^>/ {print($0)}; /^[^>]/ {print(toupper($0))}' 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
qiime feature-classifier classify-consensus-vsearch \
  --i-query dada2out/representative_sequences.qza \
  --i-reference-reads silva_132_99_16S.qza \
  --i-reference-taxonomy 16S99.taxonomy_7_levels.qza \
  --p-strand both \
  --p-threads 96 \
  --p-perc-identity 0.97 \
  --o-classification taxonomy.vsearch97.qza \
  --output-dir vsearch97

# También podemos usarlos para generar el archivo necesario para el método classify-sklearn donde necesitamos un archivo classify (ver 2.3.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-132-99-16S-nb-classifier.qza

2.3.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 o Silva v132), 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 si la base de datos tiene mas de 100mil secuencias
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 Classify the sequence variants (representative sequences)

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

Ojo con el número de trabajos a lanzar... yo he usado la opción 0 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.

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

Summarize the results

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

2.5 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 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 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.

2.6 Create data subset with the abundance average according to a sample group

Es posible agrupar las muestras según algún la propiedad que se desee. Por ejemplo, si el archivo metadata tiene la columna "Altitude" que incluye las categorías "High, Medium, Low", es posible agrupar las muestras según dichas categorías. De esta forma, el nuevo archivo ya no tendrá "N número de muestras", sinó que 3 categorías donde se han promediado las abundancias de cada muestra perteneciente a cada categoría. Eso permite hacer estadística por categoría.

qiime feature-table group --i-table dada2out/table.qza --p-axis sample --m-metadata-file metadata.txt --m-metadata-column "Altitude" --p-mode mean-ceiling --o-grouped-table dada2out/table.Altitude.qza

Ahora podemos crear un gráfico interactivo para visualizar las abundancias de cada taxa, en cada categoría mediante barras apiladas.

qiime taxa barplot --i-table dada2out/table.Altitude.qza --i-taxonomy taxonomy.qza --o-visualization taxa-bar-plots.Altitude.qzv

3.1 Alineamiento múltiple de secuencias representativas

3. Análisis de Diversidad

# 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

qiime tools export --input-path rarefaction_r7k.qzv --output-path 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

nohup 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.MEGAN.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_r7k/rarefied_table.qza --i-taxonomy taxonomy.qza --m-metadata-file metadata.txt --o-visualization exported_r7k/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
wc -l exported/*.tsv

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 [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