Qiime2-r16S: Guía de comandos Nanopore (Vsearch)

620  Jonathan Maldonado  
(last edited )

Tienes secuencias provenientes de una secuenciación masiva de amplicones r16S (r16S metabarcoding) utilizando Nanopore, 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 en formato fasta o fasta.gz

Vamos a utilizar como ejemplo un set de secuencias de ampliconres r16S secuenciados por Nanopore disponibles en NCBI. Se trata de una muestra marina recolectadad en el Océano Atlántico (PR) y en el Mar Mediterráneo (IL) en verano e invierno.

Este es el link al bioproject para obtener mas información: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA795711

Para facilitar el tutorial y usar menos recursos de cómputo sólo usaremos las muestras del Océano Atlántico, verano.

El set de datos es el siguiente con sus respectivos códigos en la base de datos SRA de NCBI:

  • 16S rRNA barcode sequences of Atlantic Ocean water in Winter:
    SRR17509858 SRR17509857 SRR17509856
  • 16S rRNA barcode sequences of Atlantic Ocean plastics in Winter:
    SRR17509855 SRR17509854 SRR17509853
  • 16S rRNA barcode sequences of Atlantic Ocean water in Summer:
    SRR17509852 SRR17509851 SRR17509849
  • 16S rRNA barcode sequences of Atlantic Ocean plastics in Summer:
    SRR17509848 SRR17509847 SRR17509846

Lo primero que haremos será descargar estas secuencias en nuestro disco duro utilizando un ambiente conda.

# Creamos el ambiente conda
conda create -n nanopore1 -c conda-forge -c bioconda -c cyclus parallel-fastq-dump porechop fastp kraken2 seqtk flye unicycler canu NanoPlot megahit NanoFilt java-jdk fastqc multiqc prokka pigz

# Activamos el ambiente conda
conda activate nanopore1

# Descargamos y descomprimimos la herramienta de NCBI que permite manejar archivos de secuencias
wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/sratoolkit.current-ubuntu64.tar.gz
tar xfz sratoolkit.current-ubuntu64.tar.gz

# Ejecutamos el comando para descargar los archivos
# El parámetro -pe indica el número de hilos de procesamiento a utilizar. Deben cambiar ese parámtro según corresponda.

sratoolkit.3.1.1-ubuntu64/bin/prefetch --max-size 1GB SRR17509856 && sratoolkit.3.1.1-ubuntu64/bin/fasterq-dump -pe 10 ./SRR17509856
sratoolkit.3.1.1-ubuntu64/bin/prefetch --max-size 1GB SRR17509857 && sratoolkit.3.1.1-ubuntu64/bin/fasterq-dump -pe 10 ./SRR17509857
sratoolkit.3.1.1-ubuntu64/bin/prefetch --max-size 1GB SRR17509858 && sratoolkit.3.1.1-ubuntu64/bin/fasterq-dump -pe 10 ./SRR17509858

sratoolkit.3.1.1-ubuntu64/bin/prefetch --max-size 1GB SRR17509853 && sratoolkit.3.1.1-ubuntu64/bin/fasterq-dump -pe 10 ./SRR17509853
sratoolkit.3.1.1-ubuntu64/bin/prefetch --max-size 1GB SRR17509854 && sratoolkit.3.1.1-ubuntu64/bin/fasterq-dump -pe 10 ./SRR17509854
sratoolkit.3.1.1-ubuntu64/bin/prefetch --max-size 1GB SRR17509855 && sratoolkit.3.1.1-ubuntu64/bin/fasterq-dump -pe 10 ./SRR17509855

## Estos son otros comandos que en teoría también permiten hacer la descarga pero en la práctica no me han funcionado.
# prefetch --max-size 1GB SRR17509856 && fasterq-dump -pe 4 ./SRR17509856
# prefetch --max-size 1GB SRR17509857 && fasterq-dump -pe 4 ./SRR17509857
# prefetch --max-size 1GB SRR17509858 && fasterq-dump -pe 4 ./SRR17509858
# prefetch --max-size 1GB SRR17509853 && fasterq-dump -pe 4 ./SRR17509853
# prefetch --max-size 1GB SRR17509854 && fasterq-dump -pe 4 ./SRR17509854
# prefetch --max-size 1GB SRR17509855 && fasterq-dump -pe 4 ./SRR17509855

