Divergencia de la microbiota radicular en diferentes hábitats basada en redes de correlación ponderadas

Resumen

Se aplicó el análisis de redes para evaluar la asociación de varias comunidades microbianas ecológicas, como el suelo, el agua y la rizosfera. Aquí se presenta un protocolo sobre cómo usar el algoritmo WGCNA para analizar diferentes redes de co-ocurrencia que pueden ocurrir en las comunidades microbianas debido a diferentes entornos ecológicos.

Resumen

El microbioma radicular juega un papel importante en el crecimiento de las plantas y la adaptación ambiental. El análisis de redes es una herramienta importante para estudiar las comunidades, que pueden explorar eficazmente la relación de interacción o el modelo de co-ocurrencia de diferentes especies microbianas en diferentes entornos. El propósito de este manuscrito es proporcionar detalles sobre cómo utilizar el algoritmo de red de correlación ponderada para analizar diferentes redes de co-ocurrencia que pueden ocurrir en comunidades microbianas debido a diferentes entornos ecológicos. Todo el análisis del experimento se realiza en el paquete WGCNA. WGCNA es un paquete R para el análisis de redes de correlación ponderada. Los datos experimentales utilizados para demostrar estos métodos fueron datos de la comunidad microbiana de la base de datos NCBI (Centro Nacional de Información Biotecnológica) para tres nichos del sistema radicular del arroz (Oryza sativa). Utilizamos el algoritmo de red de correlación ponderada para construir redes de coabundancia de comunidad microbiana en cada uno de los tres nichos. Luego, se identificaron redes diferenciales de coabundancia entre el suelo de endosfera, rizóplano y rizosfera. Además, los géneros centrales en la red se obtuvieron mediante el paquete "WGCNA", que desempeña un importante papel regulado en las funciones de red. Estos métodos permiten a los investigadores analizar la respuesta de la red microbiana a la perturbación ambiental y verificar diferentes teorías de respuesta ecológica microbiana. Los resultados de estos métodos muestran que las importantes redes microbianas diferenciales identificadas en la endosfera, el rizóplano y la rizosfera del suelo del arroz.

Introducción

La investigación del microbioma tiene implicaciones importantes para la comprensión y manipulación de los procesos de los^{ecosistemas 1,2}. Las poblaciones microbianas están interconectadas por redes ecológicas que interactúan, cuyas características pueden afectar la respuesta de los microorganismos a los cambios ambientales^3,4. Además, las propiedades de estas redes afectan a la estabilidad de las comunidades microbianas, y están estrechamente asociadas con la función del suelo⁵. El análisis de redes de correlación de genes ponderados ahora se ha aplicado ampliamente para la investigación sobre la relación entre los genes y las comunidades microbianas⁶. Estudios previos se han centrado principalmente en las asociaciones entre redes de diferentes genes o poblaciones y el mundo exterior⁷. Sin embargo, las diferencias en las redes de correlación formadas por poblaciones microbianas bajo diferentes condiciones ambientales apenas se han investigado. El propósito de la investigación presentada en este documento es proporcionar información y detalles sobre la rápida implementación del algoritmo WGCNA para construir una red de co-ocurrencia de muestras de microbioma recolectadas bajo diferentes condiciones ambientales. Sobre la base de los resultados del análisis, evaluamos la composición y las diferencias de la población y discutimos más a fondo la relación entre las diferentes poblaciones microbianas. Se aplicó el siguiente flujo básico del algoritmo de red de correlación ponderada^8. En primer lugar, era necesario construir una matriz de similitud calculando el coeficiente de correlación de Pearson entre los perfiles de expresión de unidades taxonómicas operativas (OTU). Luego, los parámetros de las funciones de adyacencia (las funciones de potencia o de adyacencia sigmoide) se adoptaron con un criterio de topología libre de escala, la matriz de similitud se transformó en una matriz de adyacencia y cada red de co-ocurrencia correspondió a una matriz de adyacencia. Utilizamos la agrupación jerárquica de enlace promedio junto con la disimilitud basada en TOM para agrupar OTU con perfiles de expresión coherentes en módulos. Además, calculamos la relación entre la estadística conservadora y los módulos de análisis de parámetros relacionados, identificando finalmente el hub OTU en el módulo. Estos métodos son particularmente adecuados para el análisis de las diferencias en las estructuras de red entre varias poblaciones microbianas en condiciones ambientales divergentes. En este manuscrito, hemos descrito en detalle el método de desarrollo de redes de coexpresión, el análisis de las diferencias entre los módulos, y hemos proporcionado una breve visión general de los pasos en el procedimiento aplicado para obtener las especies centrales en diferentes redes de módulos.

