Skip to content

Home

Jump to bottom
Silvia Talavera Marcos edited this page Apr 22, 2022 · 1 revision

Scripts desarrollados

Este Anexo recoge los scripts desarrollados para automatizar los análisis correspondientes al Trabajo de Fin de Máster "Modelización metabólica de comunidades microbianas estables crecidas con fuentes de carbono y energía únicas y simples" (Silvia Talavera Marcos, Máster en Bioinformática y Biología Computacional, Universidad Autónoma de Madrid, 2020-2021).

Esta wiki también está disponible en PDF y como archivo HTML (formato mejorado).

A continuación se incluye una breve descripción de todos ellos:

  • modelado.R incluye el alineamiento con Nucmer, la creación de modelos con CarveMe y (opcionalmente) el análisis con Smetana para una pareja de nodos dada.

  • annotate.R se encarga de la anotación funcional con eggNOG-mapper y, opcionalmente, la creación de un modelo consenso. También es capaz de llamar a Nucmer para iniciar el proceso desde el principio.

Las funciones definidas para estos scripts se encuentran en utils.R.

  • RefrFBA.py es un script que se encarga de ejecutar simulaciones de FBA para todos los modelos de un PCG dado. Devuelve las tasas de crecimiento en forma de archivos .csv y por la salida estándar.

  • parser.R se encarga de generar informes en formato de texto plano que resumen los resultados de Smetana.

  • Por último, se han desarrollado dos scripts en Python para crear los consensos: consenso.py crea un modelo SBML a partir de otros modelos SBML y consenso_EGG.py crea una tabla de anotaciones consenso a partir de múltiples archivos de anotaciones.

Todos los scripts aquí incluidos se pueden encontrar también en el repositorio https://github.com/silvtal/TFM. En todos los casos se incluye el código del script seguido de una explicación del funcionamiento del mismo.

modelado.R

Uso: modelado.R [options]

Argumentos:

Este script, como todos los demás, se ejecuta llamándolo desde el terminal. Los argumentos se recogen con la función OptionParser de la librería optparse. La opción --help devuelve una explicación en inglés de todos ellos.

  • -n, --nodes. Nombre del archivo de texto con información sobre los dos nodos de interés (e. g. 'my_node_data.txt'). El formato del archivo debe ser el que se muestra a continuación. La taxonomía debe incluir al menos dos filos, correspondiendo el segundo al filo:

      Node35562	k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacteriales;f__Enterobacteriaceae;

  Node27828	k__Bacteria;p__Proteobacteria

  • -m, --medium. Medio (e.g. 'M9', 'M9[glc]').

  • --mediadb. Archivo con la base de datos de medios, por defecto "my_media.tsv" incluido en la carpeta original del script.

  • -c, --checking. Si TRUE, los modelos generados y analizados se limitan a aquellos pares de especies que coinciden en una o más muestra de un experimento especificado (corresponde a la Estrategia 1). FALSE (Estrategia 2) por defecto.

  • -e, --experiment. Nombre del archivo de texto con información del experimento de interés (e.g. 'table.from_biom.tsv'). Es un archivo separado por tabulaciones con el siguiente formato:

      #OTU ID	S1	S2	S3	(...)

  O1	000000	000001	000002	(...)

  • --coupling. Si TRUE, Smetana hace un análisis adicional sin su opción --no-coupling.

  • --nucmer. Ruta del ejecutable de Nucmer (e.g. './MUMmer3.23/nucmer', '~/my_apps/nucmer').

  • --showcoords. Ruta del ejecutable del show-coords de MUMmer (e.g. './MUMmer3.23/show-coords', '~/my_apps/show_coords').

  • --tree. Nombre del archivo del árbol filogenético 16S. Los nombres de los nodos deben modificarse para que coincidan con los del archivo de información de nodos (-n, --nodes). Esta modificación se puede hacer fácilmente por ejemplo usando la función makeNodeLabel del paquete ape de R. El árbol por defecto es '99_otus_nodes.tree' (Greengenes gg_13_5), en la carpeta original del script.

  • --fasta. Archivo multifasta de secuencias 16S de las hojas del árbol utilizado. El archivo por defecto es'99_otus_nodes.tree' (Greengenes gg_13_5) en la carpeta original del script.

  • --db16s. Base de datos de secuencias 16S contra la cual Nucmer va a alinear las secuencias de las hojas del árbol. Por defecto se usa 'bac120_ssu_reps_r95.fna' (GTDB https://data.ace.uq.edu.au/public/gtdb/data/releases/release95/95.0/genomic_files_reps/).

  • --dbproteins. Base de datos de secuencias de aminoácidos. CarveMe tomará de aquí los archivos que correspondan a los hits de Nucmer y creará modelos metabólicos a partir de ellos. Por defecto es 'protein_faa_reps/bacteria/' (GTDB https://data.ace.uq.edu.au/public/gtdb/data/releases/release95/95.0/genomic_files_reps/).

  • --run_smetana. Si TRUE, ejecuta un análisis (global y detallado) de Smetana para parejas internodo. Si FALSE, se salta este paso.

  • --cores. Número de núcleos usados en procesos de paralelización (mclapply y mcmapply). Por defecto son 4.

Descripción: modelado.R se encarga de generar con CarveMe modelos de los nodos/PCGs especificados en la entrada (--nodos) y, opcionalmente, los análisis de Smetana correspondientes.

En primer lugar, si --checking es TRUE, se seleccionan las hojas de cada nodo que sí coincidían con las hojas del otro nodo en al menos una muestra del experimento especificado (guardado en la variable exp). Se hace con la función check definida en utils.R.

En segundo lugar, se hace un alineamiento con Nucmer de las secuencias 16S (GreenGenes) de las hojas deseadas para cada nodo y una base de datos de genomas completos, GTDB. Solo se seleccionan aquellas asignaciones que superan un filtro fijo (identidad > 97% y cobertura (de la referencia/hit) > 90%). A las hojas cuyo alineamiento pase el filtro se les asigna un archivo con una lista completa de todas las proteínas de su genoma. Ese archivo es a partir del cual CarveMe generará modelos SBML. La función principal de este proceso es find_alignment_hits.

En tercer lugar, se generan los modelos con CarveMe. Los modelos se guardan con el nombre de la hoja correspondiente. Es uno de los pasos más largos, pero está paralelizado con mcapply para ahorrar tiempo de ejecución. Además, si ya existe un modelo con el nombre de esa hoja no se genera de nuevo, permitiendo ahorrar tiempo si la ejecución de modelado.R se ha reiniciado. Cada nodo tiene su propia carpeta con el nombre del nodo (el que se indique en el archivo dado con --nodos) que está dentro, a su vez, de la carpeta "models". Esta carpeta se guarda automáticamente en el directorio de trabajo desde el que se ha ejecutado el script (que no tiene por qué coincidir con el directorio donde se encuentra el mismo).

Por último, se ejecuta el paso de Smetana, run_smetana==TRUE. Primero se crea una lista con todas las parejas inter-nodo posibles de modelos o, si checking==TRUE, solo aquellas parejas inter-nodo que coincidieran en la misma muestra del experimento. A continuación hay, para cada pareja, dos o tres ejecuciones según si coupling es TRUE o no. En cada ejecución, Smetana va comprobando a) que los modelos SBML existen y b) que el análisis no existiera ya. Asimismo, estas ejecuciones están paralelizadas con mcapply. Los análisis se guardan en la carpeta "smetana_results", que también se crea en el directorio de trabajo, con subcarpetas para cada tipo de análisis.

Durante la ejecución se imprimen en la salida estándar mensajes que indican la actividad del script. Al final se imprime el tiempo de ejecución.

Salida:

  • Modelos SBML para OTUs/hojas de los nodos/PCGs especificados.
  • Si --run_SMETANA es TRUE, análisis globales y detallados de Smetana.
  • Si --coupling es TRUE, se generan también informes detallados sobre los metabolitos acoplados al crecimiento.
  • También se generan archivos con el output del alineamiento de Nucmer, incluyendo un archivo queries_v_hits para cada PCG.

Requiere: R (3.5.0), CarveMe (1.4.0+), Smetana (1.2.0), NUCmer (cualquier versión; usamos la 3.1; MUMmer 3.23), ape (se instala automáticamente), parallel, optparse