# comprimimos todos los archivos para usar menos recursos de disco
pigz *.fastq

Ahora debemos remover las secuencias adaptadoras lo cual realizaremos en el ambiente nanopore1.

# Activamos el ambiente conda
conda activate nanopore1

# Utilizamos el programa porechop para remover adaptadores
# El parámetro -t indica el número de hilos de procesamiento a utilizar. Deben cambiar ese parámtro según corresponda.

porechop -i SRR17509856.fastq.gz --check_reads 1000 -t 4 -v 1 | gzip > SRR17509856-trim1.fastq.gz
porechop -i SRR17509857.fastq.gz --check_reads 1000 -t 4 -v 1 | gzip > SRR17509857-trim1.fastq.gz
porechop -i SRR17509858.fastq.gz --check_reads 1000 -t 4 -v 1 | gzip > SRR17509858-trim1.fastq.gz

porechop -i SRR17509853.fastq.gz --check_reads 1000 -t 4 -v 1 | gzip > SRR17509853-trim1.fastq.gz
porechop -i SRR17509854.fastq.gz --check_reads 1000 -t 4 -v 1 | gzip > SRR17509854-trim1.fastq.gz
porechop -i SRR17509855.fastq.gz --check_reads 1000 -t 4 -v 1 | gzip > SRR17509855-trim1.fastq.gz

Podemos revisar la calidad de las secuencias con fastqc y multiqc. Por ahora vamos a usar programa que está diseñado para secuencias Nanopore llamado "NanoPlot". Sólo lo haré para la primera secuencia como ejemplo.

# Activamos el ambiente conda
conda activate nanopore1

# El parámetro -t indica el número de hilos de procesamiento a utilizar. Deben cambiar ese parámtro según corresponda.

NanoPlot --fastq SRR17509856-trim1.fastq.gz -t 4 --plots kde -o Plot.all.SRR17509856-trim1

# Si el gráfico es muy extenso hacia la derecha (ie. existen unos pocos reads largos) podemos realizar una transformación logarítmica o elegir un tamaño máximo para mostrar:
#NanoPlot --fastq SRR17509856-trim1.fastq.gz -t 4 --plots kde --maxlength 50000 -o Plot50000.SRR17509856-trim1
#NanoPlot --fastq SRR17509856-trim1.fastq.gz -t 4 --plots kde --maxlength 10000 -o Plot10000.SRR17509856-trim1

El siguiente paso será revisar la calidad de las secuencias.

# Activamos el ambiente conda
conda activate nanopore1

# Utilizamos el programa porechop para remover adaptadores
# Nos quedaremos con secuencias que tengan al menos 1000 pb (-l 1000 y no mas de 2000 pb (--length_limit 2000 con una calidad promedio mayor o igual a 10 (-q 10)
# El parámetro -w indica el número de hilos de procesamiento a utilizar. Deben cambiar ese parámtro según corresponda.

fastp -i SRR17509856-trim1.fastq.gz --stdout -A -G -q 10 -l 1000 --length_limit 2000 -w 4 | gzip > SRR17509856-trim2.fastq.gz
fastp -i SRR17509857-trim1.fastq.gz --stdout -A -G -q 10 -l 1000 --length_limit 2000 -w 4 | gzip > SRR17509857-trim2.fastq.gz
fastp -i SRR17509858-trim1.fastq.gz --stdout -A -G -q 10 -l 1000 --length_limit 2000 -w 4 | gzip > SRR17509858-trim2.fastq.gz

