• WGCNA教程 官网教程 总结版 代码tutorials 1

    1. Data input and cleaning: PDF document, R script
    2. Network construction and module detection
      1. Automatic, one-step network construction and module detection: PDF document, R script
      2. Step-by-step network construction and module detection: PDF document, R script
      3. Dealing with large datasets: block-wise network construction and module detection: PDF document, R script
    3. Relating modules to external clinical traits and identifying important genes: PDF document, R script
    4. Interfacing network analysis with other data such as functional annotation and gene ontology PDF document, R script
    5. Network visualization using WGCNA functions: PDF document, R script
    6. Export of networks to external software: PDF document, R script

    1. Data input and cleaning: PDF document, R script

    1. #=====================================================================================
    2. #
    3. #  Code chunk 1
    4. #
    5. #=====================================================================================
    6.  
    7.  
    8. # Display the current working directory
    9. getwd();
    10. # If necessary, change the path below to the directory where the data files are stored. 
    11. # "." means current directory. On Windows use a forward slash / instead of the usual \.
    12. workingDir = ".";
    13. setwd(workingDir); 
    14. # Load the WGCNA package
    15. library(WGCNA);
    16. # The following setting is important, do not omit.
    17. options(stringsAsFactors = FALSE);
    18. #Read in the female liver data set
    19. femData = read.csv("LiverFemale3600.csv");
    20. # Take a quick look at what is in the data set:
    21. dim(femData);
    22. names(femData);
    23.  
    24.  
    25. #=====================================================================================
    26. #
    27. #  Code chunk 2
    28. #
    29. #=====================================================================================
    30.  
    31.  
    32. datExpr0 = as.data.frame(t(femData[, -c(1:8)]));
    33. names(datExpr0) = femData$substanceBXH;
    34. rownames(datExpr0) = names(femData)[-c(1:8)];
    35.  
    36.  
    37. #=====================================================================================
    38. #
    39. #  Code chunk 3
    40. #
    41. #=====================================================================================
    42.  
    43.  
    44. gsg = goodSamplesGenes(datExpr0, verbose = 3);
    45. gsg$allOK
    46.  
    47.  
    48. #=====================================================================================
    49. #
    50. #  Code chunk 4
    51. #
    52. #=====================================================================================
    53.  
    54.  
    55. if (!gsg$allOK)
    56. {
    57.   # Optionally, print the gene and sample names that were removed:
    58.   if (sum(!gsg$goodGenes)>0) 
    59.      printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
    60.   if (sum(!gsg$goodSamples)>0) 
    61.      printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
    62.   # Remove the offending genes and samples from the data:
    63.   datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
    64. }
    65.  
    66.  
    67. #=====================================================================================
    68. #
    69. #  Code chunk 5
    70. #
    71. #=====================================================================================
    72.  
    73.  
    74. sampleTree = hclust(dist(datExpr0), method = "average");
    75. # Plot the sample tree: Open a graphic output window of size 12 by 9 inches
    76. # The user should change the dimensions if the window is too large or too small.
    77. sizeGrWindow(12,9)
    78. #pdf(file = "Plots/sampleClustering.pdf", width = 12, height = 9);
    79. par(cex = 0.6);
    80. par(mar = c(0,4,2,0))
    81. plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, 
    82.      cex.axis = 1.5, cex.main = 2)
    83.  
    84.  
    85. #=====================================================================================
    86. #
    87. #  Code chunk 6
    88. #
    89. #=====================================================================================
    90.  
    91.  
    92. # Plot a line to show the cut
    93. abline(h = 15, col = "red");
    94. # Determine cluster under the line
    95. clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
    96. table(clust)
    97. # clust 1 contains the samples we want to keep.
    98. keepSamples = (clust==1)
    99. datExpr = datExpr0[keepSamples, ]
    100. nGenes = ncol(datExpr)
    101. nSamples = nrow(datExpr)
    102.  
    103.  
    104. #=====================================================================================
    105. #
    106. #  Code chunk 7
    107. #
    108. #=====================================================================================
    109.  
    110.  
    111. traitData = read.csv("ClinicalTraits.csv");
    112. dim(traitData)
    113. names(traitData)
    114.  
    115. # remove columns that hold information we do not need.
    116. allTraits = traitData[, -c(31, 16)];
    117. allTraits = allTraits[, c(2, 11:36) ];
    118. dim(allTraits)
    119. names(allTraits)
    120.  
    121. # Form a data frame analogous to expression data that will hold the clinical traits.
    122.  
    123. femaleSamples = rownames(datExpr);
    124. traitRows = match(femaleSamples, allTraits$Mice);
    125. datTraits = allTraits[traitRows, -1];
    126. rownames(datTraits) = allTraits[traitRows, 1];
    127.  
    128. collectGarbage();
    129.  
    130.  
    131. #=====================================================================================
    132. #
    133. #  Code chunk 8
    134. #
    135. #=====================================================================================
    136.  
    137.  
    138. # Re-cluster samples
    139. sampleTree2 = hclust(dist(datExpr), method = "average")
    140. # Convert traits to a color representation: white means low, red means high, grey means missing entry
    141. traitColors = numbers2colors(datTraits, signed = FALSE);
    142. # Plot the sample dendrogram and the colors underneath.
    143. plotDendroAndColors(sampleTree2, traitColors,
    144.                     groupLabels = names(datTraits), 
    145.                     main = "Sample dendrogram and trait heatmap")
    146.  
    147.  
    148. #=====================================================================================
    149. #
    150. #  Code chunk 9
    151. #
    152. #=====================================================================================
    153.  
    154.  
    155. save(datExpr, datTraits, file = "FemaleLiver-01-dataInput.RData")
    156.  
    157.

     Tutorials for WGCNA R package

    step2.

    Network construction and module detection

    1. Automatic, one-step network construction and module detection: PDF document, R script

    1. #=====================================================================================
    2. #
    3. #  Code chunk 1
    4. #
    5. #=====================================================================================
    6.  
    7.  
    8. # Display the current working directory
    9. getwd();
    10. # If necessary, change the path below to the directory where the data files are stored. 
    11. # "." means current directory. On Windows use a forward slash / instead of the usual \.
    12. workingDir = ".";
    13. setwd(workingDir); 
    14. # Load the WGCNA package
    15. library(WGCNA)
    16. # The following setting is important, do not omit.
    17. options(stringsAsFactors = FALSE);
    18. # Allow multi-threading within WGCNA. This helps speed up certain calculations.
    19. # At present this call is necessary for the code to work.
    20. # Any error here may be ignored but you may want to update WGCNA if you see one.
    21. # Caution: skip this line if you run RStudio or other third-party R environments. 
    22. # See note above.
    23. enableWGCNAThreads()
    24. # Load the data saved in the first part
    25. lnames = load(file = "FemaleLiver-01-dataInput.RData");
    26. #The variable lnames contains the names of loaded variables.
    27. lnames
    28.  
    29.  
    30. #=====================================================================================
    31. #
    32. #  Code chunk 2
    33. #
    34. #=====================================================================================
    35.  
    36.  
    37. # Choose a set of soft-thresholding powers
    38. powers = c(c(1:10), seq(from = 12, to=20, by=2))
    39. # Call the network topology analysis function
    40. sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
    41. # Plot the results:
    42. sizeGrWindow(9, 5)
    43. par(mfrow = c(1,2));
    44. cex1 = 0.9;
    45. # Scale-free topology fit index as a function of the soft-thresholding power
    46. plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    47.      xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
    48.      main = paste("Scale independence"));
    49. text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    50.      labels=powers,cex=cex1,col="red");
    51. # this line corresponds to using an R^2 cut-off of h
    52. abline(h=0.90,col="red")
    53. # Mean connectivity as a function of the soft-thresholding power
    54. plot(sft$fitIndices[,1], sft$fitIndices[,5],
    55.      xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
    56.      main = paste("Mean connectivity"))
    57. text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
    58.  
    59.  
    60. #=====================================================================================
    61. #
    62. #  Code chunk 3
    63. #
    64. #=====================================================================================
    65.  
    66.  
    67. net = blockwiseModules(datExpr, power = 6,
    68.                        TOMType = "unsigned", minModuleSize = 30,
    69.                        reassignThreshold = 0, mergeCutHeight = 0.25,
    70.                        numericLabels = TRUE, pamRespectsDendro = FALSE,
    71.                        saveTOMs = TRUE,
    72.                        saveTOMFileBase = "femaleMouseTOM", 
    73.                        verbose = 3)
    74.  
    75.  
    76. #=====================================================================================
    77. #
    78. #  Code chunk 4
    79. #
    80. #=====================================================================================
    81.  
    82.  
    83. # open a graphics window
    84. sizeGrWindow(12, 9)
    85. # Convert labels to colors for plotting
    86. mergedColors = labels2colors(net$colors)
    87. # Plot the dendrogram and the module colors underneath
    88. plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
    89.                     "Module colors",
    90.                     dendroLabels = FALSE, hang = 0.03,
    91.                     addGuide = TRUE, guideHang = 0.05)
    92.  
    93.  
    94. #=====================================================================================
    95. #
    96. #  Code chunk 5
    97. #
    98. #=====================================================================================
    99.  
    100.  
    101. moduleLabels = net$colors
    102. moduleColors = labels2colors(net$colors)
    103. MEs = net$MEs;
    104. geneTree = net$dendrograms[[1]];
    105. save(MEs, moduleLabels, moduleColors, geneTree, 
    106.      file = "FemaleLiver-02-networkConstruction-auto.RData")
    107.  
    108.

    1. Step-by-step network construction and module detection: PDF document, R script

    1. #=====================================================================================
    2. #
    3. #  Code chunk 1
    4. #
    5. #=====================================================================================
    6.  
    7.  
    8. # Display the current working directory
    9. getwd();
    10. # If necessary, change the path below to the directory where the data files are stored. 
    11. # "." means current directory. On Windows use a forward slash / instead of the usual \.
    12. workingDir = ".";
    13. setwd(workingDir); 
    14. # Load the WGCNA package
    15. library(WGCNA)
    16. # The following setting is important, do not omit.
    17. options(stringsAsFactors = FALSE);
    18. # Allow multi-threading within WGCNA. This helps speed up certain calculations.
    19. # At present this call is necessary for the code to work.
    20. # Any error here may be ignored but you may want to update WGCNA if you see one.
    21. # Caution: skip this line if you run RStudio or other third-party R environments. 
    22. # See note above.
    23. enableWGCNAThreads()
    24. # Load the data saved in the first part
    25. lnames = load(file = "FemaleLiver-01-dataInput.RData");
    26. #The variable lnames contains the names of loaded variables.
    27. lnames
    28.  
    29.  
    30. #=====================================================================================
    31. #
    32. #  Code chunk 2
    33. #
    34. #=====================================================================================
    35.  
    36.  
    37. # Choose a set of soft-thresholding powers
    38. powers = c(c(1:10), seq(from = 12, to=20, by=2))
    39. # Call the network topology analysis function
    40. sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
    41. # Plot the results:
    42. sizeGrWindow(9, 5)
    43. par(mfrow = c(1,2));
    44. cex1 = 0.9;
    45. # Scale-free topology fit index as a function of the soft-thresholding power
    46. plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    47.      xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
    48.      main = paste("Scale independence"));
    49. text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    50.      labels=powers,cex=cex1,col="red");
    51. # this line corresponds to using an R^2 cut-off of h
    52. abline(h=0.90,col="red")
    53. # Mean connectivity as a function of the soft-thresholding power
    54. plot(sft$fitIndices[,1], sft$fitIndices[,5],
    55.      xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
    56.      main = paste("Mean connectivity"))
    57. text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
    58.  
    59.  
    60. #=====================================================================================
    61. #
    62. #  Code chunk 3
    63. #
    64. #=====================================================================================
    65.  
    66.  
    67. net = blockwiseModules(datExpr, power = 6,
    68.                        TOMType = "unsigned", minModuleSize = 30,
    69.                        reassignThreshold = 0, mergeCutHeight = 0.25,
    70.                        numericLabels = TRUE, pamRespectsDendro = FALSE,
    71.                        saveTOMs = TRUE,
    72.                        saveTOMFileBase = "femaleMouseTOM", 
    73.                        verbose = 3)
    74.  
    75.  
    76. #=====================================================================================
    77. #
    78. #  Code chunk 4
    79. #
    80. #=====================================================================================
    81.  
    82.  
    83. # open a graphics window
    84. sizeGrWindow(12, 9)
    85. # Convert labels to colors for plotting
    86. mergedColors = labels2colors(net$colors)
    87. # Plot the dendrogram and the module colors underneath
    88. plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
    89.                     "Module colors",
    90.                     dendroLabels = FALSE, hang = 0.03,
    91.                     addGuide = TRUE, guideHang = 0.05)
    92.  
    93.  
    94. #=====================================================================================
    95. #
    96. #  Code chunk 5
    97. #
    98. #=====================================================================================
    99.  
    100.  
    101. moduleLabels = net$colors
    102. moduleColors = labels2colors(net$colors)
    103. MEs = net$MEs;
    104. geneTree = net$dendrograms[[1]];
    105. save(MEs, moduleLabels, moduleColors, geneTree, 
    106.      file = "FemaleLiver-02-networkConstruction-auto.RData")
    107.  
    108.

    Dealing with large datasets: block-wise network construction and module detection: PDF document, R script

    1. #=====================================================================================
    2. #
    3. #  Code chunk 1
    4. #
    5. #=====================================================================================
    6.  
    7.  
    8. # Display the current working directory
    9. getwd();
    10. # If necessary, change the path below to the directory where the data files are stored. 
    11. # "." means current directory. On Windows use a forward slash / instead of the usual \.
    12. workingDir = ".";
    13. setwd(workingDir); 
    14. # Load the WGCNA package
    15. library(WGCNA)
    16. # The following setting is important, do not omit.
    17. options(stringsAsFactors = FALSE);
    18. # Allow multi-threading within WGCNA. This helps speed up certain calculations.
    19. # At present this call is necessary.
    20. # Any error here may be ignored but you may want to update WGCNA if you see one.
    21. # Caution: skip this line if you run RStudio or other third-party R environments.
    22. # See note above.
    23. enableWGCNAThreads()
    24. # Load the data saved in the first part
    25. lnames = load(file = "FemaleLiver-01-dataInput.RData");
    26. #The variable lnames contains the names of loaded variables.
    27. lnames
    28.  
    29.  
    30. #=====================================================================================
    31. #
    32. #  Code chunk 2
    33. #
    34. #=====================================================================================
    35.  
    36.  
    37. # Choose a set of soft-thresholding powers
    38. powers = c(c(1:10), seq(from = 12, to=20, by=2))
    39. # Call the network topology analysis function
    40. sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
    41. # Plot the results:
    42. sizeGrWindow(9, 5)
    43. par(mfrow = c(1,2));
    44. cex1 = 0.9;
    45. # Scale-free topology fit index as a function of the soft-thresholding power
    46. plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    47.      xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
    48.      main = paste("Scale independence"));
    49. text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    50.      labels=powers,cex=cex1,col="red");
    51. # this line corresponds to using an R^2 cut-off of h
    52. abline(h=0.90,col="red")
    53. # Mean connectivity as a function of the soft-thresholding power
    54. plot(sft$fitIndices[,1], sft$fitIndices[,5],
    55.      xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
    56.      main = paste("Mean connectivity"))
    57. text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
    58.  
    59.  
    60. #=====================================================================================
    61. #
    62. #  Code chunk 3
    63. #
    64. #=====================================================================================
    65.  
    66.  
    67. bwnet = blockwiseModules(datExpr, maxBlockSize = 2000,
    68.                        power = 6, TOMType = "unsigned", minModuleSize = 30,
    69.                        reassignThreshold = 0, mergeCutHeight = 0.25,
    70.                        numericLabels = TRUE,
    71.                        saveTOMs = TRUE,
    72.                        saveTOMFileBase = "femaleMouseTOM-blockwise",
    73.                        verbose = 3)
    74.  
    75.  
    76. #=====================================================================================
    77. #
    78. #  Code chunk 4
    79. #
    80. #=====================================================================================
    81.  
    82.  
    83. # Load the results of single-block analysis
    84. load(file = "FemaleLiver-02-networkConstruction-auto.RData");
    85. # Relabel blockwise modules
    86. bwLabels = matchLabels(bwnet$colors, moduleLabels);
    87. # Convert labels to colors for plotting
    88. bwModuleColors = labels2colors(bwLabels)
    89.  
    90.  
    91. #=====================================================================================
    92. #
    93. #  Code chunk 5
    94. #
    95. #=====================================================================================
    96.  
    97.  
    98. # open a graphics window
    99. sizeGrWindow(6,6)
    100. # Plot the dendrogram and the module colors underneath for block 1
    101. plotDendroAndColors(bwnet$dendrograms[[1]], bwModuleColors[bwnet$blockGenes[[1]]],
    102.                     "Module colors", main = "Gene dendrogram and module colors in block 1", 
    103.                     dendroLabels = FALSE, hang = 0.03,
    104.                     addGuide = TRUE, guideHang = 0.05)
    105. # Plot the dendrogram and the module colors underneath for block 2
    106. plotDendroAndColors(bwnet$dendrograms[[2]], bwModuleColors[bwnet$blockGenes[[2]]],
    107.                     "Module colors", main = "Gene dendrogram and module colors in block 2", 
    108.                     dendroLabels = FALSE, hang = 0.03,
    109.                     addGuide = TRUE, guideHang = 0.05)
    110.  
    111.  
    112. #=====================================================================================
    113. #
    114. #  Code chunk 6
    115. #
    116. #=====================================================================================
    117.  
    118.  
    119. sizeGrWindow(12,9)
    120. plotDendroAndColors(geneTree,
    121.                     cbind(moduleColors, bwModuleColors),
    122.                     c("Single block", "2 blocks"),
    123.                     main = "Single block gene dendrogram and module colors",
    124.                     dendroLabels = FALSE, hang = 0.03,
    125.                     addGuide = TRUE, guideHang = 0.05)
    126.  
    127.  
    128. #=====================================================================================
    129. #
    130. #  Code chunk 7
    131. #
    132. #=====================================================================================
    133.  
    134.  
    135. singleBlockMEs = moduleEigengenes(datExpr, moduleColors)$eigengenes;
    136. blockwiseMEs = moduleEigengenes(datExpr, bwModuleColors)$eigengenes;
    137.  
    138.  
    139. #=====================================================================================
    140. #
    141. #  Code chunk 8
    142. #
    143. #=====================================================================================
    144.  
    145.  
    146. single2blockwise = match(names(singleBlockMEs), names(blockwiseMEs))
    147. signif(diag(cor(blockwiseMEs[, single2blockwise], singleBlockMEs)), 3)
    148.  
    149.

    Relating modules to external clinical traits and identifying important genes: PDF document, R script

    1. #=====================================================================================
    2. #
    3. #  Code chunk 1
    4. #
    5. #=====================================================================================
    6.  
    7.  
    8. # Display the current working directory
    9. getwd();
    10. # If necessary, change the path below to the directory where the data files are stored. 
    11. # "." means current directory. On Windows use a forward slash / instead of the usual \.
    12. workingDir = ".";
    13. setwd(workingDir); 
    14. # Load the WGCNA package
    15. library(WGCNA)
    16. # The following setting is important, do not omit.
    17. options(stringsAsFactors = FALSE);
    18. # Load the expression and trait data saved in the first part
    19. lnames = load(file = "FemaleLiver-01-dataInput.RData");
    20. #The variable lnames contains the names of loaded variables.
    21. lnames
    22. # Load network data saved in the second part.
    23. lnames = load(file = "FemaleLiver-02-networkConstruction-auto.RData");
    24. lnames
    25.  
    26.  
    27. #=====================================================================================
    28. #
    29. #  Code chunk 2
    30. #
    31. #=====================================================================================
    32.  
    33.  
    34. # Define numbers of genes and samples
    35. nGenes = ncol(datExpr);
    36. nSamples = nrow(datExpr);
    37. # Recalculate MEs with color labels
    38. MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
    39. MEs = orderMEs(MEs0)
    40. moduleTraitCor = cor(MEs, datTraits, use = "p");
    41. moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
    42.  
    43.  
    44. #=====================================================================================
    45. #
    46. #  Code chunk 3
    47. #
    48. #=====================================================================================
    49.  
    50.  
    51. sizeGrWindow(10,6)
    52. # Will display correlations and their p-values
    53. textMatrix =  paste(signif(moduleTraitCor, 2), "\n(",
    54.                            signif(moduleTraitPvalue, 1), ")", sep = "");
    55. dim(textMatrix) = dim(moduleTraitCor)
    56. par(mar = c(6, 8.5, 3, 3));
    57. # Display the correlation values within a heatmap plot
    58. labeledHeatmap(Matrix = moduleTraitCor,
    59.                xLabels = names(datTraits),
    60.                yLabels = names(MEs),
    61.                ySymbols = names(MEs),
    62.                colorLabels = FALSE,
    63.                colors = greenWhiteRed(50),
    64.                textMatrix = textMatrix,
    65.                setStdMargins = FALSE,
    66.                cex.text = 0.5,
    67.                zlim = c(-1,1),
    68.                main = paste("Module-trait relationships"))
    69.  
    70.  
    71. #=====================================================================================
    72. #
    73. #  Code chunk 4
    74. #
    75. #=====================================================================================
    76.  
    77.  
    78. # Define variable weight containing the weight column of datTrait
    79. weight = as.data.frame(datTraits$weight_g);
    80. names(weight) = "weight"
    81. # names (colors) of the modules
    82. modNames = substring(names(MEs), 3)
    83.  
    84. geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
    85. MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
    86.  
    87. names(geneModuleMembership) = paste("MM", modNames, sep="");
    88. names(MMPvalue) = paste("p.MM", modNames, sep="");
    89.  
    90. geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = "p"));
    91. GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
    92.  
    93. names(geneTraitSignificance) = paste("GS.", names(weight), sep="");
    94. names(GSPvalue) = paste("p.GS.", names(weight), sep="");
    95.  
    96.  
    97. #=====================================================================================
    98. #
    99. #  Code chunk 5
    100. #
    101. #=====================================================================================
    102.  
    103.  
    104. module = "brown"
    105. column = match(module, modNames);
    106. moduleGenes = moduleColors==module;
    107.  
    108. sizeGrWindow(7, 7);
    109. par(mfrow = c(1,1));
    110. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
    111.                    abs(geneTraitSignificance[moduleGenes, 1]),
    112.                    xlab = paste("Module Membership in", module, "module"),
    113.                    ylab = "Gene significance for body weight",
    114.                    main = paste("Module membership vs. gene significance\n"),
    115.                    cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
    116.  
    117.  
    118. #=====================================================================================
    119. #
    120. #  Code chunk 6
    121. #
    122. #=====================================================================================
    123.  
    124.  
    125. names(datExpr)
    126.  
    127.  
    128. #=====================================================================================
    129. #
    130. #  Code chunk 7
    131. #
    132. #=====================================================================================
    133.  
    134.  
    135. names(datExpr)[moduleColors=="brown"]
    136.  
    137.  
    138. #=====================================================================================
    139. #
    140. #  Code chunk 8
    141. #
    142. #=====================================================================================
    143.  
    144.  
    145. annot = read.csv(file = "GeneAnnotation.csv");
    146. dim(annot)
    147. names(annot)
    148. probes = names(datExpr)
    149. probes2annot = match(probes, annot$substanceBXH)
    150. # The following is the number or probes without annotation:
    151. sum(is.na(probes2annot))
    152. # Should return 0.
    153.  
    154.  
    155. #=====================================================================================
    156. #
    157. #  Code chunk 9
    158. #
    159. #=====================================================================================
    160.  
    161.  
    162. # Create the starting data frame
    163. geneInfo0 = data.frame(substanceBXH = probes,
    164.                       geneSymbol = annot$gene_symbol[probes2annot],
    165.                       LocusLinkID = annot$LocusLinkID[probes2annot],
    166.                       moduleColor = moduleColors,
    167.                       geneTraitSignificance,
    168.                       GSPvalue)
    169. # Order modules by their significance for weight
    170. modOrder = order(-abs(cor(MEs, weight, use = "p")));
    171. # Add module membership information in the chosen order
    172. for (mod in 1:ncol(geneModuleMembership))
    173. {
    174.   oldNames = names(geneInfo0)
    175.   geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]], 
    176.                          MMPvalue[, modOrder[mod]]);
    177.   names(geneInfo0) = c(oldNames, paste("MM.", modNames[modOrder[mod]], sep=""),
    178.                        paste("p.MM.", modNames[modOrder[mod]], sep=""))
    179. }
    180. # Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance
    181. geneOrder = order(geneInfo0$moduleColor, -abs(geneInfo0$GS.weight));
    182. geneInfo = geneInfo0[geneOrder, ]
    183.  
    184.  
    185. #=====================================================================================
    186. #
    187. #  Code chunk 10
    188. #
    189. #=====================================================================================
    190.  
    191.  
    192. write.csv(geneInfo, file = "geneInfo.csv")
    193.  
    194.

    Interfacing network analysis with other data such as functional annotation and gene ontology PDF document, R script

    1. #=====================================================================================
    2. #
    3. #  Code chunk 1
    4. #
    5. #=====================================================================================
    6.  
    7.  
    8. # Display the current working directory
    9. getwd();
    10. # If necessary, change the path below to the directory where the data files are stored. 
    11. # "." means current directory. On Windows use a forward slash / instead of the usual \.
    12. workingDir = ".";
    13. setwd(workingDir); 
    14. # Load the WGCNA package
    15. library(WGCNA)
    16. # The following setting is important, do not omit.
    17. options(stringsAsFactors = FALSE);
    18. # Load the expression and trait data saved in the first part
    19. lnames = load(file = "FemaleLiver-01-dataInput.RData");
    20. #The variable lnames contains the names of loaded variables.
    21. lnames
    22. # Load network data saved in the second part.
    23. lnames = load(file = "FemaleLiver-02-networkConstruction-auto.RData");
    24. lnames
    25.  
    26.  
    27. #=====================================================================================
    28. #
    29. #  Code chunk 2
    30. #
    31. #=====================================================================================
    32.  
    33.  
    34. # Read in the probe annotation
    35. annot = read.csv(file = "GeneAnnotation.csv");
    36. # Match probes in the data set to the probe IDs in the annotation file 
    37. probes = names(datExpr)
    38. probes2annot = match(probes, annot$substanceBXH)
    39. # Get the corresponding Locuis Link IDs
    40. allLLIDs = annot$LocusLinkID[probes2annot];
    41. # $ Choose interesting modules
    42. intModules = c("brown", "red", "salmon")
    43. for (module in intModules)
    44. {
    45.   # Select module probes
    46.   modGenes = (moduleColors==module)
    47.   # Get their entrez ID codes
    48.   modLLIDs = allLLIDs[modGenes];
    49.   # Write them into a file
    50.   fileName = paste("LocusLinkIDs-", module, ".txt", sep="");
    51.   write.table(as.data.frame(modLLIDs), file = fileName,
    52.               row.names = FALSE, col.names = FALSE)
    53. }
    54. # As background in the enrichment analysis, we will use all probes in the analysis.
    55. fileName = paste("LocusLinkIDs-all.txt", sep="");
    56. write.table(as.data.frame(allLLIDs), file = fileName,
    57.             row.names = FALSE, col.names = FALSE)
    58.  
    59.  
    60. #=====================================================================================
    61. #
    62. #  Code chunk 3
    63. #
    64. #=====================================================================================
    65.  
    66.  
    67. GOenr = GOenrichmentAnalysis(moduleColors, allLLIDs, organism = "mouse", nBestP = 10);
    68.  
    69.  
    70. #=====================================================================================
    71. #
    72. #  Code chunk 4
    73. #
    74. #=====================================================================================
    75.  
    76.  
    77. tab = GOenr$bestPTerms[[4]]$enrichment
    78.  
    79.  
    80. #=====================================================================================
    81. #
    82. #  Code chunk 5
    83. #
    84. #=====================================================================================
    85.  
    86.  
    87. names(tab)
    88.  
    89.  
    90. #=====================================================================================
    91. #
    92. #  Code chunk 6
    93. #
    94. #=====================================================================================
    95.  
    96.  
    97. write.table(tab, file = "GOEnrichmentTable.csv", sep = ",", quote = TRUE, row.names = FALSE)
    98.  
    99.  
    100. #=====================================================================================
    101. #
    102. #  Code chunk 7
    103. #
    104. #=====================================================================================
    105.  
    106.  
    107. keepCols = c(1, 2, 5, 6, 7, 12, 13);
    108. screenTab = tab[, keepCols];
    109. # Round the numeric columns to 2 decimal places:
    110. numCols = c(3, 4);
    111. screenTab[, numCols] = signif(apply(screenTab[, numCols], 2, as.numeric), 2)
    112. # Truncate the the term name to at most 40 characters
    113. screenTab[, 7] = substring(screenTab[, 7], 1, 40)
    114. # Shorten the column names:
    115. colnames(screenTab) = c("module", "size", "p-val", "Bonf", "nInTerm", "ont", "term name");
    116. rownames(screenTab) = NULL;
    117. # Set the width of R's output. The reader should play with this number to obtain satisfactory output.
    118. options(width=95)
    119. # Finally, display the enrichment table:
    120. screenTab

    Network visualization using WGCNA functions: PDF document, R script

    1. #=====================================================================================
    2. #
    3. #  Code chunk 1
    4. #
    5. #=====================================================================================
    6.  
    7.  
    8. # Display the current working directory
    9. getwd();
    10. # If necessary, change the path below to the directory where the data files are stored. 
    11. # "." means current directory. On Windows use a forward slash / instead of the usual \.
    12. workingDir = ".";
    13. setwd(workingDir); 
    14. # Load the WGCNA package
    15. library(WGCNA)
    16. # The following setting is important, do not omit.
    17. options(stringsAsFactors = FALSE);
    18. # Load the expression and trait data saved in the first part
    19. lnames = load(file = "FemaleLiver-01-dataInput.RData");
    20. #The variable lnames contains the names of loaded variables.
    21. lnames
    22. # Load network data saved in the second part.
    23. lnames = load(file = "FemaleLiver-02-networkConstruction-auto.RData");
    24. lnames
    25. nGenes = ncol(datExpr)
    26. nSamples = nrow(datExpr)
    27.  
    28.  
    29. #=====================================================================================
    30. #
    31. #  Code chunk 2
    32. #
    33. #=====================================================================================
    34.  
    35.  
    36. # Calculate topological overlap anew: this could be done more efficiently by saving the TOM
    37. # calculated during module detection, but let us do it again here.
    38. dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
    39. # Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
    40. plotTOM = dissTOM^7;
    41. # Set diagonal to NA for a nicer plot
    42. diag(plotTOM) = NA;
    43. # Call the plot function
    44. sizeGrWindow(9,9)
    45. TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
    46.  
    47.  
    48. #=====================================================================================
    49. #
    50. #  Code chunk 3
    51. #
    52. #=====================================================================================
    53.  
    54.  
    55. nSelect = 400
    56. # For reproducibility, we set the random seed
    57. set.seed(10);
    58. select = sample(nGenes, size = nSelect);
    59. selectTOM = dissTOM[select, select];
    60. # There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
    61. selectTree = hclust(as.dist(selectTOM), method = "average")
    62. selectColors = moduleColors[select];
    63. # Open a graphical window
    64. sizeGrWindow(9,9)
    65. # Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing 
    66. # the color palette; setting the diagonal to NA also improves the clarity of the plot
    67. plotDiss = selectTOM^7;
    68. diag(plotDiss) = NA;
    69. TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
    72. #=====================================================================================
    73. #
    74. #  Code chunk 4
    75. #
    76. #=====================================================================================
    79. # Recalculate module eigengenes
    80. MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
    81. # Isolate weight from the clinical traits
    82. weight = as.data.frame(datTraits$weight_g);
    83. names(weight) = "weight"
    84. # Add the weight to existing module eigengenes
    85. MET = orderMEs(cbind(MEs, weight))
    86. # Plot the relationships among the eigengenes and the trait
    87. sizeGrWindow(5,7.5);
    88. par(cex = 0.9)
    89. plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
    90. = 90)
    93. #=====================================================================================
    94. #
    95. #  Code chunk 5
    96. #
    97. #=====================================================================================
    100. # Plot the dendrogram
    101. sizeGrWindow(6,6);
    102. par(cex = 1.0)
    103. plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),
    104.                       plotHeatmaps = FALSE)
    105. # Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
    106. par(cex = 1.0)
    107. plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2),
    108.                       plotDendrograms = FALSE, xLabelsAngle = 90)

    6.Export of networks to external software: PDF document, R script

    1. #=====================================================================================
    2. #
    3. #  Code chunk 1
    4. #
    5. #=====================================================================================
    6.  
    7.  
    8. # Display the current working directory
    9. getwd();
    10. # If necessary, change the path below to the directory where the data files are stored. 
    11. # "." means current directory. On Windows use a forward slash / instead of the usual \.
    12. workingDir = ".";
    13. setwd(workingDir); 
    14. # Load the WGCNA package
    15. library(WGCNA)
    16. # The following setting is important, do not omit.
    17. options(stringsAsFactors = FALSE);
    18. # Load the expression and trait data saved in the first part
    19. lnames = load(file = "FemaleLiver-01-dataInput.RData");
    20. #The variable lnames contains the names of loaded variables.
    21. lnames
    22. # Load network data saved in the second part.
    23. lnames = load(file = "FemaleLiver-02-networkConstruction-auto.RData");
    24. lnames
    25.  
    26.  
    27. #=====================================================================================
    28. #
    29. #  Code chunk 2
    30. #
    31. #=====================================================================================
    32.  
    33.  
    34. # Recalculate topological overlap
    35. TOM = TOMsimilarityFromExpr(datExpr, power = 6);
    36. # Read in the annotation file
    37. annot = read.csv(file = "GeneAnnotation.csv");
    38. # Select module
    39. module = "brown";
    40. # Select module probes
    41. probes = names(datExpr)
    42. inModule = (moduleColors==module);
    43. modProbes = probes[inModule];
    44. # Select the corresponding Topological Overlap
    45. modTOM = TOM[inModule, inModule];
    46. dimnames(modTOM) = list(modProbes, modProbes)
    47. # Export the network into an edge list file VisANT can read
    48. vis = exportNetworkToVisANT(modTOM,
    49.   file = paste("VisANTInput-", module, ".txt", sep=""),
    50.   weighted = TRUE,
    51.   threshold = 0,
    52.   probeToGene = data.frame(annot$substanceBXH, annot$gene_symbol) )
    53.  
    54.  
    55. #=====================================================================================
    56. #
    57. #  Code chunk 3
    58. #
    59. #=====================================================================================
    60.  
    61.  
    62. nTop = 30;
    63. IMConn = softConnectivity(datExpr[, modProbes]);
    64. top = (rank(-IMConn) <= nTop)
    65. vis = exportNetworkToVisANT(modTOM[top, top],
    66.   file = paste("VisANTInput-", module, "-top30.txt", sep=""),
    67.   weighted = TRUE,
    68.   threshold = 0,
    69.   probeToGene = data.frame(annot$substanceBXH, annot$gene_symbol) )
    70.  
    71.  
    72. #=====================================================================================
    73. #
    74. #  Code chunk 4
    75. #
    76. #=====================================================================================
    77.  
    78.  
    79. # Recalculate topological overlap if needed
    80. TOM = TOMsimilarityFromExpr(datExpr, power = 6);
    81. # Read in the annotation file
    82. annot = read.csv(file = "GeneAnnotation.csv");
    83. # Select modules
    84. modules = c("brown", "red");
    85. # Select module probes
    86. probes = names(datExpr)
    87. inModule = is.finite(match(moduleColors, modules));
    88. modProbes = probes[inModule];
    89. modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)];
    90. # Select the corresponding Topological Overlap
    91. modTOM = TOM[inModule, inModule];
    92. dimnames(modTOM) = list(modProbes, modProbes)
    93. # Export the network into edge and node list files Cytoscape can read
    94. cyt = exportNetworkToCytoscape(modTOM,
    95.   edgeFile = paste("CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
    96.   nodeFile = paste("CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
    97.   weighted = TRUE,
    98.   threshold = 0.02,
    99.   nodeNames = modProbes,
    100.   altNodeNames = modGenes,
    101.   nodeAttr = moduleColors[inModule]);
    102.  
    103.

  • 相关阅读:
    mysql统计整个数据库记录条数
    C++ 统计程序运行时间
    EF列表分页查询(单表、多表)，排除参数为空的条件
    【数据结构】| 并查集及其优化实现
    Dify后端源码目录结构和蓝图
    Vue 路由 使用
    CSS - 响应式布局（二）响应式栅格系统
    linux发展历程
    企望制造ERP存在远程命令执行漏洞 附POC
    从SAP ECC升级到SAP S4HANA, 几个Key Points
  • 原文地址：https://blog.csdn.net/qq_52813185/article/details/127651588