Protocolo

1. Descarga de datos

1. Descargue los datos de la adhesión PRJNA386367 de la base de datos ncbi. A partir de los datos de la adhesión PRJNA386367, seleccione los datos del microbioma de rizosfera, rizóplano y endosfera de plantas de arroz cultivadas durante 14 semanas en un campo de arroz sumergido en Arbuckle, California, en 2014.
   NOTA: Los datos del microbioma de rizosfera, rizóplano y endosfera fueron presentados por la tabla de OTU en la adhesión PRJNA386367.

2. Determinación óptima del valor de potencia

NOTA: El paquete WGCNA contiene todos los siguientes parámetros funcionales. WGCNA es un paquete R para el análisis de redes de correlación ponderada. Las líneas de comandos clave hacen referencia al Suplemento S1.

1. En el entorno del lenguaje R, abra el software Rstudio e instale el paquete WGCNA.
2. Cargue los datos y utilice la función goodSamplesGenes para comprobar la exactitud de los datos. Ejecute las líneas de comandos:
   "gsg = goodSamplesGenes(datExpr0, verbose = 3)
   gsg$allOK "
   Haga clic en Ejecutar.
3. Compruebe si hay valores atípicos y almacene muestras que cumplan con los requisitos. Cuando el resultado de la comprobación sea TRUE, continúe con el siguiente paso. Guarde el resultado.
4. Utilice la función PickSoftThreshold para calcular el índice R² sin escala de los dos grupos de datos con diferentes valores de potencia. Ejecute la línea de comandos:
   "sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)"
   Haga clic en Ejecutar.
5. Visualizar los resultados (Figura 1). Ejecute la línea de comandos:
   "plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
   xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
   main = paste("ES_Scale independencia"));
   text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
   labels=powers,cex=cex1,col="red");
   abline(h=0.9,col="rojo")
   plot(sft$fitIndices[,1], sft$fitIndices[,5],
   xlab="Umbral suave (potencia)",ylab="Conectividad media", type="n",
   main = paste("conectividad ES_Mean"))
   text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")"
   Haga clic en Ejecutar.
   NOTA: La premisa del algoritmo de red de correlación ponderada es que la estructura de red de coexpresión establecida se ajusta a los estándares del criterio de topología libre de escala, aumentando su robustez. Un índice sin escala más cercano a 1 indica una estructura de red que está más cerca de la red sin escala.
6. Seleccione el valor de potencia cuando el índice sin escala R² cuadreó más de 0,9 y continúe con el siguiente paso del análisis.
   NOTA: Cuando el índice sin escala está cerca de 1, la estructura de red está más cerca de la red sin escala. Al analizar dos o más redes, es necesario elegir hacer que cada red se acerque al valor de potencia de la red libre de escala para satisfacer la comparabilidad entre las redes co-expresadas.

3. Construcción de una red de coexpresión e identificación de módulos

NOTA: Sobre la base del valor de potencia calculado anteriormente, se construye la red de co-ocurrencia. Las líneas de comandos clave hacen referencia al Suplemento S2.

1. Utilice la función de adyacencia del paquete WGCNA para agregar parámetros firmados para la construcción de una red de co-ocurrencia simbólica. Ejecute la línea de comandos:
   "adyacencia = adyacencia(datExpr0, poder = softPower)"
   Haga clic en Ejecutar.
2. Aplique la función TOM-similarity para desarrollar una red topológica superpuesta y calcule la red de disimilitud. Ejecute la línea de comandos:
   "TOM = TOMsimilarity(adyacencia);
   dissTOM = 1-TOM"
   Haga clic en Ejecutar.
   Nota : el parámetro firmado se agregó para establecer el tipo de red de superposición de topología.
3. Utilice la función hclust para seleccionar el método de agrupación jerárquica de vinculación media para la agrupación en clústeres jerárquicas. Ejecute la línea de comandos:
   "geneTree = hclust(as.dist(dissTOM), method = "average");"
   Haga clic en Ejecutar.