fastp -i SRR17509853-trim1.fastq.gz --stdout -A -G -q 10 -l 1000 --length_limit 2000 -w 4 | gzip > SRR17509853-trim2.fastq.gz
fastp -i SRR17509854-trim1.fastq.gz --stdout -A -G -q 10 -l 1000 --length_limit 2000 -w 4 | gzip > SRR17509854-trim2.fastq.gz
fastp -i SRR17509855-trim1.fastq.gz --stdout -A -G -q 10 -l 1000 --length_limit 2000 -w 4 | gzip > SRR17509855-trim2.fastq.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 (un nombre que nos ayude a distinguir las muestras) y en las siguientes, los datos "extra", por ejemplo, el agrupamiento (cada dato separado por un tabulador)
Ejemplo:

ID Ocean Season Type SeasonType Replicate SRR
WW.r1 Atlantic Winter Water WW 1 SRR17509858
WW.r2 Atlantic Winter Water WW 2 SRR17509857
WW.r3 Atlantic Winter Water WW 3 SRR17509856
WP.r1 Atlantic Winter Plastic WP 1 SRR17509855
WP.r2 Atlantic Winter Plastic WP 2 SRR17509854
WP.r3 Atlantic Winter Plastic WP 3 SRR17509853

nano metadata.txt
Pegar las líneas y luego salir con control+X

Ejecutar
sed -i 's/ /\t/g' metadata.txt

nano labels.txt
Pegar las líneas y luego salir con control+X

Nuestras secuencias están preparadas. A partir de acá hay varias formas de proseguir y les mostraré la que siempre uso.

Genero un archivo de texto "labels.txt" que contiene el nombre que asigné a cada secuencia en el archivo metadata.txt

WW.r1
WW.r2
WW.r3
WP.r1
WP.r2
WP.r3

Nano

Comenzamos con la preparación de las muestras:

# Extraer/transformar secuencias en formato FASTA
# Para disminuir el tiempo de cómputo vamos a quedarnos con un máximo de 10.000 reads.
# Uso el comando "head" asumiendo que la secuencia de los archivos FASTA están en una sola línea.
# Si deseas trabajar con todas las secuencias debes eliminar esta parte del código "| head -n 20000"

seqtk seq -A SRR17509858-trim2.fastq.gz | head -n 20000 > SRR17509858-trim.fa
seqtk seq -A SRR17509857-trim2.fastq.gz | head -n 20000 > SRR17509857-trim.fa
seqtk seq -A SRR17509856-trim2.fastq.gz | head -n 20000 > SRR17509856-trim.fa

seqtk seq -A SRR17509855-trim2.fastq.gz | head -n 20000 > SRR17509855-trim.fa
seqtk seq -A SRR17509854-trim2.fastq.gz | head -n 20000 > SRR17509854-trim.fa
seqtk seq -A SRR17509853-trim2.fastq.gz | head -n 20000 > SRR17509853-trim.fa

# Renombrar los IDs de nuestras secuencias a un código que distinga cada set de datos.
# Desde luego, usaré el mismo nombre que definí en el archivo metadata.
sed -i 's/>.*/>WW.r1/' SRR17509858-trim.fa
sed -i 's/>.*/>WW.r2/' SRR17509857-trim.fa
sed -i 's/>.*/>WW.r3/' SRR17509856-trim.fa

sed -i 's/>.*/>WP.r1/' SRR17509855-trim.fa
sed -i 's/>.*/>WP.r2/' SRR17509854-trim.fa
sed -i 's/>.*/>WP.r3/' SRR17509853-trim.fa

# Concatenamos todos los archivos en un único archivo FASTA
cat *.fa > All-rename.fasta

En este momento todos los reads de nuestro archivo tienen el nombre de cada muestra, pero no se distinguen unos de otros. Vamos a agregar un conteo distintivo a cada read, lo cual nos permitirá "volver" al read si lo necesitamos.

# Descargamos un script que realiza el trabajo y lo ejecutamos.
wget https://github.com/jomaldon/scripts_bioinfo/raw/master/reads_add_number.pl
perl reads_add_number.pl All-rename.fasta labels.txt > All-rename2.fasta

