Benzo[a]pireno y artritis reumatoide: una investigación computacional integrada

Resumen

Este estudio empleó un enfoque integrado de toxicología computacional para investigar sistemáticamente la relación entre la exposición al benzo[a]pireno y la artritis reumatoide. El análisis identificó cinco genes objetivo principales, reveló su enriquecimiento en vías inmunitarias clave y validó la unión estable a la proteína BaP, elucidando los posibles mecanismos moleculares de la artritis reumatoide (AR) inducida por contaminantes ambientales.

Resumen

Los hidrocarburos aromáticos policíclicos (HAP), contaminantes ambientales ubicuos, se consideran factores ambientales significativos que contribuyen a la patogénesis de la AR. El benzo[a]pireno (BaP), un componente clave de las HAP, puede estar asociado con el inicio de la AR; sin embargo, el mecanismo toxicológico subyacente aún no está completamente esclarecido. En este estudio, abordamos sistemáticamente esta brecha de conocimiento utilizando un enfoque integrado que combina toxicología de redes, aprendizaje automático y acoplamiento molecular. Inicialmente, se realizó un análisis toxicológico en red basado en la estructura molecular de BaP. Mediante la integración y el cribado de información objetivo de múltiples bases de datos, se identificaron finalmente 15 genes objetivo relacionados con la AR, y se construyó su red de interacción. Los análisis de enriquecimiento GO y KEGG revelaron que estos genes estaban significativamente enriquecidos en procesos biológicos como la migración de leucocitos y la transducción de señales celulares inmunitarias, y estaban asociados con las vías de señalización de los receptores NF-κB y T, entre otros. Análisis topológicos posteriores utilizando la base de datos STRING y el software Cytoscape filtraron cinco genes centrales (LCK, ZAP70, ITK, GZMA e ITGAL), cuya importancia fue validada aún más mediante aprendizaje automático. Los resultados del acoplamiento molecular y la simulación de dinámica molecular indicaron que BaP muestra una fuerte afinidad de unión con los productos proteicos de estos genes objetivo, lo que resulta en la formación de complejos conformacionalmente estables. En resumen, este estudio emplea un enfoque computacional integrado para esclarecer los posibles mecanismos por los cuales BaP puede contribuir al desarrollo de la AR, ofreciendo así una base teórica para futuras investigaciones sobre la prevención y tratamiento de la AR asociada a contaminantes ambientales.

Introducción