4. Utilice la función cutreeDynamic para realizar el corte dinámico por rama y establezca el parámetro minClusterSize en 30. Obtener el resultado de reconocimiento del módulo. Ejecute la línea de comandos:
   "dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize);"
   Haga clic en Ejecutar.
   NOTA: El tamaño mínimo del módulo no puede ser inferior a 30.
5. Calcule el eigen del módulo de cada módulo otus mediante la función moduleEigengenes. Ejecute la línea de comandos:
   "MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
   MEs = MEList$eigengenes"
   Haga clic en Ejecutar.
   NOTA: El eigen del módulo representaba el nivel de expresión OTU general en el módulo. No era una OTU específica, sino el primer componente principal de cada clúster obtenido por descomposición de valores de red singulares.
6. Realice la función de clúster basada en el coeficiente de correlación del módulo eigen. Utilice la función mergeCloseModules para combinar los módulos con un valor inferior a 0,25. Ejecute la línea de comandos:
   "merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)"
   Haga clic en Ejecutar.
7. Finalmente, utilice la función plotDendroAndColors para la visualización para obtener el diagrama de visualización de asignación de módulos de cada red de coexpresión (Figura 2). Utilice la función de tabla para extraer la atribución del módulo correspondiente de cada OTin la tabla de asignación del módulo. Ejecute la línea de comandos:
   "plotDendroAndColors(geneTree, mergedColors, "Merged dynamic",dendroLabels = FALSE,
   hang = 0.03,addGuide = TRUE, guideHang = 0.05,
   main = "ES_Gene dendrograma y colores de módulo")"
   Haga clic en Ejecutar.
   NOTA: En el diagrama de asignación de módulos de la red de co-expresión, diferentes colores representan diferentes módulos, y el gris representa OTU que no se pueden clasificar en ningún módulo. Un mayor número de OTU en el módulo gris indica que la calidad de preprocesamiento en etapa inicial de la matriz de expresión es deficiente.

4. Comparación de módulos

NOTA: Este método se puede utilizar para comparar los módulos de red de dos comunidades microbianas ecológicas. En este artículo, compare las diferencias de los módulos de red microbiana entre endosfera y rizóplano, endosfera y rizosfera, rizosfera y rizóplano.

1. Prueba de conservación
   1. Cargue los parámetros y resultados de los dos conjuntos de datos guardados en los pasos anteriores.
   2. Establezca el resultado de asignación del módulo de red de un grupo de datos microbianos como grupo de referencia, mientras que el otro grupo como grupo de prueba.
   3. Utilice la función modulePreservation para calcular los valores de los parámetros estadísticos de conservadurismo Z_summary y medianRank. Ejecute la línea de comandos:
      "system.time({mp=modulePreservation(multiExpr,
      multiColor,referenceNetworks=1,
      nPermutation=100, randomSeed=1,quickCor=0,verbose=3)})"
      Haga clic en Ejecutar.
      NOTA: Este resultado puede cuantificar el conservadismo entre módulos. Z_summary>10 indica que dos módulos están altamente conservados, mientras que Z_summary<2 denota módulos no conservados. medianRank expresa la preservación relativa del módulo evaluado por ranking. Los valores de medianRank más altos denotan módulos no conservados. (Las líneas de comandos clave se refieren al Suplemento S3.)
   4. Utilice la función de trazado para visualizar los resultados (Figura 3). Obtenga los parámetros Z_summary y medianRank (Tabla 1).
      NOTA: Los módulos de red que satisfacen tanto el valor Z_summary inferior a 2 como el valor medio Rank en la parte superior, es el módulo más no conservado en las dos comunidades microbianas ecológicas.
   5. Basado en los resultados de los dos parámetros estadísticos antes mencionados para identificar el módulo con el módulo más altamente no conservado de las dos redes.