Lo último es instalar un ambiente conda relacionado con Qiime2. Usaremos la versión 2024.2 debido a que es la última de las versionde de Qiime2 que permite descargar bases de datos pre-indexadas.
Más información en el siguiente link https://docs.qiime2.org/2024.2/install/native/

# Creamos el ambiente Qiime2
wget https://data.qiime2.org/distro/amplicon/qiime2-amplicon-2024.2-py38-linux-conda.yml
conda env create -n qiime2-amplicon-2024.2 --file qiime2-amplicon-2024.2-py38-linux-conda.yml

2. Manos a la obra

Todo lo que sigue es usando el ambiente Qiime2 por lo cual no olvides activarlo.

# Desactivamos el ambiente conda anterior
conda deactivate

# Activamos el ambiente qiime2
conda activate qiime2-amplicon-2024.2

2.1 Importar las secuencias

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

2.2 Dereplicar secuencias (agrupar 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
qiime tools export --input-path tableviz.qzv --output-path tableviz

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

# 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 nuestras secuencias provienen de una secuenciación Nanopore que posee un nivel de error mayor a que posee la secuenciación Illumina no podemos usar Dada2 (el cual requiere secuencias con alto nivel de confiabilidad).

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 closed-reference a un 90% de similitud con lo cual obtenemos ~2.100 secuencias representativas o OTUs. Lo habitual es usar un 97% o 99% de similitud, pero en este set de datos obtendríamos muchos OTUs lo cual hace mas intensivo el requerimiento de cómputo en los pasos siguientes.

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. Con el presente set de datos y usando 64 hilos debiera demorar ~7 minutos.

# Descargamos la base de datos Silva v138 en formato qiime2 (94 Mb)
wget https://data.qiime2.org/2024.2/common/silva-138-99-seqs.qza

# Ejecutamos la búsqueda con vsearch
qiime vsearch cluster-features-closed-reference \
  --p-threads 4 \
  --i-table table.qza \
  --i-sequences rep-seqs.qza \
  --i-reference-sequences silva-138-99-seqs.qza \
  --p-perc-identity 0.90 \
  --p-strand both \
  --o-clustered-table table-cr-90.qza \
  --o-clustered-sequences rep-seqs-cr-90.qza \
  --o-unmatched-sequences nomatch-seqs-cr-90.qza

# resumir table-cr-90.qza en qzv
qiime feature-table summarize --i-table table-cr-90.qza --o-visualization table-cr-90.qzv --m-sample-metadata-file metadata.txt
qiime tools export --input-path table-cr-90.qzv --output-path table-cr-90

# en la tabla podemos observar que la muestra WP.r2 tiene un total de 704 reads mapeados.
# luego le sigue la muestra WP.r1 con 1.669 reads y las siguientes tienen mas que eso.

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

2.4 Import the latest r16S database 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-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/2024.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.4.2 OPCION 2: 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 2a: Descargar archivos de la v138 listos para usar desde Qiime2

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

Opción 2b: 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.5 Classify the sequence variants (representative sequences)

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

El número de hilos a usar se controla con el comando --p-n-jobs. Yo usaré 64 hilos con lo cual tarda ~5 mins. Modificar según corresponda.

# Descargamos la base de datos Silva r16S v138 indexada para clasificación (508 Mb)
wget https://data.qiime2.org/2024.2/common/silva-138-99-nb-classifier.qza

qiime feature-classifier classify-sklearn \
  --i-classifier silva-138-99-nb-classifier.qza \
  --i-reads rep-seqs-cr-90.qza \
  --o-classification table-cr-90.tax.sklearn.qza \
  --p-n-jobs 4

Summarize the results

qiime metadata tabulate --m-input-file table-cr-90.tax.sklearn.qza --o-visualization table-cr-90.tax.sklearn.qzv
qiime tools export --input-path table-cr-90.tax.sklearn.qzv --output-path table-cr-90.tax.sklearn

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-cr-90.qza \
  --i-taxonomy table-cr-90.tax.sklearn.qza \
  --m-metadata-file metadata.txt \
  --o-visualization taxa-bar-plots-90-sklearn.qzv

qiime tools export --input-path taxa-bar-plots-90-sklearn.qzv --output-path taxa-bar-plots-90-sklearn

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). Toma ~15 minutos usando 64 hilos.
qiime phylogeny align-to-tree-mafft-fasttree \
        --p-n-threads 64 \
        --i-sequences rep-seqs-cr-90.qza \
        --o-alignment aligned-rep-seqs-cr-90.qza \
        --o-masked-alignment masked-aligned-rep-seqs-cr-90.qza \
        --o-tree unrooted-tree.qza \
        --o-rooted-tree rooted-tree.qza