La AR es una enfermedad autoinmune común caracterizada por la sinopolitis crónica y la formación de panus, que provocan la destrucción progresiva del cartílago y la erosión ósea. Estos cambios patológicos provocan disfunción articular, un mayor riesgo de fracturas patológicas y, en última instancia, discapacidad, perjudicando gravemente la calidad de vida de los^{afectados 1}. Aunque la etiología de la AR no se comprende completamente, su patogénesis se atribuye generalmente a los efectos combinados de factores genéticos (por ejemplo, subtipos de genes HLA-DR4 y HLA-DR1) y ambientales (por ejemplo, tabaquismo, consumo de alcohol, infección por el virus de Epstein-Barr, exposición a la contaminación del aire⁾². Los estudios epidemiológicos han establecido una asociación significativa entre los HAP derivados de la contaminación atmosférica y un mayor riesgo de aparición de AR^.

Los HAP son contaminantes atmosféricos comunes que se originan en la combustión incompleta de sustancias como el carbón, el petróleo, el gas natural y el tabaco, y se consideran importantes mediadores ambientales de la patogénesis de la AR. Pueden unirse al complejo receptor de hidrocarburos arílicos (RSA) en las células inmunitarias, lo que conduce a la exposición de señales de localización nuclear y a la consecuente translocación del complejo ligando–RSA en el núcleo. En el contexto de la patogénesis de la AR, la RHA unida a HAP actúa como un sensor ambiental central que impulsa la desregulación inmunitaria a través de vías interconectadas. Al activarse, la RHA forma heterodímeros con ARNT y modula la expresión del gen CYP, iniciando una cascada^{inflamatoria 4}. Esta activación de la RSA altera simultáneamente el equilibrio entre las células T proinflamatorias y reguladoras: por un lado, activa el eje de señalización de la RCA/Jag1/Notch, potenciando^{la liberación} 5 de citocinas de células Th17. Por otro lado, la RSA se une directamente al promotor de GOT1 para aumentar la expresión de GOT1, lo que induce la hipermetilación del locus FOXP3 y suprime la diferenciación de^{Treg 6}. El desequilibrio resultante entre Th17 y Treg favorece un ambiente proinflamatorio. Además, la activación de la RSA amplifica aún más las respuestas tipo Th2 al aumentar la expresión de CCR8 y citocinas como IL-4 e IL-13, que colectivamente perpetúan la inflamación sinovial y el daño^{tisular 7}. Así, la activación de la RSA sirve como un centro que vincula la exposición ambiental a la HAP con la desregulación de Th17/Treg/Th2, estableciendo un puente mecanicista entre el genotipo, los factores ambientales y la patología de la AR.

En la atmósfera, los HAP existen como mezclas complejas, siendo el BaP un componente clave. En macrófagos humanos, BaP puede inducir la producción de CXCL8 (IL-8) promoviendo la unión de la AHR al promotor CXCL8, induciendo posteriormente la expresión de factores quimiotácticos^{neutrófilos 8}. Además, BaP puede aumentar la expresión de Slug en sinoviocitos similares a fibroblastos (FLS) de pacientes con AR de forma dependiente de la dosis, lo que agrava la progresión de la^artritis. En ratones de tipo salvaje, BaP promueve la activación del osteoclast (OC) mediada por el factor receptor activador del factor nuclear kappaB Ligando (RANKL) induciendo la actividad enzimática CYP1A1, lo que finalmente conduce a la pérdida^{ósea 10}. Sin embargo, el mecanismo preciso subyacente al papel de la toxicidad por BaP en la patogénesis de la AR sigue sin estar claro. Planteamos la hipótesis de que BaP promueve la patogénesis de la AR interactuando directamente con proteínas clave relacionadas con el sistema inmunológico y alterando múltiples vías de señalización implicadas en la activación de las células T, el equilibrio Th17/Treg y la producción de citoquinas inflamatorias, vinculando así la exposición ambiental a BaP con la inflamación sinovial y la destrucción articular. La toxicología de redes es más adecuada que los métodos experimentales tradicionales de vía única para este estudio, ya que el BaP probablemente actúa sobre múltiples objetivos inmunológicos y vías que se intersectan. En comparación con los enfoques experimentales convencionales que normalmente examinan una vía o varios objetivos a la vez, la toxicología de redes permite una visión holística de las interacciones multi-objetivo y los efectos sistémicos, aunque sus predicciones dependen de la base de datos y requieren validación experimental.

La investigación existente sobre el papel de BaP en la AR se limita principalmente a vías individuales o descripciones mecanicistas lineales, careciendo de análisis integrado de las características regulatorias de redes multiobjetivo y multinivel. La investigación existente sobre el papel de BaP en la AR se limita principalmente a vías individuales o descripciones mecanicistas lineales, careciendo de análisis integrado de las características regulatorias de redes multiobjetivo y multinivel. Por ello, se necesitan enfoques holísticos como la toxicología de redes para desentrañar el complejo vínculo entre la exposición a BaP y la patogénesis de la^AR 11,12,13. No obstante, existe una escasez de estudios toxicológicos en red relacionados con enfermedades inducidas por contaminantes ambientales. La novedad de este estudio radica en integrar la toxicología de redes, el aprendizaje automático y el acoplamiento molecular para investigar sistemáticamente la patogénesis de la AR mediada por BaP, en lugar de centrarse en una sola vía o objetivos aislados. Identifica genes clave de núcleo mediante análisis topológico combinado con aprendizaje automático y, por primera vez, proporciona una validación a nivel molecular de los modos de unión y la estabilidad termodinámica entre BaP y los productos de cada gen central. Al identificar sistemáticamente los posibles mecanismos moleculares a través de los cuales BaP puede promover la aparición y el desarrollo de la AR, este estudio computacional pretende proporcionar una base teórica para comprender los desencadenantes ambientales de la AR y desarrollar estrategias terapéuticas específicas. En comparación con los análisis tradicionales de vía única, este enfoque integrado permite una evaluación sistemática de interacciones multi-objetivo, ofreciendo una aplicabilidad más amplia en el estudio de mecanismos complejos de enfermedades ambientales. Sin embargo, cabe señalar que nuestro método prioriza los genes hub de alta confianza mediante la intersección y el análisis topológico, lo que puede excluir inadvertidamente a genes candidatos biológicamente relevantes que no cumplen simultáneamente los umbrales de selección. Estudios futuros podrían explorar estrategias complementarias, como aplicar aprendizaje automático al conjunto de objetivos predichos de unión, integrar otros datos ómicos o realizar validaciones experimentales dirigidas, para confirmar y ampliar aún más nuestros hallazgos.

Protocolo

Declaración ética
Este estudio no involucró directamente a ningún participante humano ni a sujetos animales.

Adquisición de objetivos de BaP
BaP se caracterizaba por integrar datos de múltiples bases de datos. La base de datos PubChem (https://pubchem.ncbi.nlm.nih.gov/) fue consultada usando la palabra clave "Benzo[a]pyrene" para obtener su estructura química y la estructura canónica 2D (cadena SMILES: C1=CC=C2C3=C4C(=CC2=C1)C=CC5=C4C(=CC=C5)C=C3)¹⁴. Los posibles objetivos de BaP se recuperaron de las bases de datos ChEMBL (https://www.ebi.ac.uk/chembl/), SEA (https://sea.bkslab.org/) y PharmMapper (http://lilab-ecust.cn/pharmmapper)^{15, 16 y 17}. Todos los objetivos previstos se restringieron al proteoma de Homo sapiens. La lista completa de objetivos BaP predichos (n = 474) se proporciona en la Tabla Suplementaria S1. El flujo de trabajo analítico completo se representa esquemáticamente en la Figura 1.

Figura 1. Diagrama de flujo del análisis de conjuntos de datos en este artículo, que ilustra el flujo de trabajo general incluyendo adquisición de datos, preprocesamiento, análisis de expresiones diferenciales, construcción de redes y pasos de validación. Por favor, haz clic aquí para ver una versión ampliada de esta figura.

Adquisición de objetivos relacionados con la AR
En este estudio, se adquirieron cinco conjuntos de datos de ARA de la base de datos Omnibus de Expresión Génica (GEO) del NCBI (https://www.ncbi.nlm.nih.gov/gds/) utilizando las palabras clave "Artritis reumatoide" y "Homo sapiens^"18. Basándose en el tamaño del conjunto de datos y el diseño experimental, GSE77298 (RA: 16 muestras; Control: 7 muestras), GSE1919 (RA: 5 muestras; Control: 5 muestras), y GSE55235 (RA: 10 muestras; Control: 10 muestras) formaron el conjunto de entrenamiento para identificar genes expresados diferencialmente (DEGs), mientras que GSE12021 (RA: 24 muestras; Control: 13 muestras) y GSE55457 (RA: 13 muestras; Control: 10 muestras) sirvieron como conjunto de validación. Más detalles sobre estos conjuntos de datos, como plataformas, muestras y series GSE, pueden encontrarse en la Tabla 1.

Los datos se estandarizaron utilizando la herramienta online GEO2R, generando matrices de expresión transformadas en log2 para análisis posteriores. Para eliminar interferencias de diferentes lotes experimentales, se corrigieron sesgos sistemáticos entre conjuntos de datos utilizando la función ComBat del paquete SVA basada en un marco empírico paramétrico de Bayes. Posteriormente, se utilizó el Análisis de Componentes Principales (PCA) para verificar el efecto de corrección, mostrando una mejora significativa en el agrupamiento entre lotes de muestras y, por tanto, la eliminación efectiva de los efectos por lotes. La matriz de datos fusionada y corregida se utilizó para el análisis diferencial posterior.

Serie GSE	Muestras	Andén	Grupo
GSE77298	16 controles RA y 7	GPL570	Cohorte de formación
GSE1919	5 RA y 5 controles	GPL91	Cohorte de formación
GSE55235	10 RA y 10 controles	GPL96	Cohorte de formación
GSE12021	24 controles RA y 13	GPL96	Cohorte de validación
GSE55457	13 controles RA y 10	GPL9	Cohorte de validación

Tabla 1: Resumen de los cinco conjuntos de datos GEO utilizados en este estudio.
La tabla proporciona el número de acceso GEO (serie GSE), la composición muestral (número de pacientes con artritis reumatoide y controles sanos), el identificador de plataforma (GPL) para cada conjunto de datos y la asignación a la cohorte de formación o a la cohorte de validación.

Análisis ponderado de la red de coexpresión génica (WGCNA)
Se utilizó WGCNA para evaluar las características de la red de coexpresións de los DEGs asociados con la^{AR 19}. Basándose en la matriz de expresión corregida por efecto de lotes, se realizó primero el preprocesamiento de datos: se eliminaron genes de baja varianza con desviación estándar inferior a 0,5, mientras que la calidad de la muestra y el gen se evaluaron utilizando una función para evaluar buenas muestras y genes. Posteriormente, se aplicó agrupamiento jerárquico para identificar y eliminar muestras atípicas. Para construir una red ponderada de coexpresiones, se empleó una función para la evaluación sistemática de valores de potencia de umbral suave para evaluar sistemáticamente valores de potencia de umbral suave que van de 1 a 20. Se seleccionó Power = 12 como umbral blando óptimo (topología libre de escala, índice de ajuste R2 = 0,90), asegurando que la topología de la red cumpliera con un criterio libre de escala. A partir de este valor de potencia, se construyó una matriz de adyacencia y se calculó la matriz de solapamiento topológico (TOM). Los genes se agrupaban jerárquicamente y se utilizaba un algoritmo dinámico de corte de árboles para identificar los módulos génicos iniciales. Posteriormente, módulos similares se fusionaron mediante la agrupación de los eigengenes de módulos, dando lugar a una red robusta de módulos génicos. Todos los análisis se realizaron con un paquete R dedicado para el análisis ponderado de redes de coexpresión, con el fin de asegurar la fiabilidad y reproducibilidad de la construcción de la red. Se realizó un análisis de la intersección entre los genes central de DEGs/WGCNA y los objetivos predichos de BaP para identificar objetivos centrales de BaP asociados a la patogénesis de AR, que se visualizaron mediante el software de diagrama de Venn.

Identificación de objetivos asociados a BaP asociados con la patogénesis de la AR
El análisis de intersección se realizó utilizando un paquete R para diagramas de Venn para identificar objetivos de BaP que se solapan con la patogénesis de AR. Estos se importaron a la base de datos STRING para construir una red de interacción proteína-proteína (PPI), con la especie configurada en "Homo sapiens" y la puntuación de confianza en interacción en > 0,7 para asegurar una alta fiabilidad de la red²⁰. Este umbral se seleccionó porque corresponde a un nivel de "alta confianza" en la base de datos STRING, que equilibra la retención de interacciones biológicamente relevantes mientras minimiza los falsos positivos típicamente asociados a puntuaciones de confianza más bajas. Un corte de > 0,7 ha sido ampliamente adoptado en estudios de toxicología en redes para priorizar asociaciones proteicas robustas y reproducibles. El archivo TSV resultante se descargaba de la base de datos de interacción proteína-proteína (STRING) e importaba al software de visualización de red (Cytoscape) para la visualización de la red. Las proteínas centrales de la red se identificaron basándose en los resultados de clasificación generados por el algoritmo Degree en el plugin CytoHubba y se utilizaron para análisis posteriores.

Análisis de enriquecimiento KEGG y GO
Las abreviaturas de los genes asociados tanto con la modulación BaP como con la patogénesis de la AR se convirtieron en identificadores de Entrez usando el "org. Paquete de anotación Hs.eg.db" en R. Posteriormente, se realizó un análisis de enriquecimiento de vías KEGG utilizando la herramienta clusterProfer, con el umbral de significación fijado en 0,05. Mientras tanto, la anotación funcional GO abarcaba las tres principales categorías de GO: Proceso Biológico (BP), Componente Celular (CC) y Función Molecular (MF), y se realizaba utilizando la función enrichGO, con cortes tanto de valores P como de valor q establecidos en 0,05. Cabe señalar que no se aplicó corrección mediante múltiples pruebas, ya que el objetivo principal de este análisis exploratorio era maximizar el descubrimiento de posibles vías biológicas relevantes y términos funcionales, generando así un conjunto más amplio de hipótesis comprobables para futuras validaciones experimentales. Finalmente, los resultados del análisis de enriquecimiento se mostraron gráficamente usando las funciones de gráfico de barras y puntos del paquete enrichplot.

Validación basada en aprendizaje automático de genes centrales
Para evaluar la capacidad predictiva de los genes centrales asociados con BaP y AR, y mantener la transparencia del modelo, implementamos un flujo de trabajo sistemático de aprendizaje automático. Utilizando los perfiles de expresión de los genes centrales seleccionados, se construyeron modelos predictivos con 11 algoritmos de aprendizaje automático distintos: regresión de lazo (LR), Máquina de Vectores de Soporte (SVM), Bosque Aleatorio (RF), glmBoost, Modelo Lineal Generalizado Escalón (GLM), regresión de crestas, red elástica (Enet), Máquina de Aumento de Gradiente (GBM), Análisis Discriminante Lineal (LDA), Aumento de Gradiente EXtreme (XGBoost) y Bayes naïve. Los hiperparámetros se optimizaron mediante una validación cruzada de cinco vías, utilizando muestreo estratificado para dividir los datos en conjuntos de entrenamiento y validación interna. Se utilizó una semilla aleatoria fija (set.seed(123)) a lo largo del flujo de trabajo de aprendizaje automático para garantizar la reproducibilidad de la división de datos, los pliegues de validación cruzada y el entrenamiento de modelos. Los hiperparámetros clave para cada algoritmo se proporcionan en la Tabla Suplementaria S2. El rendimiento del modelo se evaluó utilizando múltiples métricas, incluyendo área bajo la curva (AUC), precisión y puntuación F1. Para abordar las limitaciones inherentes a los enfoques de modelo único, aplicamos una estrategia de ensamble apilado que integró predicciones de los modelos base con mejor rendimiento. Reconociendo la naturaleza de "caja negra" de muchos modelos de aprendizaje automático, empleamos el algoritmo SHapley Aditivive ExPlanations (SHAP) para cuantificar la contribución de cada gen a las predicciones. La magnitud y dirección de los valores de SHAP se utilizaron para interpretar la importancia génica en las decisiones de clasificación, mejorando así la interpretabilidad de los resultados del modelo.

Acoplamiento molecular de BaP con objetivos centrales
Para investigar las características de unión entre BaP y los productos génicos centrales, se realizaron simulaciones de acoplamiento molecular. La estructura tridimensional de BaP (ligando) se obtuvo en formato SDF a partir de la base de datos PubChem. Las estructuras proteicas correspondientes a los objetivos centrales se recuperaron del Banco de Datos de Proteínas (https://www.rcsb.org/) RCSB en formato PDB, seleccionadas según sus identificadores UniProt, con preferencia por estructuras que contienen ligandos cocristalizados o coordenadas de alta resolución. Antes del acoplamiento, la preparación de proteínas se realizó usando PyMol, durante la cual se eliminaron moléculas de agua, ligandos cocristalizados y componentes no proteicos como iones para evitar^{interferencias 21}. Para proteínas con ligandos cocristalizados en sus estructuras originales de PDB, el centro del sitio activo se definió usando las coordenadas atómicas del ligando unido. Para proteínas sin ligandos cocristalizados, el centro del sitio activo se determinó en función de las coordenadas de residuos clave que la literatura reporta como críticos para la actividad catalítica o la unión a inhibidores. La cuadrícula de acoplamiento estaba centrada en las coordenadas definidas del sitio activo, con una caja cúbica de dimensiones de 25 × 25 × 25 Å aplicada a cada objetivo. Este tamaño estándar de caja de 25 Å garantiza una cobertura completa de cada sitio activo con margen suficiente para el muestreo de ligandos, evitando al mismo tiempo un coste computacional excesivo. Todos los cálculos de acoplamiento se ejecutaban con AutoDock Vina (versión 1.2.5). Se seleccionó la conformación que presentaba la puntuación de Vina más favorable como modo de enlace representativo, y se registró la energía de enlace correspondiente. Se generaron poses de unión tridimensional usando PyMol (versión 2.5.7), y se produjeron diagramas de interacción bidimensional usando Discovery Studio (versión 2021) para visualizar interacciones clave, incluyendo enlaces de hidrógeno y contactos hidrofóbicos.

Simulación de dinámica molecular
Se realizaron simulaciones de dinámica molecular con Gromacs 2025.3, utilizando los complejos derivados del acoplamiento como estructuras iniciales. Los átomos de la proteína se modelaron con el campo de fuerza AMBER14SB, y las moléculas de agua se representaron usando el modelo TIP3P. Cada complejo proteína-ligando se solvía en una caja cúbica de agua, con una distancia mínima de 1 nm entre la superficie de la proteína y el límite de la caja. Se añadieron iones de sodio o cloruro según fuera necesario para lograr la electroneutralidad del sistema. Se realizó una minimización inicial de energía utilizando una combinación de algoritmos de descenso más pronunciado y gradiente conjugado, cada uno con hasta 10.000 escalones. Las interacciones electrostáticas de largo alcance se calcularon mediante el método Particle-Mesh Ewald (PME), mientras que se aplicó una distancia de corte de 1,0 nm tanto a van der Waals como a interacciones electrostáticas de corto alcance. Tras la minimización de energía, los sistemas se equilibraron gradualmente bajo condiciones NVT (volumen y temperatura constantes) y NPT (presión y temperatura constantes). A continuación, se realizaron tiradas de producción de 100 ns bajo temperatura y presión constantes, con un paso de tiempo de 0,002 ps (2 fs) y un total de 50.000.000 de pasos. Cada simulación se realizó una vez (sin replicaciones), ya que el objetivo principal era evaluar la estabilidad de los complejos de unión bajo condiciones estándar. La temperatura se mantenía mediante el termostato V-rescale y la presión se controlaba con el barostato de Parrinello–Rahman. Durante toda la simulación, se aplicó consistentemente un corte de 1,0 nm para las interacciones no enlazadas. Para evaluar la estabilidad y flexibilidad estructural, calculamos la desviación cuadrática media (RMSD) de las posiciones atómicas, la fluctuación cuadrática media (RMSF) por residuo, el radio de giro (Rg) como medida de la compacidad estructural y la superficie accesible al disolvente (SASA). Todos los gráficos se generaron usando QtGrace.

Resultados

Adquisición de objetivos de BaP
Los datos sobre la estructura molecular de BaP se obtuvieron de la base de datos PubChem (Figura 2A). Los posibles objetivos biológicos de BaP se predijeron sistemáticamente integrando información de tres bases de datos complementarias—ChEMBL, PharmMapper y SEA—resultando en la identificación de 474 objetivos potenciales (Figura 2B).

Discusión

La AR es una enfermedad autoinmune compleja que resulta de la interacción entre la susceptibilidad genética y los factores ambientales. Entre numerosos factores de riesgo ambientales, los HAP, uno de los contaminantes atmosféricos más comunes, se consideran un vínculo importante que conecta la exposición ambiental con el inicio de la AR. Estudios anteriores revelaron preliminarmente que los HAP pueden influir en el equilibrio de la diferenciación de células inmunitarias a través de la ví...

Materiales

Lista de materiales utilizados en este artículo

Nombre	Empresa	Número de catálogo	Comentarios
AutoDock Vina	https://vina.scripps.edu	1.2.5 (SCR_011958)	Software de acoplamiento molecular
clusterProfiler (paquete R)	https://bioconductor.org/packages/clusterProfiler	4.10.0 (SCR_016884)	Paquete R para análisis de enriquecimiento
CytoHubba (plugin de Cytoscape)	https://apps.cytoscape.org/apps/cytohubba	0.1	Plugin para la identificación de genes hub (algoritmo de grado)
Cytoscape	https://cytoscape.org	3.10.1 (SCR_003032)	Software de visualización de redes
Discovery Studio	Dassault Systè mes BIOVIA	2021	Software para la generación de diagramas de interacción 2D
enrichplot (paquete R)	https://bioconductor.org/packages/enrichplot	1.22.0 (SCR_021165)	Paquete R para la visualización de resultados de enriquecimiento
GROMACS	https://www.gromacs.org	2025.3 (SCR_014565)	Software de simulación de dinámica molecular
limma (paquete R)	https://bioconductor.org/packages/limma	3.58.1 (SCR_010943)	Paquete R para análisis de expresiones diferenciales
org. Hs.eg.db (paquete R)	https://bioconductor.org/packages/org. Hs.eg.db	3.18.0 (SCR_006442)	Paquete de anotación R para identificadores genéticos humanos
PyMol	Schrö dinger, Inc	2.5.7 (SCR_000305)	Software de visualización molecular
QtGrace	https://sourceforge.net/projects/grace/	0.2.6	Herramienta de trazado para el análisis de trayectorias
R (entorno de programación)	https://www.r-project.org	4.3.1 (SCR_001905)	Software de computación estadística
Base de datos STRING	https://string-db.org	12 (SCR_005223)	Base de datos de interacción proteína-proteína
venn (paquete R)	https://cran.r-project.org/package=venn	1.11	Paquete R para la generación de diagramas de Venn
WGCNA (paquete R)	https://cran.r-project.org/package=WGCNA	1,72 (SCR_003302)	Paquete dedicado R para análisis de redes de coexpresiones ponderadas