Desarrollo del flujo de trabajo para filtrar taxones de microeucariotas de muestras tomadas por shotgun

# Desarrollo del flujo de trabajo para filtrar taxones de microeucariotas de muestras tomadas por shotgun --- El presente documento desarrolla los pasos y presenta las herramientas utilizadas de manera secuencial para el análisis de datos provenientes de muestras tomadas por metagenómica shotgun en tres estaciones del Océano Pacífico Tropical (OPTN) Nororiental en la Zona de Minimo Oxígeno (ZMO) frente a las costas de Colima. De cada programa utilizado se describe el archivo de entrada requerido y los archivos de salida generados. Cada uno de estos se analizó por separado, no se explica detalladamente en este documento pero sí se aborda la información rescatada de cada resultado en cada programa o paso. Los datos finales obtenidos con este flujo de trabajo, definido a continuación, se presentarán como resultados en la presentación de la evaluación de comité. La consulta de los script ejecutados en cada programa se puede realizar ingresando [aquí](https://github.com/SoFY-414/An-lisis-metagen-mico-de-microeucariotas-en-ZMO-ETNP/tree/main). --- ## Identificación de las muestras Las muestras de las 3 estaciones a analizar se divididen en dos Bioproyectos del NCBI (National Center for Biotechnology Information): [PRJNA632347](https://www.ncbi.nlm.nih.gov/bioproject/PRJNA632347) y [PRJNA254808](https://www.ncbi.nlm.nih.gov/bioproject/254808). Estas 3 estaciones son identificadas como **sta_1**, **sta_06** y **sta_2** en este informe. El Bioproyecto [PRJNA254808](https://www.ncbi.nlm.nih.gov/bioproject/254808) contiene todas las muestras de la sta_10 (4 muestras) y una parte de la sta_06 (5 muestras). Cada una de estas lecturas permanece guardada en un archivo de lectura de secuencias, identificado como *Sequence Read Archive* (SRA). | SRA | Estación | Muestras totales | | :---: | :---: | :---: | | SRR1509790, SRR1509792, SRR1509793, SRR1509794, SRR1509796| 06 | 5 | |SRR1509797, SRR1509798, SRR1509799, SRR1509800| 10 | 4 | El Bioproyecto [PRJNA632347](https://www.ncbi.nlm.nih.gov/bioproject/PRJNA632347) contiene todas las muestras de la sta_02 (7 muestras) y la otra parte de la sta_6 (8 muetras). Cada una de estas lecturas permanece guardada en un archivo de lectura de secuencias o *Sequence Read Archive* (SRA). | SRA | Estación | Muestras totales | | :---: | :---: | :---: | | SRR12426463, SRR12426461, SRR12426459, SRR11786911, SRR12426445, SRR12425871, SRR12425620| 02 | 7 | |SRR12431144, SRR12431036, SRR12430884, SRR12430814, SRR12426458, SRR11787837, SRR11787098, SRR12424656| 06 | 8 | Los archivos en conjunto de estos dos repositorios superan los >5 G por lo que la página del NCBI automáticamente sugiere la descarga y el uso de la herramienta SRA Toolkit, sin embargo, desde el portal del ENA (European Nucleotide Archive) se pueden buscar los bioproyectos con el mismo número de acceso: [PRJNA254808](https://www.ebi.ac.uk/ena/browser/view/PRJNA254808) [PRJNA632347](https://www.ebi.ac.uk/ena/browser/view/PRJNA632347) Para descargarlos es necesario seleccionar las muestras e indicar la generación de un archivo. Este archivo es un script que almacena las descargas de los archivos seleccionados y es ejecutable desde el servidor. Ejemplo: *ena-file-download-selected-files-20250917-2245.sh* --- ## Descripción de las muestras --- En total se recuperan 27 archivos asociados a las 27 muestras. Estos archivos contienen secuencias de extremos pareados (o paired-end reads). En metagenómica los extremos pareados son un método de secuenciación de ADN de alto rendimiento. En lugar de secuenciar el fragmento completo de principio a fin, la tecnología de secuenciación (como Illumina) lee una porción de la secuencia del extremo inicial (lectura 1 o Forward) y una porción del extremo final (lectura 2 o Reverse) del mismo fragmento. Como se puede observar en la imagen de abajo las 27 muestras están asociadas a metadatos que permiten ubicar su posición en la columna de agua de cada estación y en la ZMO del OPTN. ![Asignación taxonómica](https://hackmd.io/_uploads/rk-4MIY1-l.png) --- ## Descarga en conjunto de las muestras Consiste en recuperar el script de descarga que genera la misma plataforma del ENA (European Nucleotide Archive). Este script se descarga en el equipo de computo personal y posteriormente se sube a la cuenta del servidor del CIBNOR con ayuda del programa [Cyberduck](https://cyberduck.io/). Se ejecuta el script en la cuenta del servidor para iniciar con la descarga de las muestras seleccionadas. --- ## Ejecución de [FastQC](https://github.com/s-andrews/FastQC) El programa FastQC permite conocer y analizar la calidad de las lecturas en una muestra. Genera un informe descriptivo con gráficos y estadísticos. No modifica los archivos de lectura. --- #### Archivos de entrada Archivos crudos de extremo pareado de cada muestra en formato fastq (generalmente comprimido: ___.fastq.gz). #### Archivos de salida: 1. **Informe HTML** (*__fastqc.html). Formato web que se abre en un navegador y contiene gráficos, tablas e indicadores de Pasa/Falla/Advertencia sobre la calidad de las lecturas. 2. **Archivo ZIP** (*_fastqc.zip): Contiene los datos sin procesar del informe (gráficos e información en formato TXT). #### Script ejecutado Mediante el script **run_reads_fastqc.sh** se ejecuta el programa de manera automatizada para cada uno de los archivos. #### ¿Qué realizar después de la ejecución del programa FastQC? Los archivos `.html` los podemos descargar a nuestro equipo local con ayuda del programa Cyberduck y una vez descargados abrirlos para visualizar su contenido. Otra opción es descomprimir los archivos `.zip` y ver los resultados del reporte, esta manera no es tan visual pero podemos identificar las pruebas que pasaron y fallaron en el análisis con FastQC. --- ## Ejecución de [MultiQC](https://github.com/MultiQC/MultiQC) A diferencia del programa FastQC, el programa MultiQC realiza un solo reporte con todos los archivos que queremos analizar. También es un reporte con gráficas y estadísticos que nos permite comparar la información contenida en cada archivo. No modifica los archivos de lectura. --- #### Archivos de entrada: Archivos zip/txt generados por FastQC. #### Archivos de salida: 1. **Informe HTML** (multiqc_report.html): Reporte único con todos los resultados de todas las muestras. 2. Directorio de Datos (multiqc_data/): Contiene distintos archivos con datos brutos y metadatos utilizados para construir el informe HTML En la carpeta `reads_multiqc` se verifica que los archivos se hayan generado y contengan información válida, después se descarga el reporte `.html` al equipo local con ayuda de Cyberduck. #### Script ejecutado Mediante el script **run_reads_multiqc.sh** se ejecuta el programa de manera automatizada para cada uno de los archivos. #### **¿Qué analizar en el informe generado por MultiQC?** >tambien checar lo de la profundidad en el reporte 1. **Estadísticas generales (General Statistics)** Millones de lecturas por muestra, lontidud promedio de las lecturas, porcentaje de lecturas duplicadas o si falló un módulo. 2. **Calidad de la Secuencia por Base (Per Base Sequence Quality)** Muestra la distribución de las puntuaciones de calidad Phred (Q-scores) a lo largo de la longitud de las lecturas. La línea media debe mantenerse idealmente por encima de Q30 (99.9% de precisión de base) o Q20 en la mayor parte de la longitud. Es normal que la calidad caiga hacia el final de las lecturas (3' end). 3. **Contenido de Secuencia por Base (Per Base Sequence Content)** Muestra la proporción de las cuatro bases (A, T, G, C) por cada posición en las lecturas. Para una librería aleatoria (shotgun), las líneas de A, T, G y C deben ser casi paralelas y planas a lo largo de toda la lectura (idealmente entre 20-30% cada una). 4. **Secuencias Sobre-representadas (Overrepresented Sequences)** Esta prueba identifica secuencias que aparecen un número inusualmente alto de veces. Idealmente, esta lista debería estar vacía o mostrar secuencias con un porcentaje muy bajo. 5. **Contenido de Adaptadores (Adapter Content)** La contaminación por secuencias de adaptadores (que provienen de los kits de secuenciación) es uno de los problemas más serios y más comunes, ya que estas secuencias son artificiales y no biológicas. Verificar que el porcentaje de adaptadores detectados sea bajo o nulo. --- ## Ejecución de [Trimmomatic](https://github.com/usadellab/Trimmomatic) Trimmomatic es una herramienta que permite recortar y filtar aquellas lecturas de mejor calidad, así que es uno de los pasos escenciales ya que determinará los resultados de la asignación taxonómica y/o el ensamblaje. --- #### Archivos de entrada: Archivos FastQ de lecturas pareadas (Paired-End, PE)que pueden o no estar comprimidos. #### Archivos de salida: Para archivos de lecturas pareadas se generan 5 archivos. 1. Dos archivos en formato fastq (Pareados Sobrevivientes _1P y _2P). 1. Dos archivos en formato fastq (Lecturas Huérfanas/No Pareadas _1U y _2U). 1. Archivo de Log (un registro detallado de cada lectura recortada, indicando el nombre de la lectura, la longitud restante, y la cantidad de bases recortadas del inicio y del final). #### Parámetros modificados en Trimmomatic Una vez verificada la condición de las secuencias con MUultiQC antes del recorte y filtrado, y tomando en consideranción que las muestras provienen de un ecosistema marino de dificil acceso, se consideran los siguientes parámetros reportados en [Alexander, et al 2023](https://journals.asm.org/doi/10.1128/mbio.01676-23). Estos parámetros forman parte del script **(inserte nombre del script)** en donde: ~~~~ ILLUMINACLIP:TruSeq3-PE.fa:2:30:7 \ LEADING:2 \ TRAILING:2 \ SLIDINGWINDOW:4:2 \ MINLEN:50 ~~~~ Los anteriores parámetros representan la parte que le dice al programa cómo recortar y filtrar las secuencias. **Explicación** * `ILLUMINACLIP:TruSeq3-PE.fa:2:30:7` Elimina adaptadores de las lecturas usando un archivo de referencia (TruSeq3-PE.fa) que contiene las secuencias adaptadoras. * `TruSeq3-PE.fa` Archivo FASTA con las secuencias adaptadoras comunes de Illumina para datos paired-end. * `2` Número máximo de errores permitidos en la semilla de 16 bases para buscar el adptador/ alineación con la secuencia adaptadora (mismatches). * `30` Puntaje de corte para el alineamiento del adaptador. Si el alineamiento tiene score ≥ 30, se considera verdadero positivo y se recorta. * `7` Mínima longitud restante del fragmento (read) después del adaptador. Si lo que queda tras recortar es menor a 7 bases, se descarta la lectura. * `LEADING:2` Corta base de baja calidad al inicio de la lectura, 5', cuando no cumpla con el umbral de calidad ≥2. * `TRAILING:2` Corta bases de baja calidad al final de la lectura. Elimina desde el 3' hasta encontrar una base con calidad ≥2. * `SLIDINGWINDOW:4:20` Recorte dinámico de calidad usando una ventana deslizante: 4 tamaño de la ventana (número de bases): 20 umbral de calidad promedio. * `MINLEN:50` Descarta lecturas recortadas que sean demasiado cortas. `50` longitud mínima permitida después de todos los recortes (calidad y daptadores). Las lectaras más cortas se eliminan del archivo de salida. **¡Importante! :** Lecturas muy cortas pueden mapear de forma ambigua o ser útiles para ensamblado o binning. **¡Atención!** Antes de ejecutar Trimmomatic asegurarse de contar con el archivo `TruSeq3-PE.fa` o en dado caso descargarlo de la página de GitHub https://raw.githubusercontent.com/timflutre/trimmomatic/master/adapters/TruSeq3-PE.fa #### Script ejecutado Mediante el script **run_trimmo.sh** se ejecuta el programa de manera automatizada para cada uno de los archivos. #### ¿Qué realizar después de la ejecución del programa Trimmomatic ? Al igual que con los archivos de secuencias crudas, se realizan los reportes FastQC y MultiQC para comparar las secuencias después del recorte y filtrado con las secuencias crudas. Es importante comparar ambos reportes y poner especial atención a los apartados del reporte, ya mencionados en **Ejecución de MultiQC**, al igual que el valor total de las secuencias después del recorte y filtrado. Se rescata los dos archivos fastq con las secuencias pareados sobrevivientes (1P y 2P) al recorte y filtrado. * Poner un esquema de las lecturas sobrevivientes --- ## Ejecución de [MEGAHIT ](https://github.com/voutcn/megahit) Un ensamblador de genomas como MEGAHIT es un programa bioinformático diseñado para unir fragmentos cortos de ADN (llamados lecturas o reads) generados por tecnologías de secuenciación de nueva generación (NGS). Su objetivo es reconstruir la secuencia de ADN original de un organismo, produciendo secuencias más largas y continuas denominadas **contigs**. --- #### Parámetros modificados en MEGAHIT ``` MIN_CONTIG_LENGTH="" ##Se trabajó con 500 pb, 1000 pb y 3000 pb. ``` El anterior parámetro determina la longitud mínima de los contigs que se recuperarán depués del ensamble. #### Archivos de entrada: Archivos FastQ de lecturas pareadas (Paired-End, PE) que pueden o no estar comprimidos. #### Archivos de salida: Se generan 5 archivos y un directorio. 1. **Un archivo en formato fastq (final.contigs.fa)**. Contiene el resultado final del ensamblaje. Cada entrada en este archivo representa un contig, que es una secuencia de ADN continua y más larga reconstruida a partir de las lecturas cortas superpuestas Paired-End. Los encabezados de cada contig suelen contener información generada por MEGAHIT, como un identificador único, la longitud del contig (len) y la cobertura media (cov). 2. **Un archivo de texto (__.log)**. Contiene un registro detallado de todos los pasos que MEGAHIT ejecutó, incluyendo los parámetros utilizados, las estadísticas del proceso y cualquier advertencia o error. Se puede consultar la información del ensamblaje para cada muestra, por ejemplo: ``` 2025-08-22 16:34:03 - 14799 contigs, total 23678234 bp, min 1000 bp, max 30913 bp, avg 1599 bp, N50 1510 bp 2025-08-22 16:34:04 - ALL DONE. Time elapsed: 8427.881909 seconds ``` Esta información da una idea de la calidad de un ensamblaje genómico y es la información que se rescata para todas las muetras: el tiempo que demoró, el número de contigs, la longitud total de todas las secuencias de ADN combinadas, la longitud del contig más corto y del más largo en el ensamblaje, la longitud promedio de todos los contigs en el ensamblaje y la métrica N50. >Considerar esto para discusión en el reporte. N50: Te dice cuán grandes son los pedazos de ADN en los que está esa información.El N50 siempre es el mejor indicador de la calidad y continuidad de tu ensamblaje. Nos interesa recuperar un alto número de contigs y a su vez obtener 3. **Un archivo en un formato estructurado (JSON)**. Contiene todos los parámetros y opciones, tanto los que especificaste tú como los que MEGAHIT usó por defecto. Sirve para documentar y reproducir el ensamblaje de manera precisa. 1. **Un archivo vacío de nombre "done"** Funciona como una "bandera" o señal. Es un archivo de 0 bytes que indica que MEGAHIT terminó el proceso sin errores fatales. 1. **Un archivo checkpoints.txt** Contiene una lista de los puntos de control o etapas clave del ensamblaje que ya se han completado. Si el ensamblaje se detiene MEGAHIT leerá este archivo para saltarse los pasos que ya había terminado, ahorrando mucho tiempo y recursos computacionales. 1. **Un directorio intermediate_contigs**. Contiene varios archivos FASTA con los contigs generados en cada una de las iteraciones del algoritmo de MEGAHIT (que usa diferentes tamaños de k-mer). Son importantes para que el programa funcione, pero una vez que tienes tu ensamblaje final y has verificado que todo salió bien, este directorio se puede borrar para ahorrar espacio en disco, que suele ser considerable. #### ¿Qué realizar después de la ejecución del programa MEGAHIT ? Los archivos final.contigs.fa (un archivo por cada muestra) contendrá los contigs derivados del ensamblaje, pasarán a la asignación taxonómica. De los archivos .log generados se rescata la siguiente información. >esquema de la información analizada #### Script ejecutado Mediante el script **run_megahit.sh** se ejecuta el programa de manera automatizada para cada uno de los archivos. --- ## ASIGNACIÓN TAXONÓMICA La asignación taxonómica es el proceso de identificar el origen taxonómico de las secuencias de ADN (lecturas o reads) obtenidas de la muestra. Dado que la muestra contiene ADN de múltiples organismos, la asignación taxonómica busca responder a la pregunta: "¿A qué organismo pertenece cada fragmento de ADN secuenciado?" Esto permite determinar la composición de la comunidad microbiana (el microbioma). De acuerdo con la información encontrada y la documentación revisada, la asignación taxonómica puede realizarse de 3 distintas maneras: **a nivel nucletido** (CCMetagen), **a nivel proteína** (MetaEuk+Eukulele y EukMetaSanity), por la identificación de **genes marcadores** (Berrnap). Cada una de estas maneras tiene sus consideraciones, ventajas y desventajas, en algunas ocasiones se pueden admitir secuencias sin ensamblar o otras es mejor trabajas con secuencias ensambladas. Todos estos puntos se desarrollan y discuten en el informe y aquí solo se describen los programas utilizados. --- ## Ejecución de [CCMETAGEN](https://github.com/vrmarcelino/CCMetagen) --- CCMetagen (acrónimo de ConClave-based Metagenomics) es un flujo de trabajo de clasificación metagenómica **a nivel nucleótido**. Su objetivo principal es la identificación taxonómica completa y exacta de eucariotas y procariotas (incluyendo bacterias, hongos y otros microorganismos) presentes en una muestra de metagenoma o metatranscriptoma. Las bases de datos que se utilizan pueden contener desde genomas completos (RefSeq_bf.zip) o la colección completa de nucleótidos (nt_NCBI). Este programa contiene bases de datos listas para usarse y también esta la opción de crear una base de datos personalizada con las secuencias de los microorganismos de interés. Se utilizó la base de datos nt_NCBI disponible para su [descarga](https://figshare.unimelb.edu.au/projects/CCMetagen_databases/195755) y también se creo una base de datos a partir de los genomas disponibles que pertenecen a microorganismos marinos (29 genomas en total) (*realizar un gráfico con estos genomas y a qué grupo pertenecen) El programa CCMetagen se divide en dos partes principales: --- #### KMA La herramienta KMA (k-mer alignment) se ejecuta primero para identificar y cuantificar las lecturas que se alinean a las secuencias de referencia utilizando el algoritmo ConClave para mapeos precisos. ##### Parámetros modificados en KMA ``` ## Ruta completa a la base de datos de referencia (el índice de KMA) ref_db="" ``` La ruta de referencia anterior determina la ubicación de la base de datos será utilizada. Dentro de las bases de datos que el programa CCMetagen tiene a la disposición de ser descargadas y utilizadas se encuentra: NCBI(nt) y Refseq. Primero mapea y luego alinea las lecturas. 1. **Mapeo (Mapping)**: KMA utiliza k-mers (subcadenas cortas y fijas de la secuencia) para identificar rápidamente la plantilla o región de la base de datos más probable a la que pertenece cada lectura. Este es un paso de búsqueda rápida y selección. 1. **Alineación (Alignment)**: Una vez que se encuentra la plantilla candidata, KMA realiza una alineación exacta base por base. Esto determina la relación precisa entre cada base de la lectura y la secuencia de referencia, incluyendo mismatches, indels (inserciones/deleciones). Se consultaron las guías y espeficicaciones de la [página oficial](https://bitbucket.org/genomicepidemiology/kma/src/master/) para su aplicación en el flujo de trabajo. #### Archivos de entrada de KMA 1. **Archivos en formato fastq (Pareados Sobrevivientes _1P y _2P)**: secuencias sin ensamblar. 1. **Archivo en formato fastq (final.contigs.fa)**: secuencias ensambladas en contigs. #### Archivos de salida de KMA 1. **Archivo __.aln**: Alineamiento consenso, muestra como la secuencia consenso se alinea con la plantilla de referencia. Identificaicón de variaciones genéticas. 2. **Archivo __.frag.gz**:Reporte detallada sobre cada una de las lecturas individuales (nombre, a cuántas plantillas se parece, el puntaje asignado a la similitud, el punto de inicio y de fin de la alineacion de la platillas, nombre de la plantilla con la que alineo). 3. **Archivo __.fsa**: Contiene las secuencias consenso, la secuencia más probable que representa a todas y cada una de las lecturas. 4. **Archivo __.res**: Resultados con las estadísticas para cada mapeo. | Columnas en el archivo __.res | Significado | | -------- | -------- | | Nombre del templado| Text | | Score | Puntuacion ConClave de todas las lecturas que se aceptaron, a mayor puntaje mayor similitud. | |Template_length | Longitud total de la platilla o la referencia en la base de datos. | | Expected | Puntaje esperado si todas tus secuencias de ADN estuvieran distribuidas al azar en la base de datos. | | Template_Identity | Porcentaje que indica cuántas bases de ADN en la plantilla son exactamente igual a las de la secuencia. | | Template_Coverage | Porcentaje que indica qué tanto de la plantilla esta cubierta por la secuencia| | Query_Identity | Porcentaje que me indica qué tanto de mis bases son idénticas a las de la plantilla. | | Query_Coverage | Qué tanto de la secuencia está cubierta por la plantilla. | |Depth | Indica cuántas veces, en promedio, cada posición de la plantilla ha sido cubierta por las secuencias. 20X indica que cada base en la plantilla fue leída por 20 secuencias diferentes. | | Q_value | ¿Hay diferencia entre el score y el expected? ...Un valor alto indica que la diferencia es significativa. | | P-value | Un valor inferior a 0.05 indica que mi resultado no es aleatorio, casual o accidente. | #### Script ejecutado Mediante el script **run_kma.sh** se ejecuta el programa de manera automatizada para cada uno de los archivos. #### CCMetagen_py CCMetagen_py se encarga de las siguientes tareas esenciales: 1. **Lectura y Procesamiento de Resultados de KMA**: Lee el archivo de resultados (.res) generado por KMA. Este archivo contiene la información detallada de cómo se mapeó cada lectura de la muestra contra la base de datos de referencia. 2. **Filtrado de Calidad y Puntuación**: Aplica filtrado de calidad y de alineación adicional para refinar los resultados. Puede usar criterios como puntuaciones de alineación (scores) o cobertura para asegurar que solo las asignaciones de alta confianza sean consideradas. 3. **Asignación Taxonómica Jerárquica**: Utiliza la información de alineación para asignar la clasificación taxonómica final (desde el reino hasta la especie o cepa) a las lecturas. Realiza esta asignación utilizando el esquema de clasificación ConClave para resolver ambigüedades, lo que aumenta la precisión al considerar todos los mapeos posibles para una lectura. 4. **Generación de Reportes y Visualización**: Produce los archivos de salida finales que contienen la información clasificada. Genera un archivo de salida que puede ser usado directamente para crear gráficos Krona, una herramienta interactiva para visualizar la composición taxonómica de la muestra de forma jerárquica. #### Archivos de entrada de CCMetagen_py 1. **Archivo __.res** con las estadísticas de las secuencias mapeadas y alineadas con KMA. #### Archivos de salida de CCMetagen_py 1. **Archivo __.ccm.csv** (Comma-Separated Values - Valores Separados por Comas) : Este es el archivo de salida principal y más detallado de CCMetagen. Cada fila representa una secuencia de referencia mapeada (una contig o entrada del NCBI), y los valores asociados son las estadísticas de las lecturas que se mapearon a esa referencia. A continuación, se desglosa el formato y el significado de las columnas más importantes. | Columna | Ejemplo | Significado | | --- | ----------------------------------------------------------------- | ---------------------------------------------------------------------------------------------------------------- | | Closest_match | MG667112.1 Uncultured Synechococcus sp | Identificador de Referencia. | | Score | 29.0 | Número de Lecturas Clasificadas (Reads) que se asignaron a esta referencia. | | Expected | 0.0 | Número de Lecturas No Clasificadas o ambiguas. | | Template_length | 890.0 | Longitud de la Secuencia de Referencia (en nucleótidos). | | Template_Identity | 27.87 | Cobertura de la Secuencia (porcentaje de la secuencia de referencia que fue cubierta por las lecturas mapeadas). | | Template_Coverage | 28.09 | Profundidad de Lectura Promedio (Average Depth of Coverage). | | Query_Identity | 99.2 | Identidad de la Secuencia (Porcentaje de identidad promedio de las lecturas mapeadas). | | Query_Coverage | 356.0 | Score de Alineación (Una puntuación agregada de la calidad de la alineación). | | Depth | 0.28 | Porcentaje de Ambigüedad (Mide qué tan específica fue la asignación). | | q_value | 28.97 | Score de ConClave (Puntuación refinada usada por CCMetagen para la asignación final). | | p_value | 7.4e-08 | Valor E (E-value), indicando la significancia estadística de la alineación. | | LCA_TaxId | 154535 | ID Taxonómico (TaxID): El número oficial del NCBI del Último Ancestro Común (LCA). | | Superkingdom hasta Species | unk_sk, Bacillati, Cyanobacteriota, Cyanophyceae, Synechococcales |Jerarquía taxonómica. | 2. **Archivo __alignment.html** : Contiene la visualización jerárquica de la composición taxonómica de la muestra, gráficos Krona con los que se puede interactuar para ver la distribución de las lecturas a cada nivel. Es ideal para una exploración rápida y una visión general de la diversidad microbiana de la muestra. 4. **Archivo __alignment.tsv**(Tab-Separated Values - Valores Separados por Tabuladores) : Contiene información similar a la del archivo .ccm.csv, pero está diseñado específicamente para ser compatible con Krona. #### Script ejecutado Mediante el script **run_ccmetagen.sh** se ejecuta el programa de manera automatizada para cada uno de los archivos. --- ## Ejecución de [MetaEuk](https://github.com/soedinglab/metaeuk) + [Eukulele](https://github.com/AlexanderLabWHOI/EUKulele) La combinación de los programas MetaEuk y Eukulele permite realizar una clasificación metagenómica **a nivel proteína**, la principal diferencia con los clasificadores a nivel nucleótido es la traducción de las secuencias nucleotídicas a sus respectivos aminoáciácidos. En la clasificación a nivel proteína se pierden regiones no codificantes, es decir, únicamente se realiza la asignación taxonómica a aquellas secuencias identificadas como posibles proteínas putativas. Las bases de datos que se utilizan contienen proteínas predichas a escala genómica, esto quiere decir que la estructura primaria (la secuencia de aminoácidos) ha sido inferida computacionalmente a partir de datos genómicos (genomas o transcritos), sin necesidad de haberlas aislado o caracterizado experimentalmente. Las bases de datos son descargables desde el programa, mismo que les da una estructura para ser utilizadas. Se trabaja con la base de datos [MMETSP](https://zenodo.org/records/1212585#.Xw3PoJNKhTZ) ( Marine Microbial Eukaryotic Transcriptome Sequencing Project) porque se basa en muestras cultivadas de especies eucariotas marinas pelágicas y endosimbióticas que representan más de 40 filos. También la base de datos [EukProt](https://figshare.com/articles/dataset/EukProt_a_database_of_genome-scale_predicted_proteins_across_the_diversity_of_eukaryotic_life/12417881/3) se utiliza por contener una diversidad eucariota que actualmente incluye 993 especies de todos los principales supergrupos, así como taxones huérfanos. A diferencia de la base de datos MMETSP, ésta incluye organismos terrestres, marinos, unicelulares y pluricelulares. En su composición una parte significativa proviene de MMETSP y la otra de NCBI y otros repositorios. --- #### Ejecución de MetaEuk MetaEuk es el primer programa de este flujo de trabajo. Esta diseñado para identificar, localizar y anotar secuencias de ADN que son genes en grandes volúmenes de datos genéticos de organismos eucariotas, obtenidos directamente de muestras ambientales. Este programa incorpora un algoritmo que le permite manejar estructuras complejas (multi-exón) por lo que es ideal en el análisis de genomas eucariotas. MetaEuk toma los contigs y predice proteínas codificadas. Esto genera archivos .faa (proteínas) que luego se usan con EUKulele. #### Parámetros modificados en MetaEuk ``` metaeuk easy-predict \ --min-aln-len 40 \ --min-length 100 \ --cov-mode 2 -c 0.3 \ -e 1e-5 \ --alignment-mode 2 \ -s 7.5 \ --use-all-table-starts 1 \ --compressed 1 \ -v 3 \ --sort-results 1 \ --filter-hits 1 \ --remove-tmp-files 0 ``` ¿Qué significa cada parámetro? * `-s 7.5` Sensibilidad de la Búsqueda. Sensitivity: 1.0 faster; 4.0 fast; 7.5 sensitive [Predeterminado 4.000]. El valor 7.5 es el más alto ("sensitive"). Esto significa que el programa dedicará más tiempo a buscar alineamientos más débiles o distantes. Incrementa la posibilidad de encontrar homólogos, pero también aumenta significativamente el tiempo de ejecución y puede generar más falsos positivos. * `--min-aln-len 40`. Longitud Mínima del Alineamiento. Requiere que el alineamiento (la coincidencia de secuencia entre la consulta y la base de datos) tenga un mínimo de 40 aminoácidos. Esto es un umbral razonable para asegurar que los resultados sean biológicamente significativos y no solo coincidencias cortas al azar. * `-e 1e-5` Umbral de Valor E (E-value). Solo se listarán las coincidencias con un Valor E inferior o igual a $10^{-5}$. Este es un umbral estricto, lo que significa que solo se aceptarán coincidencias con una muy baja probabilidad de ser aleatorias. Este valor ayuda a reducir los falsos positivos a pesar de la alta sensibilidad (-s 7.5). * `--min-length 100` Longitud Mínima de Codones en ORF. Establece que los marcos de lectura abiertos (ORF) deben tener al menos 100 codones (300 nucleótidos) para ser considerados. 15 es el valor predeterminado. Esto filtra ORFs muy cortos que probablemente no codifican proteínas funcionales. * `--cov-mode 2 -c 0.3` Modo de Cobertura (--cov-mode 2): La cobertura se calcula únicamente sobre la secuencia de consulta (query). Cobertura Mínima (-c 0.3): Requiere que al menos el 30% de la secuencia de consulta esté cubierta por el alineamiento. Esto asegura que el alineamiento sea sustancial respecto a la longitud de la secuencia que está buscando. * `--alignment-mode 2` Cálculo del Alineamiento: Indica que, además de la puntuación y la posición final, se calculen la posición inicial y la cobertura para el alineamiento. Esto proporciona información más detallada sobre la coincidencia. * `--use-all-table-starts 1` Uso de Codones de Inicio Alternativos. En lugar de solo el codón estándar (ATG/AUG), se considerarán todos los codones de inicio alternativos definidos en la tabla de código genético que esté utilizando. MetaEuk solo se centra en la predicción de regiones codificantes (CDS) y, por lo tanto, no incluye regiones no codificantes como promotores o UTRs : #### Archivos de entrada: 1. **Archivo generado por el ensamblaje con MEGAHIT final.contigs.fa** (específico para cada muestra SRA). #### Archivos de salida: 1. **Archivo en formato fasta __.fas** : Alberga las secuencias de proteínas predichas. 2. **Archivo de formato fasta __.codon.fas** : Contiene las secuencias de codones correspondientes a las predicciones (CDS). 3. **Archivo __.headersMap.tsv** : Contiene el mapeo de encabezados (headersMap) en valores separados por tabulaciones (.tsv) . 4. **Archivo en formato GFF (General Feature Format) __.gff** : Contiene las anotaciones de genes y exones predichas. El archivo **__.fas** es el que pasa como archivo de entrada al siguente programa, Eukulele. Teniendo presente de que las secuencias predichas están escritas con los aminoácidos correspondientes de las secuencias de nucleótidos que fueron evaluadas por MetaEuk. El archivo **__.fas** contiene la siguiente información. Ejemplo: >A0A1W6BG25|k141_28420|-|620|1.256e-176|1|2|1030|1030[1030]:2[2]:1029[1029] | Ejemplo | Descripcción | | ---------- | ------------------------------------------------------------------------------------------------------------------------------------------ | | A0A1W6BG25 | Identificador de la proteína de referencia que se utilizó para encontrar esta predicción de gen. Las proteínas de referencia forman parte de la base de datos MMETSP | | k141_28420 | Identificador de la secuencia ensamblada (contig) donde se predijo el gen. | | - | Indica la cadena del ADN donde se encuentra el gen. Puede ser + (cadena positiva) o - (cadena negativa). | | 620 | Puntuación de similitud calculada entre la predicción y la proteína target de referencia. Un valor más alto indica una mejor coincidencia. | | 1.256e-176 | Indica la probabilidad de que una coincidencia con esta puntuación (o mejor) ocurra por casualidad. Los valores más cercanos a cero (como este 1.256e-176) indican una predicción altamente significativa y confiable. | | 1 | El número de exones predichos en el gen. En este caso, es un gen con un solo exón (predicción de gen completo o fragmento sin intrones). | | 2 | La coordenada más baja (más cercana al inicio) del gen predicho en el contig. | | 1030 | La coordenada más alta (más cercana al final) del gen predicho en el contig.| |1030[1030]:2[2]:1029[1029]| Estructura detallada del exón/gen. Para una predicción de un solo exón, como esta: 1) **1030[1030]**: Coordenada de fin del gen en el contig (y coordenada del fragmento original), 2). **2[2]**: Coordenada de inicio del gen en el contig (y coordenada del fragmento original) y 3) **1029[1029]**: Longitud del exón/gen en nucleótidos (bases). | #### Script ejecutado **run_metaeuk.sh** #### Ejecución de Eukulele El programa Eukulele tiene dos modos para ejecutar el análisis de las secuencias dependiendo de la naturaleza de los datos, el primero es **MAG** (Genomas ensamblados de metagenoma): son las proteínas predichas de un MAG derivado del agrupamiento de contigs similares de un ensamblaje de metagenoma basado en la frecuencia de tetranucleótidos, abundancia, etc. El segundo modo es **MET** (Metatranscriptomas): son contigs generados a partir de la secuenciación y el ensamblaje de datos metatranscriptómicos (ARN) de una comunidad mixta. Estos contigs se pueden proporcionar como secuencias de nucleótidos o secuencias de proteínas predichas de estos contigs. El modo que elegido es MET, teniendo siempre presente que los contigs con los que contamos no provienen de metatranscriptomas, porque el programa también permite ejecutar esta opcción en contigs derivados del metagenoma (no MAG). Si se desea analizar contigs metagenómicos como se describe para MET se recomienda encarecidamente proporcionar proteínas predichas. **¡Atención!** Dentro de la [documentación oficial de Eukulele](https://eukulele.readthedocs.io/en/latest/) se recomienda encarecidamente proporcionar proteínas predichas a partir de sus contigs metagenómicos, ya que este puede ser un proceso complejo (particularmente para metagenomas eucariotas) y no está dentro del alcance del programa. Las bases de datos, como se describío en un inicio, utilizadas son EukProt y MMETSP por la facilidad de información oficial descriptiva asociada a éstas bases. El programa contiene otras como PhyloDB y MMETSP+MarRef pero no pude acceder a información sobre su contenido o si están curadas. #### Parámetros modificados en Eukulele ``` --mets_or_mags mets \ --database eukprot \ --protein_extension .fas \ --proteins_or_nt proteins \ ``` * `--mets_or_mags` Modo mets para realizar análisis con contigs o metatranscriptomas y Modo mags para realizar análisis a genomas ensamblados de metagenomas (MAGs) * `--database` Se indica la base de datos a descargar y utilzar por el programa. MMETSP y EukProt fueron las elegidas. * `--protein_extension ` formato de las secuencias, si se encuentran como nucleotidos o aminoácidos. * `--proteins_or_nt` también para especificar el tipo de archivos de entrada. * #### Archivos de entrada: 1. **Archivo en formato fasta __.fas** : Secuencias de proteínas predichas. #### Archivos de salida: 1. **Directorio LOG** : Útil para: Verificar si hubo errores, ¿qué pasos se ejecutaron?, y ¿cómo se procesaron tus muestras? 3. **Archivo README_DB.txt** : Contiene información sobre la base de datos usada (por ejemplo, MMETSP). Útil para saber qué organismos, taxonomías o versiones están incluidas en la base de datos. 4. **Directorio taxonomy_estimation** : Contiene un archivo con los resultados de la asignación taxonómica por proteína. Este archivo es el que se analiza porque contiene todos los contigs que recibieron asignación taxonómica a nivel proteína. 5. **Directorio taxonomy_counts** : Contiene distintos archivos cada uno corresponde a un nivel de organización taxonómico y sus respectivos conteos. 6. **Directorio taxonomy_visualization** : Archivos para visualizar los resultados. No son muy claras las gráficas que se generan automáticamente (opté por generar y examinar los datos en Rstudio) El archivo **taxonomy_estimation** contiene la siguiente información: | Ejemplo | Columna | Significado | | ------- | ------- | ----------- | | A0A075FGF4|k141_38377|-|145|1.225e-33|1|2|298|298[298]:2[2]:297[297] | transcript_name | Se ubica de manera secuencial de izquierda a derecha: ID de la proteína en la base de datos, ID del conting donde se predijo la proteína, sentido de la lectura, longitud de la proteína, e-value del alineamiento, posiciones de alineamiento |inicio| fin| longitud| y la posición de la proteína dentro del contig. | | class | classification_level | El nivel taxonómico específico en el que se clasificó el contig. | | Eukaryota;Fungi;Zoopagomycota | full_classification | El desarrollo de la clasificación taxonómica completa. | | Zoopagomycota | classification | La clasificación específica obtenida. | | 70.9 | max_score | El porcentaje máximo de identidad notificado por el programa de alineación para la coincidencia. | | 0 | ambiguous | 1 si hubo varias coincidencias en desacuerdo para el contig, resueltas por anotación de consenso utilizando el límite proporcionado por el usuario (el valor predeterminado es 75%) o por los enfoques del último ancestro común (LCA), de lo contrario 0. | | EP00863 | database | Identificador en la base de datos | #### Script ejecutado **run_eukulele.sh** --- ## Ejecución de [EukMetaSanity](https://github.com/cjneely10/EukMetaSanity) El pipeline Run predice los loci génicos y los límites exón-intrón. Por cada genoma de entrada genera un conjunto de archivos esenciales para la anotación funcional posterior: **Archivos de entrada**: Contigs en formato FASTA. **Archivos de salida**: 1.- **Archivos de Formato de Búsqueda de Genes GFF/GFF3 (gene-finding format files)** : Estos archivos contienen las coordenadas de las predicciones de los loci génicos, exones e intrones. Archivos de salida específicos para cada herramienta individual (MetaEuk, Augustus, GeneMark-EP/ES). Archivos que contienen las predicciones consolidadas de las diferentes estrategias de fusión (Tier 1, Tier 2, o Tier 3), que reducen las anotaciones a un único conjunto de predicciones génicas por locus 2. **Archivos de Secuencia FASTA** de regiones Codificantes (CDS/cDNA): Contiene las regiones codificantes de genes (gene coding regions) del ADN genómico. 3. **FASTA de Secuencias de Proteínas Derivadas**: Contiene las secuencias proteicas derivadas (derived protein sequences). 4. **Archivos Relacionados con Repeticiones** : Archivos generados por la identificación y enmascaramiento de repeticiones (por ejemplo, archivos de repeticiones). El pipeline Report utiliza las salidas de Run para realizar la anotación funcional: **Archivos de entrada**: Secuencias de Proteínas Predichas: Los genes predichos en la etapa Run. **Archivos de salida**: Archivos de resumen delimitados por tabulaciones (tab-delimited summary files) para cada llamada génica. #### Script ejecutado **run_eukmetasanity_local.sh** No se completó el análisis con este programa. --- ## Ejecución de [Berrnap](https://github.com/tseemann/barrnap) Programa de alineación a nivel gen marcador