# Exportamos el archivo 
qiime tools export --input-path rooted-tree.qza --output-path rooted-tree
# En varios pasos (alternativa)
qiime alignment mafft \
  --i-sequences rep-seqs-cr-90.qza \
  --p-n-threads 64 \
  --o-alignment aligned-rep-seqs-cr-90.qza

qiime alignment mask \
  --i-alignment aligned-rep-seqs-cr-90.qza \
  --o-masked-alignment masked-aligned-rep-seqs-cr-90.qza

qiime phylogeny fasttree \
  --i-alignment masked-aligned-rep-seqs-cr-90.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 WP.r2 con 704 la cual descartaré por tener pocos reads comparada con las otras muestras. La que sigue es WP.r1 con 1669 reads, entoces mi punto de corte de rarefacción será 1600 reads.

# Este comando toma ~7 minutos usando 64 hilos de procesamiento.

qiime diversity alpha-rarefaction \
  --i-phylogeny rooted-tree.qza \
  --i-table table-cr-90.qza \
  --p-max-depth 1600 \
  --p-metrics observed_features chao1 shannon faith_pd goods_coverage \
  --m-metadata-file metadata.txt \
  --o-visualization rarefaction_r1600

# 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

# Este comando toma ~2 minutos usando 64 hilos de procesamiento.

qiime diversity core-metrics-phylogenetic \
  --i-phylogeny rooted-tree.qza \
  --i-table table-cr-90.qza \
  --p-sampling-depth 1600 \
  --p-n-jobs-or-threads 64 \
  --m-metadata-file metadata.txt \
  --output-dir core-metrics-phylogenetic-metrics_r1600

3.4 Obtain the rarefaction table (tsv and biom format)

Vamos a extraer la información de contero de datos rarefaccionda y la taxonomía en una misma tabla que podremos abrir en Excel.

Lo primero es exportar el archivo que contiene las taxonomias (table-cr-90.tax.sklearn.qza) a un archivo de texto separado por tabulaciones (.tsv). Dicho archivo se llamará taxonomy2.tsv.

qiime tools export --input-path table-cr-90.tax.sklearn.qza --output-path exported_r1600
sed -e 's/Taxon/taxonomy/' -e 's/^Feature ID/#OTUID/' exported_r1600/taxonomy.tsv > exported_r1600/taxonomy2.tsv

Ahora podemos exportar la tabla de conteo de datos y mezclarla con la tabla de taxonomía.


cd core-metrics-phylogenetic-metrics_r1600

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

cd ../exported_r1600

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

# 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
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 Extraer estadística básica desde el archivo BIOM

Con estos comandos obtendremos dos archivos de texto con el resumen de datos de nuestras muestras rarefaccionadas.

# 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

3.6 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_r1600/rarefied_table.qza --i-taxonomy table-cr-90.tax.sklearn.qza --m-metadata-file metadata.txt --o-visualization taxa-bar-plots_r1600.qzv

qiime tools export --input-path taxa-bar-plots_r1600.qzv --output-path taxa-bar-plots_r1600

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

Este sitio usa Akismet para reducir el spam. Aprende cómo se procesan los datos de tus comentarios.

es Español
en Englishfr Françaisde Deutsches Españolpt Portuguêssk Slovenčinaru Русскийzh-CN 简体中文