2. Análisis de correlación de la pertenencia al módulo
   1. Establezca los resultados de asignación de módulos de las dos redes que se establecieron como la referencia y el grupo de prueba, respectivamente.
      NOTA: La configuración debe ser la misma que la prueba de preservación.
   2. Utilice la función corPvalueStudent para extraer el valor kME (pertenencia a módulos) de cada OTU en varios módulos candidatos.
      Ejecute la línea de comandos:
      "Pvalue = as.data.frame(corPvalueStudent(as.matrix)
      (ModuleMembership), Samples))"
      Haga clic en Ejecutar.
      NOTA: kME significa el grado de pertenencia al módulo. ME significa módulo eigen, que representa el nivel general de expresión OTU en el módulo. kME es el coeficiente de correlación entre cada OTU y el ME. Cuantificar la importancia de OTU en la red por el valor kME de OTU. (Las líneas de comandos clave se refieren al Suplemento S4.)
   3. A continuación, utilice la función verboseScatterplot para calcular el coeficiente de correlación del valor kME de las OTU correspondientes en las dos redes y dibuje el diagrama de análisis de correlación(Figura 4).
      Ejecute la línea de comandos:
      "verboseScatterplot(abs(TModuleMembership
      [TmoduleGenes, Tcolumn]),
      abs(NModuleMembership[NmoduleGenes, Ncolumn]),
      xlab = paste("kME in", "ES"),
      ylab = paste("kME in", "RP"),
      main = pasta("lightyellow"),
      cex.main = 1.7, cex.lab = 1.6, cex.axis = 1.6, col = modulecolor)"
      Haga clic en Ejecutar.
   4. Seleccione el módulo con el coeficiente de correlación más pequeño del valor kME de la OTU de las dos redes. Considere que este módulo tiene la mayor diferencia de las dos redes.

5. Análisis del módulo de red diferencial microbiana

1. Obtener datos de los filos de bacterias dominantes mediante el análisis estadístico del conjunto de secuencias OTU del módulo con mayor diferencia.
   NOTA: El conjunto de secuencias OTU del módulo con la mayor diferencia se suma a la taxonomía de phyla. Los filos de bacterias dominantes representaron más del 10%.
2. A continuación, utilice la función exportNetworkToCytoscape para obtener el archivo que contiene la información de relación de interacción de la OTU en el módulo diferencial más grande.
   Ejecute la línea de comandos:
   "cyt = exportNetworkToCytoscape(modTOM,
   edgeFile = paste("NEW-ES_CytoscapeInput-edges-", modules , ".txt", sep=""),
   nodeFile = paste("NEW-ES_CytoscapeInput-nodes-", modules, ".txt", sep=""),
   ponderado = TRUE,umbral = 0,5, nodeNames = modProbes,
   altNodeNames = modGenes, nodeAttr = moduleColors[inModule])"
   Haga clic en Ejecutar.
3. Importe el archivo a Cytoscape. Establezca el umbral en 0,5 y ajuste otros parámetros según sea necesario.
4. Construir una red de co-ocurrencia de microorganismos diferenciales (Figura 5).
5. Obtuve la información del género núcleo que tiene el rol regulador más importante en la red.
   NOTA: De acuerdo con el valor kME de OUT, se puede definir el género central.
6. Finalmente, se evaluaron las funciones del género núcleo y se analizó su influencia en toda la red de diferencias.

Resultados

Los resultados representativos de este artículo se descargaron de los datos del microbioma de raíz de arroz Abaker de California 2014 en la base de datos NCBI (PRJNA386367)⁹. Los datos incluyen muestras de rizosfera, rizóplano y microbioma de endosfera de plantas de arroz cultivadas durante 14 semanas en un campo de arroz sumergido. Utilizamos el algoritmo WGCNA para seleccionar el valor de potencia que satisfacía las tres redes que estaban cerca de la red libre de escala (

Discusión

Las redes de correlación se han utilizado cada vez más en aplicaciones bioinformáticas. WGCNA es un método de biología de sistemas para el análisis descriptivo de las relaciones entre varios elementos de un sistema biológico¹². El paquete de software R se utilizó en trabajos anteriores sobre WGCNA^13,14,15. El paquete incluye funciones para la construcción de redes, detección de módulos, cálcul...

Materiales

Name	Company	Catalog Number	Comments
R	The University of Auckland	version 4.0.2	R is a free software environment for statistical computing and graphics. It compiles and runs on a wide variety of UNIX platforms, Windows and MacOS.
RStdio	JJ Allaire	version 1.4.1103	The RStudio IDE is a set of integrated tools designed to help you be more productive with R and Python.
Cytoscape		version 3.7.1	Cytoscape is an open source software platform for visualizing complex networks and integrating these with any type of attribute data.
NCBI database			The National Center for Biotechnology Information advances science and health by providing access to biomedical and genomic information.