annotate.R

Uso: annotate.R [options]

El input principal es bien --nodes o bien --genomes.

Argumentos:

  • -g, --genomes. Ruta de un archivo de texto que contenga una lista de genomas a anotar. El formato debe ser el siguiente, en consonancia con el de la base de datos de genomas de proteínas dada:

      RS_GCF_000281895.1

  RS_GCF_900100495.1

  RS_GCF_900104015.1

  (...)

  • -n, --nodes. Ruta de un archivo de texto con información de los nodos de interés. El formato debe ser el siguiente:

      Node35562	k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacteriales;f__Enterobacteriaceae;

  Node27828	k__Bacteria;p__Proteobacteria

  • --skip_consensus. Si TRUE, eggNOG-mapper crea los archivos anotados pero no se crean modelos ni un consenso a partir de los mismos. Por defecto es FALSE.

  • -m, --medium. Medio (e. g. 'M9', ' M9[glc]').

  • --outputdir. Directorio de salida para los archivos de anotación y los modelos consenso.

  • -o, --outputname. Nombre del modelo consenso creado (si el input es -g) o nombre común de los modelos consenso (si el input es -n y son varios nodos).

  • --mediadb. Ruta del archivo de medios de cultivo.

  • --checking. Si TRUE, solo se anotan aquellos genomas cuyas secuencias hayan sido detectadas en el experimento especificado. Solo se aplica si se hace alineamiento.

  • -e, --experiment. Ruta del archivo con información del experimento de interés cuando --checking es TRUE (e. g. 'table.from_biom.tsv'). Solo es necesario si se hace alineamiento y --checking es TRUE. El archivo debe tener el siguiente formato:

      #OTU ID	S1	S2	S3	(...)
  O1	000000	000001	000002	(...)

  • --nucmer. Ruta del ejecutable de Nucmer (e. g. './MUMmer3.23/nucmer', '~/my_apps/nucmer'). Solo es necesario si se hace alineamiento.

  • --show-coords. Ruta del ejecutable de show-coords (e. g. './MUMmer3.23/show-coords', '~/my_apps/show_coords'). Solo es necesario si se hace alineamiento.

  • --emapper_path. Ruta del ejecutable de eggNOG-mapper (e. g. './eggnog-mapper-my_version/emapper.py').

  • --tree. Nombre del archivo del árbol filogenético 16S. Solo es necesario si se hace alineamiento. Los nombres de los nodos deben modificarse para que coincidan con los del archivo de información de nodos (-n, --nodes). Esta modificación se puede hacer fácilmente por ejemplo usando la función makeNodeLabel del paquete ape de R. El árbol por defecto es '99_otus_nodes.tree' (Greengenes gg_13_5), en la carpeta original del script.

  • --fasta. Archivo multifasta de secuencias 16S de las hojas del árbol utilizado. Solo es necesario si se hace alineamiento. El archivo por defecto es'99_otus_nodes.tree' (Greengenes gg_13_5) en la carpeta original del script.

  • --db16s. Base de datos de secuencias 16S contra la cual Nucmer va a alinear las secuencias de las hojas del árbol. Solo es necesario si se hace alineamiento. Por defecto se usa 'bac120_ssu_reps_r95.fna' (GTDB https://data.ace.uq.edu.au/public/gtdb/data/releases/release95/95.0/genomic_files_reps/).

  • --dbproteins. Base de datos de secuencias de aminoácidos. EggNOG-mapper tomará de aquí los genomas a anotar, ya sean dados por los hits de Nucmer (--nodes) o por el usuario (--genomes). Por defecto es 'protein_faa_reps/bacteria/' en la carpeta original del script (GTDB https://data.ace.uq.edu.au/public/gtdb/data/releases/release95/95.0/genomic_files_reps/).

  • --dmnd_db. Ruta de una base de datos de Diamond personalizada. Especialmente útil cuando se utiliza una base de datos diferente a la estándar de eggNOG-mapper, como puede ser una creada con create_compatible_database.R.

  • --cores. Número de cores a usar en la paralelización de procesos (mclapply). Por defecto son 4.

Descripción: En primer lugar se obtiene la lista de genomas a anotar. Esta lista puede tomarse directamente de los datos de entrada (archivo de texto dado con la opción --genomes) o bien realizando un alineamiento del nodo o los nodos especificados con --nodos. El proceso en este caso es similar al explicado para modelado.R, utilizando la función find_alignment_hits.

En segundo lugar, cada lista de genomas es anotada con eggNOG-mapper mediante la función annotate definida en utils.R. Este proceso está paralelizado internamente con mcapply.

Por último, si se activó la opción de crear consenso, se hace una llamada a consenso_EGG.py y finalmente a CarveMe.

Salida: Archivos de anotaciones de genomas para cada genoma de interés. Si --skip_consensus es FALSE, también se genera un modelo SBML-FCB2 consenso para cada nodo.

Requiere: R (3.5.0), CarveMe (1.4.0+), eggNOG-mapper (2.0.1), parallel, optparse. Si se utiliza la función de alineamiento con Nucmer, se requieren también NUCmer (cualquier versión; usamos la 3.1; MUMmer 3.23) y ape (se instala automáticamente).

create_compatible_database.R

Al utilizar eggNOG-mapper junto a otra herramienta que use Diamond (como CarveMe), o simplemente por tener otra versión de Diamond instalada, es posible que se importe al entorno una versión de Diamond diferente a la que viene con eggNOG-mapper. En ese caso, podría haber incompatibilidad entre dicha versión de Diamond y la base de datos de Diamond incluida con el ejecutable de eggNOG-mapper (<emapper_path>/data/eggnog_proteins.dmnd) y consecuentemente un fallo a la hora de anotar.

Este breve script sirve para obtener una base de datos Diamond en un formato distinto a la utilizada por eggNOG-mapper y compatible con la versión de Diamond que se esté utilizando o se quiera utilizar. Posteriormente se puede indicar a annotate.R que utilice esta nueva base de datos mediante --dmnd_db.

Uso: create_compatible_databases.R [options]

Argumentos:

  • -o, --outputdir. Ruta del directorio de salida para la nueva base de datos.

  • -e, --emapper_path. Ruta de la versión de eggNOG-mapper (contiene la versión antigua de Diamond y su base de datos antigua).

  • -d, --diamond. Versión de Diamond para la nueva base de datos. Si no se especifica, se usa la versión del PATH.

Descripción: Este script toma los datos de entrada y ejecuta un primer comando que convierte la base de datos antigua, en formato .dmnd, en un archivo .faa.

El segundo y último comando que ejecuta llama a la nueva versión de Diamond para crear, a partir del archivo .faa, una nueva base de datos en formato .dmnd. El archivo .faa se guarda también en la carpeta de salida.

Salida: /eggnog_proteins.faa y /eggnog_proteins_compatible.dmnd

Requiere: eggNOG-mapper (cualquier versión), Diamond (cualquier versión)

utils.R

Uso: -

Argumentos: -

Descripción: Este script recoge las funciones más importantes de modelado.R y annotate.R. nucmer, carve, emapper y smetana lanzan estos 4 programas desde el terminal, preparando el comando a lanzar en cada caso e incluyendo en ocasiones mensajes de error si falta algún argumento. annotate es una función sencilla que se encarga de llamar a emapper para cada genoma.

Las funciones carve y smetana son un poco más complejas. En primer lugar, dado que deben ejecutar cientos o incluso miles de comandos según el caso (uno por cada modelo o análisis generado), los comandos están paralelizados con mclapply y mcmapply. Además, en el caso de carve, antes de ejecutar cada comando se comprueba si ya existe un modelo con ese nombre (generado en alguna ejecución previa de modelado.R), para evitar dedicar tiempo computacional a generar el mismo modelo dos o más veces. En el caso de smetana se comprueba antes de cada comando que no se haya generado anteriormente un análisis con el mismo nombre, por la misma razón.

smetana se encarga además de implementar o no un tercer análisis, el de los metabolitos acoplados al crecimiento, dependiendo de si la opción de coupling está activada.

La función gram es llamada desde carve. Esta función se encarga de, dado un filo taxonómico, devolver si es Gram-positivo o Gram-negativo. Así, carve utilizará el modelo universal más adecuado en cada caso. Se basa en dos archivos que contienen breves listas de filos Gram-positivos y Gram-negativos y que se encuentran en la carpeta original del script.

check es la función encargada de hacer la selección de modelos para la Estrategia 1. En primer lugar obtiene todas las posibles combinaciones de parejas entre el par de nodos dado. A continuación, recorre con un bucle todas las muestras del archivo del experimento dado y selecciona aquellas parejas que coinciden en al menos una.

find_alignment_hits es la función principal encargada del alineamiento, llamada en modelado.R y en annotate.R. Lo que hace es ejecutar nucmer para cada una de las hojas de cada nodo dado y aplica el filtro de calidad de identidad > 97% y cobertura (de la referencia/hit) > 90%. Devuelve un vector con el nombre de cada hoja (cuyo alineamiento haya superado el filtro de calidad) y su correspondiente hit. Este vector es guardado a modo de tabla en un archivo, "queries_v_hits". También se generan los archivos "nucmer_unfiltered" y "nucmer_filtered", con la salida de Nucmer antes y después del filtro de cobertura, y el archivo "genomes", que consiste en una lista de todos los hits que han pasado el filtro. Este último es el tipo de archivo que se da como entrada a annotate.R usando --genomes.

Salida: -

Requiere: R (3.5.0), optparse, parallel. nucmer requiere el programa Nucmer (siendo válida cualquier versión; usamos la 3.1; MUMmer 3.23). smetana requiere Smetana (1.2.0). carve requiere CarveMe (1.4.0+). emapper requiere eggNOG-mapper (2.0.1).

RefrFBA.py

Uso: RefrFBA.py input_folder mediadb medium outputdir outputname

Los argumentos se indican sin flags en este script.

Argumentos:

  • input_folder. Una carpeta que contenga todos los modelos SBML-FBC2 que se quiera analizar. Puede ser, por ejemplo, una de las carpetas generadas con modelado.R.

  • mediadb. Base de datos de medios de cultivo a utilizar.

  • medium. Medio en el que simular el crecimiento de los modelos dados.

  • outputdir. Directorio de salida para el informe generado.

  • outputname. Nombre del informe generado (incluyendo extensión .csv).

Descripción: En primer lugar se hacen los preparativos: a partir de la carpeta dada como input_folder, se crea una lista que contiene todos los modelos, filtrando por extensión .xml. También se carga el medio de cultivo y se inicializa como "entorno" con el método clase Environment.from_compounds, importado de CarveMe.

A continuación, se inicia un bucle en el que cada modelo de la lista será cargado con la función de CarveMe load_cbmodel, se le aplicará el medio de cultivo y se incorporará a una simulación FBA con la función de ReFramed FBA. De cada simulación se guarda la tasa de crecimiento obtenida y, si es mayor que 0, se incorporará al cálculo de la tasa media de crecimiento para el total de los modelos. Todas las tasas de crecimiento se guardan en una lista de listas, results.

Por último, se imprime por la salida la tasa media de crecimiento y la desviación típica. Los datos de results se convierten a DataFrame de pandas y se guardan como .csv con el nombre indicado.

Salida: Archivo .csv con todas las tasas de crecimiento junto al nombre de cada modelo. Salida estándar: tasa media de crecimiento y desviación típica.

Requiere: Python (3.6+; en este Trabajo se usa 3.6.9). CarveMe (1.4.0+; incluye las funciones importadas, Pandas y ReFramed).

RefrFBA_E1.py

El script RefrFBA.py tiene una versión alternativa, RefrFBA_E1.py, que recibe como input una carpeta ligeramente distinta, con la siguiente estructura:

input_folder
├── cit
│   ├── Node28866
│   │   ├── 4370747.xml
│   │   └── (...)
│   └── Node35562
│       ├── 3944484.xml
│       └── (...)
├── glc
│   ├── Node27828
│   │   └── (...)
│   └── Node35562
│       └── (...)
└── leu
    ├── Node13821
    │   └── (...)
    └── Node28853
        └── (...)

RefrFBA_E1 funciona de la misma manera, solo que el bucle va cambiando el medio especificado automáticamente para cada subcarpeta (M9[cit], M9[glc] o M9[leu]). El .csv final incluye todos los datos juntos, especificando el medio de cultivo y añadiendo también una columna con la tasa de crecimimento media. Se incluye el código en el desplegable a continuación.

parser.R

Uso: parser.R -d/-g input_file -c cores --report_name reportname

Argumentos:

  • -g, --global. Carpeta que contenga informes del modo global de Smetana (e. g. ./my_results/NodeXXXXX/smetana_results/global).

  • -d, --detailed. Carpeta que contenga informes del modo detallado de Smetana (e. g. ./my_results/NodeXXXXX/smetana_results/detailed).

  • -c, --cores. Número de núcleos que usar en la paralelización.

  • --report_name. Ruta y nombre del archivo de salida. Por defecto es ./report.txt.

Descripción: Una vez leída la carpeta de entrada, se sigue un proceso distinto dependiendo de si es del modo global o del detallado.

Si es del modo global, primero se unen todos los archivos en un solo data.frame, leyéndolos con mclapply para mayor velocidad. Después se filtran y guardan (en not_growing) aquellas filas que contengan una MRO no disponible. Los nombres de las filas del data.frame coinciden con los de los modelos. Si la MRO es distinta de 0 en al menos un caso, el script continúa y obtiene el valor máximo, el medio y el mínimo de la MRO y la MIP. Finalmente se crea un informe que contiene una lista de aquellos archivos donde la MRO no estaba disponible, los valores de MRO y MIP y finalmente una lista completa de todos los resultados.

Si la carpeta dada es de resultados del modo detallado, primero se leen todos los archivos y se guardan en una lista. Esta lista se divide en dos una conteniendo los informes de Smetana vacíos y otra los que no están vacíos. Estos dos primeros pasos están paralelizados con mclapply.

Como los nombres predeterminados de los informes de Smetana (generados con modelado.R) incluyen los nombres de los dos modelos, uno de cada nodo, este script crea dos listas, first y second, que contienen los nombres de todos los modelos que pertenecen a cada uno de los nodos según su localización en el nombre del archivo. Esto se hace con la función check_node. A continuación se agrupan los metabolitos intercambiados de cada informe según el sentido en el que se intercambian, es decir, si van desde el nodo first hacia el second o viceversa (all.metabolites.by.node). También se crea una lista en la que se agrupan los metabolitos independientemente del sentido de intercambio (all.metabolites). Estas listas se ordenan de mayor a menor puntuación Smetana y se reportan en el informe final. También se incluye en el informe qué informes estaban vacíos.

Salida: Informe en formato de texto plano.

Requiere: R (3.5.0.), optparse, parallel

consenso.py

Uso: consenso.py input_folder output_folder outputname

Argumentos:

  • input_folder. Una carpeta que contenga todos los modelos SBML-FBC2 a partir de los cuales se desee generar un consenso. Puede ser, por ejemplo, una de las carpetas generadas con modelado.R.

  • output_folder. Carpeta de salida del consenso.

  • outputname. Nombre del consenso (+ consensus.xml).

Descripción: Este script genera un modelo consenso en formato SBML-FBC2 a partir de archivos en el mismo formato. Está optimizado para modelos que no hayan sido sometidos a gap-filling, ya que esto modifica los campos del SBML. El formato SBML es un tipo de XML con los siguientes campos:

  • root[0][0] --> notas del formato (notes)

  • root[0][1] --> lista de compartimentos (listOfCompartments)

  • root[0][2] --> lista de compuestos (listOfSpecies)

  • root[0][3] --> lista de parámetros (listOfParameters)

  • root[0][4] --> lista de reacciones (listOfReactions)

  • root[0][5] --> lista de objetivos (listOfObjectives)

  • root[0][6] --> lista de productos génicos (listOfGeneProducts)

Lo que hace este script es fijarse en los campos 2 (compuestos), 4 (reacciones) y 6 (productos génicos) e incluir en un archivo XML final (el consenso) todas las entradas de esos campos que estén en más de un 80% de modelos. Las reacciones del campo 4 es el que determina el crecimiento de cada modelo metabólico, mientras que el campo 6 es un campo propio del formato FBC2 que define productos génicos únicos para cada gen. Como los campos 2 y 6 dependen de las reacciones, también los tenemos en cuenta.

En primer lugar se crea un modelo XML vacío que hace de esqueleto del consenso. Este modelo está formado por todos los campos del primer modelo de la carpeta de entrada, excepto la reacción de biomasa (crecimiento) del campo 4. Paralelamente, creamos una variable (growth) donde se guardará, aparte, todas las reacciones de biomasa de todos los modelos con sus reactivos y coeficiente originales. Empezamos guardando la del primer modelo. También se crea una tercera variable, conteo, que es una lista de diccionarios donde se apuntará cuántas veces aparece cada reactivo (campo 2), cada reacción que no sea la de crecimiento (campo 4) y cada producto génico (campo 6).

En segundo lugar, un bucle va recorriendo los demás modelos de la carpeta de entrada. A cada modelo abierto se le extrae la reacción de crecimiento y se guarda en growth y posteriormente se recorren sus campos 2, 4 y 6. Los campos 0 (descripción del modelo), 1, 3 y 5 no se cambian. Cada elemento e de cada campo (reactivo, reacción o producto génico en cada caso) se cuenta en el diccionario de conteo. Si no estaba ya incluido en el consenso, se incluye. Cada modelo se elimina de la memoria después de ser recorrido en el bucle.

Finalmente, un segundo bucle recorre los campos 2, 4 y 6 del consenso y elimina aquellas entradas cuyo número de apariciones en conteo sea menor del 80% del total de modelos. Por último, se selecciona la función de crecimiento más común entre todos los modelos y se añade al consenso final, que se guardará como .xml.

Salida: Modelo consenso en formato SBML-FBC2.

Requiere: Python (3.6+; en este Trabajo se usa 3.6.9).

consenso_EGG.py

Uso: consenso_EGG.py input_folder output_folder (outputname) (percentage)

Argumentos:

  • input_folder. Una carpeta que contenga todos los archivos de anotaciones de eggNOG-mapper a partir de los cuales se desee generar un consenso. Los archivos de anotaciones deben contener annotations en el nombre.

  • output_folder. Carpeta de salida del consenso.

  • outputname. Nombre del consenso (+ .tsv).

  • percentage. Porcentaje de archivos en los que debe estar presente una anotación funcional para ser incluida en el consenso final. Por defecto es 80%. Debe indicarse con el formato "0.80".

Descripción: Este script genera un archivo consenso de anotaciones funcionales a partir de una carpeta dada que contenga múltiples archivos de anotaciones funcionales, obtenidos con eggNOG-mapper. El formato de estos archivos es variable según las opciones indicadas a eggNOG-mapper; el formato compatible en este caso es aquel que se genera con annotate.R y que CarveMe es capaz de leer con su opción --egg para generar modelos.

En primer lugar, se lee la carpeta de entrada. Si hay al menos un modelo con annotations en el título, el script continúa. En segundo lugar, se inicializa con la función load_eggnog_dataun DataFrame de pandas, llamado all_reactions, a partir del primer modelo de la lista. También se inicializa un diccionario de conteo de reacciones.

all_reactions incluirá todas las reacciones diferentes que aparezcan en los modelos de la carpeta de entrada y se llena en un bucle que recorre los demás modelos. Se abre cada modelo, se seleccionan sus reacciones únicas y se contabilizan en conteo. Como cada reacción tiene una puntuación asignada, utilizada por CarveMe para hacer modelos, cada vez que se repite una reacción se actualiza su puntuación haciendo una media ponderada. También se van anotando los índices de cada reacción para facilitar los pasos posteriores. Tras cada iteración, se elimina de memoria el último modelo abierto.

Por último, se aplica el filtro: se eliminan de all_reactions aquellas reacciones que no aparezcan en al menos el porcentaje especificado de archivos de anotaciones originales. La lista filtrada se guarda como .tsv. Posteriormente esta lista de anotaciones consensuada podrá ser utilizada por CarveMe para generar un modelo consenso.

Salida: Archivo de anotaciones funcionales .tsv.

Requiere: Python (3.6+; en este Trabajo se usa 3.6.9), CarveMe (1.4.0+, incluye pandas).

Clone this wiki locally