$$\rightleftharpoonup{xx}$$
$$\longleftharp{xx}$$,
$$\longrightharp{xx}$$,
1. Download dei dati
- Scarica i dati dell'adesione PRJNA386367 dalla banca dati NCBI. Dai dati dell'adesione PRJNA386367, selezionare i dati del microbioma della rizosfera, del rizopiano e dell'endosfera dalle piante di riso coltivate per 14 settimane in un campo di riso sommerso ad Arbuckle, in California, nel 2014.
NOTA: I dati sul microbioma della rizosfera, del rizopiano e dell'endosfera sono stati presentati dalla tabella delle OTU nell'adesione PRJNA386367.
2. Determinazione ottimale del valore di potenza
NOTA: il pacchetto WGCNA contiene tutti i seguenti parametri funzionali. WGCNA è un pacchetto R per l'analisi ponderata della rete di correlazione. Le righe di comando chiave si riferiscono al Supplemento S1.
- Nell'ambiente in linguaggio R, aprire il software Rstudio e installare il pacchetto WGCNA.
- Caricare i dati e utilizzare la funzione goodSamplesGenes per verificare la correttezza dei dati. Eseguire le righe di comando:
"gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK "
Fare clic su Esegui.
- Verificare la disponibilità di valori anomali e archiviare campioni che soddisfino i requisiti. Quando il risultato del controllo è TRUE, passare al passaggio successivo. Salvare il risultato.
- Utilizzare la funzione PickSoftThreshold per calcolare l'indice senza scala R2 dei due gruppi di dati con valori di potenza diversi. Eseguire la riga di comando:
"sft = pickSoftThreshold(datExpr0, powerVector = potenze, verbose = 5)"
Fare clic su Esegui.
- Visualizzare i risultati (Figura 1). Eseguire la riga di comando:
"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("indipendenza ES_Scale"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
abline(h=0.9;col="rosso")
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = incolla("connettività ES_Mean"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")"
Fare clic su Esegui.
NOTA: la premessa dell'algoritmo di rete di correlazione ponderata è che la struttura di rete di coespressione stabilita è conforme agli standard del criterio della topologia senza scala, aumentandone la robustezza. Un indice senza scalabilità più vicino a 1 indica una struttura di rete più vicina alla rete senza scalabilità.
- Selezionare il valore di potenza quando l'indice senza scala R2 al quadrato è maggiore di 0,9 e procedere alla fase successiva dell'analisi.
NOTA: quando l'indice senza scalabilità è vicino a 1, la struttura di rete è più vicina alla rete senza scalabilità. Quando si analizzano due o più reti, è necessario scegliere di rendere ogni rete vicina al valore di potenza della rete senza scala per soddisfare la comparabilità tra le reti co-espresse.
3. Costruzione di una rete di co-espressione e identificazione del modulo
NOTA: in base al valore di potenza calcolato sopra, viene costruita la rete di co-occorrenza. Le righe di comando chiave si riferiscono al Supplemento S2.
- Utilizzare la funzione di adiacenza nel pacchetto WGCNA per aggiungere parametri firmati per la costruzione di una rete di co-occorrenza simbolica. Eseguire la riga di comando:
"adiacenza = adiacenza(datExpr0, potenza = softPower)"
Fare clic su Esegui.
- Applicare la funzione di somiglianza TOM per sviluppare una rete di sovrapposizione topologica e calcolare la rete di dissomiglianza. Eseguire la riga di comando:
"TOM = TOMsimilarity(adijacency);
dissTOM = 1-TOM"
Fare clic su Esegui.
Nota : il parametro signed è stato aggiunto per impostare il tipo di rete di sovrapposizione della topologia.
- Utilizzare la funzione hclust per selezionare il metodo di clustering gerarchico di collegamento medio per il clustering gerarchico. Eseguire la riga di comando:
"geneTree = hclust(as.dist(dissTOM), metodo = "media");"
Fare clic su Esegui.
- Utilizzate la funzione cutreeDynamic per eseguire il taglio dinamico dei rami e impostate il parametro minClusterSize su 30. Ottenere il risultato del riconoscimento del modulo. Eseguire la riga di comando:
"dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize);"
Fare clic su Esegui.
NOTA: la dimensione minima del modulo non può essere inferiore a 30.
- Calcolare l'auto del modulo di ciascun modulo OTA mediante la funzione moduleEigenes. Eseguire la riga di comando:
"MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
MEs = MEList$eigenes"
Fare clic su Esegui.
NOTA: l'auto del modulo rappresentava il livello di espressione OTU complessivo nel modulo. Non era una OTU specifica, ma il primo componente principale di ogni cluster ottenuto dalla scomposizione del valore di rete singolare.
- Eseguire la funzione cluster in base al coefficiente di correlazione del modulo eigen. Utilizzare la funzione mergeCloseModules per unire i moduli con un valore inferiore a 0,25. Eseguire la riga di comando:
"merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)"
Fare clic su Esegui.
- Infine, utilizzare la funzione plotDendroAndColors per la visualizzazione per ottenere il diagramma di visualizzazione dell'assegnazione del modulo di ogni rete di coespressione (Figura 2). Utilizzare la funzione table per estrarre l'attribuzione del modulo corrispondente di ogni OTin la tabella di assegnazione del modulo. Eseguire la riga di comando:
"plotDendroAndColors(geneTree, mergedColors, "Merged dynamic",dendroLabels = FALSE,
hang = 0,03,addGuide = TRUE, guideHang = 0,05,
main = "ES_Gene dendrogramma e colori del modulo")"
Fare clic su Esegui.
NOTA: nel diagramma di assegnazione dei moduli della rete co-espressiva, colori diversi rappresentano moduli diversi e il grigio rappresenta le OTA che non possono essere classificate in alcun modulo. Un numero maggiore di OTA nel modulo grigio indica che la qualità di pre-elaborazione in fase iniziale della matrice di espressione è scarsa.
4. Confronto dei moduli
NOTA: Questo metodo può essere utilizzato per confrontare i moduli di rete di due comunità microbiche ecologiche. In questo articolo, confronta le differenze dei moduli di rete microbica tra endosfera e rizopiano, endosfera e rizosfera, rizosfera e rizopiano.
- Test di conservazione
- Caricare i parametri e i risultati dei due set di dati salvati nei passaggi precedenti.
- Impostare il risultato dell'assegnazione del modulo di rete di un gruppo di dati microbici come gruppo di riferimento, mentre l'altro gruppo come gruppo di test.
- Utilizzare la funzione modulePreservation per calcolare i valori dei parametri statistici di conservativezza Z_summary e medianRank. Eseguire la riga di comando:
"system.time({mp=modulePreservation(multiExpr,
multiColor,referenceNetworks=1,
nPermutation=100, randomSeed=1,quickCor=0,verbose=3)})"
Fare clic su Esegui.
NOTA: questo risultato può quantificare la conservatività tra i moduli. Z_summary>10 indica che due moduli sono altamente conservati, mentre Z_summary<2 indica moduli non conservati. medianRank esprime la conservazione relativa del modulo valutato per classifica. Valori medianRank più elevati indicano moduli non conservati. (Le righe di comando dei tasti si riferiscono al supplemento S3.)
- Utilizzare la funzione plot per visualizzare i risultati (Figura 3). Ottenere i parametri Z_summary e medianRank (Tabella 1).
NOTA: I moduli di rete che soddisfano sia il valore Z_summary inferiore a 2 che il valore Rank mediano in alto, è il modulo più altamente non conservato nelle due comunità microbiche ecologiche.
- Sulla base dei risultati dei suddetti due parametri statistici per identificare il modulo con il modulo più altamente non conservato delle due reti.
- Analisi di correlazione dell'appartenenza al modulo
- Impostare i risultati di assegnazione del modulo delle due reti sono stati impostati rispettivamente come riferimento e gruppo di test.
NOTA: le impostazioni devono essere uguali a test di conservazione.
- Utilizzare la funzione corPvalueStudent per estrarre il valore kME (appartenenza al modulo) di ogni OTU in diversi moduli candidati.
Eseguire la riga di comando:
"Pvalue = as.data.frame(corPvalueStudent(as.matrix
(ModuleMembership), Samples))"
Fare clic su Esegui.
NOTA: kME sta per il grado di appartenenza al modulo. ME sta per modulo eigen, che rappresenta il livello complessivo di espressione OTU nel modulo. kME è il coefficiente di correlazione tra ogni OTU e il ME. Quantificare l'importanza di OTU nella rete in modo dal valore kME di OTU. (Le righe di comando dei tasti si riferiscono al supplemento S4.)
- Quindi, utilizzare la funzione verboseScatterplot per calcolare il coefficiente di correlazione del valore kME delle OTA corrispondenti nelle due reti e disegnare il diagramma di analisi di correlazione (Figura 4).
Eseguire la riga di comando:
"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)"
Fare clic su Esegui.
- Selezionare il modulo con il più piccolo coefficiente di correlazione del valore kME dell'OTU delle due reti. Considera questo modulo per avere la più grande differenza delle due reti.
5. Analisi del modulo di rete differenziale microbica
- Ottenere i dati dei phyla batterici dominanti attraverso l'analisi statistica del set di sequenze OTU del modulo con la più grande differenza.
NOTA: l'insieme di sequenze OTU del modulo con la differenza maggiore è sommato dalla tassonomia dei phyla. I batteri dominanti phyla rappresentavano oltre il 10%.
- Quindi, utilizzare la funzione exportNetworkToCytoscape per ottenere il file contenente le informazioni sulla relazione di interazione dell'OTU nel modulo differenziale più grande.
Eseguire la riga di comando:
"cyt = exportNetworkToCytoscape(modTOM,
edgeFile = paste("NEW-ES_CytoscapeInput-edges-", modules , ".txt", sep=""),
nodeFile = paste("NEW-ES_CytoscapeInput-nodes-", modules, ".txt", sep=""),
ponderato = TRUE,soglia = 0,5, nodeNames = modProbes,
altNodeNames = modGenes, nodeAttr = moduleColors[inModule])"
Fare clic su Esegui.
- Importare il file in Cytoscape. Impostare la soglia su 0,5 e regolare altri parametri in base alle esigenze.
- Costruire una rete di co-occorrenza di microrganismi differenziali (Figura 5).
- Ottenute le informazioni del genere core che ha il ruolo normativo più importante nella rete.
NOTA: In base al valore kME di OUT, è possibile definire il genere principale.
- Infine, sono state valutate le funzioni del genere principale ed è stata analizzata la sua influenza sull'intera rete di differenze.