PASO1. OBTAINING THE DATA. Descargar el SAMPLE METADATA, lecturas paired-end y BARCODES
En este tutorial aprenderás a procesar archivos de secuencias paired-end que son diferentes de las lecturas single-end que se usan en el Tutorial de Moving Pictures.
Por lo que éste tutorial incluye los pasos iniciales de procesamiento de lecturas paired-end, hasta el punto en que los pasos de análisis son idénticos a los realizados para las lecturas single-end del tutorial de moving pictures.
Este tutorial incluye los pasos de importación de las secuencias a formato QZA, el desmultiplexing y el denoising de las lecturas, para obtener como resultado la tabla de frecuencias de las features, las secuencias representativas y las estadísticas de resultados del denoising, donde podrás ver las secuencias que fueron eliminadas en cada paso y las ASVs inferidas por cada muestra.
Los datos de este tutorial provienen de muestras de suelo del desierto de Atacama, en el norte de Chile, el cuál es uno de los lugares más áridos de la Tierra, con algunas áreas que reciben menos de un milímetro de lluvia por década. A pesar de esta extrema aridez, hay microorganismos viviendo en el suelo.
Para este estudio, se tomaron muestras de suelo siguiendo dos transectos de este a oeste llamados Baquedano y Yungay, donde se cavaron núcleos y se recogieron muestras a tres profundidades.
El primer paso de este tutorial es construir una carpeta con el comando de linux mkdir para almacenar todos los archivos que se generen durante el análisis, con el comando cd ingresamos a ésta carpeta (qiime2-atacama-tutorial).
Posteriormente vamos a descargar el documento de los metadatos, damos click en el enlace y una vez que el archivo se ha descargado en la carpeta de Downloads, lo movemos con el comando de linux mv a la carpeta que recién creamos
Abrimos el archivo de metadatos en excel para ver la información que contiene:
En la primer columna vemos el nombre de cada muestra, la segunda columna contiene la secuencia de los barcodes utilizados, la tercer columna la elevación de los sitios de estudio, la cuarta incluye la concentración de DNA, la quinta la concentración de amplicon del 16S rRNA, la séptima el nombre del transecto, la octava el nombre del sitio, la novena la profundidad a la que fue tomada la muestra y el resto de las columnas contienen información de las características fisicoquímicas del suelo.
Lo siguiente es descargar las secuencias multiplexadas y los barcodes, por lo que bajaremos 3 archivos .gz a la computadora, uno que corresponde a las secuencias forward, otro a las reverse y otro a los barcodes.
Para ello, primero vamos a crear una carpeta llamada emp-paired-end-sequences donde almacenaremos estos archivos.
Luego, damos click en cada uno de los enlaces y una vez que ha concluido la descarga, moveremos los archivos con el comando mv a la carpeta recién creada, aprovecharemos para utilizar un comodín * y moverlos en un solo paso.
Si hacemos un listado con el comando ls, veremos que ya tenemos los archivos en la carpeta correspondiente.
Para ver los pasos mencionados puedes revisar el siguiente video.
En microbioma lab ofrecemos servicios de análisis bioinformático de datos de éste tipo de lecturas y de metagenómica, puedes ver el catálogo de productos y servicios así como los costos.
Recuerda seguirnos en nuestras redes sociales: Facebook e Instagram donde compartimos videos de comandos de linux, tips de Rstudio, información de científicos que estudian microbioma, conceptos básicos y revistas donde puedes publicar tu investigación de microbioma.
Cursos de bioinformática disponibles
PASO 2. Importar PAIRED-END READS a QIIME2 en artefacto QZA
Ya que tenemos descargadas las secuencias, los barcodes y el archivo de metadatos hay que importar las secuencias y los barcodes a Qiime2 transformándolos en artefacto QZA.
Debemos estar ubicados en el interior de la carpeta qiime2-atacama-tutorial que creamos anteriormente, si no estás en ella te puedes mover con el comando de linux cd. Activa Qiime 2 usando el siguiente comando:
conda activate qiime2-2022.11
Verifica la versión que tienes instalada de Qiime 2 en tu computadora porque eso puede hacer que el comando cambie ligeramente al que te mostramos aquí.
Lo siguiente es usar el plugin qiime tools import. Puedes copiar la caja de comando dando click aquí y lo pegas en la terminal.
La primera línea indica el plugin tool y el método para importar, después, varias opciones, la primera es para el tipo de secuencias en este caso paired-end del Proyecto de Microbioma de la Tierra (indicado por las siglas EMP, Earth Microbiome Project), seguido de la opción de entrada, en este caso la carpeta que contiene las secuencias y barcodes que descargamos previamente, seguido de la opción de salida que es el nombre del directorio en que se generará el archivo de las secuencias importadas como un artefacto qza. Damos enter para ejecutar. ¡Listo! Ya tenemos las secuencias listas para continuar con el siguiente paso que es el demultuplexing.
Para ver detalladamente estos pasos te recomendamos ver el siguiente video:
PASO 3. DEMULTIPLEXING secuencias PAIRED-END en QIIME2
Ya que importamos las secuencias y los barcodes a Qiime2 como un artefacto QZA. Lo siguiente es hacer el demultiplexing de las secuencias, es decir asignar el barcode a cada muestra con ayuda del archivo de metadatos.
Para comenzar, copiamos el comando y pégalo en al terminal.
qiime demux emp-paired \
–m-barcodes-file sample-metadata.tsv \
–m-barcodes-column barcode-sequence \
–p-rev-comp-mapping-barcodes \
–i-seqs emp-paired-end-sequences.qza \
–o-per-sample-sequences demux-full.qza \
–o-error-correction-details demux-details.qza
Éste es el plugin demux y el método es para el demultiplexing de secuencias paired-end, nota que es distinto al de secuencias single-end.
La primera opción es para la ruta del archivo de metadatos que descargamos previamente, ya que contiene las secuencias de los barcodes utilizados.
La segunda opción, es para indicar al plugin el nombre de la columna en el archivo de metadatos que contiene la secuencia de los barcodes de cada muestra.
La tercera opción es un parámetro para indicar que las secuencias de los barcodes en el archivo de metadatos son reverso complementarias a las que se encuentran en los archivos de secuencias, esto aplica para éste conjunto de datos del suelo de Atacama, si tienes dudas si ésta opción aplica para tus datos envíanos un correo a contacto@microbioma-lab para agendar una consultoría.
La siguiente opción se refiere a la ruta del artefacto .qza que se importó anteriormente, la cual contiene los archivos de barcodes y secuencias.
La siguiente opción es para indicar la ruta del directorio para guardar el artefacto de salida, el cual contendrá las secuencias demultiplexadas. De acuerdo a este tutorial, al no indicar el directorio, se guardarán todos los archivos en la carpeta en la que estamos ubicados. Lo mismo aplica para los archivos de entrada, al no indicar una ruta, significa que los archivos están ubicados en la carpeta en la que estamos posicionados.
Y la última opción es para indicar la ruta y/o el nombre del archivo donde se almacenan los detalles de correcciones de los barcodes.
Una vez que se han generado los artefactos, el siguiente paso para este tutorial es hacer un submuestreo de las secuencias, es decir, seleccionar aleatoriamente un conjunto más pequeño de secuencias. Es importante recalcar que esto se hace para este tutorial por dos motivos, uno para acelerar el tiempo de ejecución de éste tutorial reduciendo la cantidad de información y el otro para mostrarte su funcionalidad. Si estás pensando en hacerlo con tus muestras, asegúrate de tener una justificación razonable para ello. El comando es el siguiente:
qiime demux subsample-paired \
–i-sequences demux-full.qza \
–p-fraction 0.3 \
–o-subsampled-sequences demux-subsample.qza
El primer script es el método subsample-paired del plugin demux para hacer el submuestreo de las secuencias paired-end, la primera opción es para el archivo de entrada, en este caso el artefacto que se generó del demultiplexing del paso anterior.
La siguiente opción es para indicar la fracción de secuencias del submuestreo. En este caso el 30% de las secuencias. Y la última opción es para indicar el nombre del artefacto de salida.
El siguiente script es para hacer un resumen de los recuentos por muestra, generando gráficas interactivas basadas en “n” secuencias seleccionadas al azar, donde:
qiime demux summarize \
–i-data demux-subsample.qza \
–o-visualization demux-subsample.qzv
La primera opción es para el archivo de entrada que es el artefacto del submuestreo que se acaba de generar y la siguiente opción es para indicar el nombre del archivo de visualización .qzv.
El archivo que se genera muestra un resumen de las secuencias demultiplexadas del submuestreo, vemos el total de secuencias tanto del forward como del reverse, el número máximo y el mínimo de secuencias en las muestras. También, hay un histograma del número de secuencias con respecto al número de muestras, y a diferencia de las secuencias single end, aquí se muestran tanto las del forward como del reverse, ya que estamos trabajando con secuencias paired end. El histograma nos muestra que aproximadamente 20 muestras tuvieron menos de 1000 secuencias, mientras que menos de 5 muestras tuvieron más de 4000 secuencias. Después vemos las secuencias por muestra, ordenadas de mayor a menor. Estos datos también los puedes descargar en formato .tsv dando clic en el enlace.
En la pestaña interactive quality, vemos una gráfica de la calidad de las secuencias tanto para lecturas forward como reverse, en el eje de las abscisas está la posición de la secuencia en pares de bases mientras que en las ordenadas el Quality score. Observamos que la calidad de las secuencias es distinta para las lecturas forward respecto de las reverse, por ejemplo si hacemos un zoom al comienzo de la gráfica vemos que para el caso de las lecturas forward el valor de Quality score se encuentra entre 30 y 34 mientras que en las secuencias reverse está entre 16 y 32. Así, si hacemos un zoom al nucleótido 5 vemos que la mediana del quality score en las secuencias forward es de 32 mientras que en las reverse es de 29. De manera general podemos ver que el quality score de las secuencias reverse es de menor calidad que de las forward
En la siguiente tabla vemos el resultado del largo de las lecturas en número de nucleótidos obtenidos por cada 10,000 secuencias muestreadas tanto del forward como del reverse.
Si observamos nuevamente la tabla de secuencias por muestra podemos notar que hay algunas que tienen menos de 100 lecturas. En este tutorial se realiza un filtrado de muestras para eliminar aquellas con menos de 100 lecturas. Toma en cuenta que este paso en el tutorial pretende demostrar cómo se hace un filtrado, si lo quieres hacer con tus muestras asegúrate de tener una razón justificable para hacerlo. El script es el siguiente:
qiime tools export \
–input-path demux-subsample.qzv \
–output-path ./demux-subsample/
El primer script es del plugin tools y el método export es para exportar la información del archivo de visualización qzv. La primera opción es para el archivo de entrada, en este caso el de visualización que acabamos de ver y la siguiente opción es para el archivo de salida que será la carpeta que contendrá varios documentos, tablas e imágenes.
El siguiente plugin es para realizar el filtrado de las muestras. La primera opción es para el archivo de entrada en este caso el artefacto del demultiplexing del submuestreo. Posteriormente el archivo de metadatos que contiene la información necesaria para hacer el filtrado, en este caso, es el número de secuencias por muestra, información que está contenida en uno de los archivos que se generó con el plugin anterior: per-sample-fastq-counts.tsv.
qiime demux filter-samples \
–i-demux demux-subsample.qza \
–m-metadata-file ./demux-subsample/per-sample-fastq-counts.tsv \
–p-where ‘CAST([forward sequence count] AS INT) > 100’ \
–o-filtered-demux demux.qza
La siguiente opción es para especificar qué muestras serán incluídas en el filtrado final, en este caso se especifica dentro de los corchetes el nombre de la columna que contiene la información para hacer el filtrado, es decir, la columna forward sequence count, AS INT les indica que son números enteros y > 100 que retenga las muestras con un número de lecturas mayor a 100. Finalmente, la opción para el nombre y ruta del archivo de salida, un artefacto qza que contendrá las muestras con más de 100 lecturas.
En el siguiente video puedes ver los pasos que acabamos de explicar.
Listo, ya tenemos las secuencias listas para continuar con el siguiente paso que es el denoising e inferencias de ASVs.
Paso 4. DENOISING e inferencia de ASV de secuencias PAIRED-END en QIIME2
Ya se hizo el demultiplexing de las secuencias y revisamos las gráficas de calidad, el siguiente paso es hacer el denoising y la inferencia de ASVs, para ello vamos a usar el siguiente código:
qiime dada2 denoise-paired \
–i-demultiplexed-seqs demux.qza \
–p-trim-left-f 13 \
–p-trim-left-r 13 \
–p-trunc-len-f 150 \
–p-trunc-len-r 150 \
–o-table table.qza \
–o-representative-sequences rep-seqs.qza \
–o-denoising-stats denoising-stats.qza
Éste es el plugin que se utiliza para hacer el denoising de las secuencias paired-end, nota que es diferente de aquel que se utiliza para secuencias single end.
La primera opción es para el archivo de entrada que son las secuencias demultiplexadas que obtuvimos en el paso anterior, las siguientes dos opciones son para hacer un corte al inicio de las secuencias tanto forward como reverse, éste se realiza en el nucléotido 13 con base en lo observado en el plot de calidad que vimos después de hacer el demultiplexing, y se hace considerando que al inicio de ambas secuencias hay menor calidad respecto al resto de la lectura.
Las siguientes dos opciones corresponden a la posición en las secuencias en donde se realizará el truncado, en ese caso en 150, tanto para el forward como para el reverse.
Las últimas tres opciones son para generar los archivos de salida como artefactos qza, los cuales corresponden a la tabla de frecuencias de las ASV, las secuencias representativas y los resultados de las estadísticas del denoising realizadas por dada2.
Una vez que se han generado los archivos de salida hay que obtener los archivos de visualización en formato .qzv para eso se utiliza el siguiente código:
qiime feature-table summarize \
–i-table table.qza \
–o-visualization table.qzv \
–m-sample-metadata-file sample-metadata.tsv
El código anterior es para obtener la tabla de frecuencias, para ello se necesita la tabla qza que se obtuvo con el código anterior y el documento de los metadatos con extensión .tsv, para devolvernos un archivo de visualización con extensión .qzv
qiime feature-table tabulate-seqs \
–i-data rep-seqs.qza \
–o-visualization rep-seqs.qzv
El código anterior es para obtener las secuencias representativas en un archivo de visualización qzv, para ello se necesitan las secuencias representativas generadas con el plugin de dada2.
qiime metadata tabulate \
–m-input-file denoising-stats.qza \
–o-visualization denoising-stats.qzv
Finalmente, el código anterior se utiliza para obtener el resumen de los estadísticos que generó dada2 durante el denoising, tomando como archivo de entrada el artefacto qza de las estadísticas del denoising y obteniendo un archivo de visualización qzv.
Podemos visualizar los archivos generados en nuestra computadora o abrir los del sitio web, ya que corresponden a los mismos.
Cuando abrimos el archivo de la tabla de frecuencias observamos: una tabla resumen, donde vemos que de un total del 54 muestras se obtuvieron 1,138 features con una frecuencia total de 65,879.
También vemos una tabla con la frecuencia mínima, la mediana y la máxima para las features por muestra. Hay una gráfica de histograma donde vemos que, por ejemplo hay muestras que tuvieron un mínimo de 6 y otras un máximo de 2,152 features. Luego vemos una tabla de la frecuencia de cada feature, por ejemplo una de ellas tuvo una frecuencia mínima de 2, mientras que otra tuvo una frecuencia máxima de 1,438.
En la siguiente pestaña podemos ver un mapa interactivo donde están graficados los resultados con base en las categorías del archivo de los metadatos, como la secuencia del barcode, pH, humedad relativa, carbono orgánico total, etc. Por ejemplo, si escogemos pH, podemos ver en la gráfica el número de muestras relacionadas con cada medida de pH, lo mismo para conductividad electrolítica o humedad relativa, entre otros.
También podemos mover la profundidad de muestreo a un número específico de observaciones, por ejemplo 295, y vemos en color gris las muestras que se eliminarían del análisis si retuviéramos únicamente aquellas que contengan al menos 295 observaciones. Esto lo podemos ver con más detalle en la tabla que está debajo, donde hay una lista de las muestras con su número respectivo de observaciones, ordenadas de mayor a menor, por ejemplo YUN3533.2 fue la muestra con el mayor número de observaciones, 2152, mientras que YUN3259.1.1 fue la menor con un valor de 6, además vemos en color rojo aquellas muestras que quedarían fuera de nuestro análisis si hiciéramos una rarificación a 295 observaciones.
En la siguiente pestaña vemos una lista de la clave que identifica a cada feature, que en este caso son ASVs, ordenadas de mayor a menor frecuencia y el número de muestras en donde está presente, por ejemplo la primera de la lista es la feature que tiene la mayor frecuencia y está en 12 muestras.
El archivo de las secuencias representativas, tiene una tabla con el número total de features diferentes, en este caso son 1138 con una longitud en promedio de 227 nucleótidos. Las secuencias las puedes descargar en formato fasta. También vemos una tabla con la secuencia respectiva de cada feature, si damos click en la primera, nos abre un link al NCBI para hacer un blast, en esta página podemos mover los parámetros o utilizar los de default y dar click en ver reporte. Vemos en la página de resultados del Blast que identifica ésta feature como una bacteria no cultivada.
Microbioma lab es una start up de análisis bioinformáticos del 16S rRNA, ITS y metagenoma, así como servicio de secuenciación en Illumina, visita nuestro sitio web donde puedes ver el catálogo de productos y servicios.
En el archivos de los denoisng stats, lo que observamos en la primera columna es la clave de las muestras, en la segunda el número de secuencias que ingresaron de cada una al análisis, después en la tercer la columna están las secuencias que quedaron por muestra una vez que se hizo el filtrado de calidad, en la siguiente vemos el porcentaje de secuencias que pasaron los filtros de calidad, seguido del número de secuencias que quedaron después del denoising. En la siguiente columna están las secuencias que quedaron después de hacer el merge del forward con el reverse y el valor expresado en porcentaje. La columna de non-chimeric, muestra el número de secuencias que quedaron al final una vez que se eliminaron las quimeras y en la columna contigua lo vemos expresado en porcentaje.
Finalmente estos archivos, la tabla de frecuencias de las ASV y las secuencias representativas, son los que ya puedes utilizar en análisis posteriores como de diversidad alpha o beta, entre otros.
A partir de éste paso puedes continuar con los análisis de diversidad alpha y beta del tutorial de moving pictures.
En el siguiente video puedes ver la explicación de estos pasos.
Cursos de bioinformática disponibles
