{"id":101,"date":"2021-10-20T20:44:49","date_gmt":"2021-10-20T12:44:49","guid":{"rendered":"https:\/\/linguopeng.top\/?p=101"},"modified":"2021-10-20T20:51:53","modified_gmt":"2021-10-20T12:51:53","slug":"wgcna%e5%ae%9e%e6%88%98","status":"publish","type":"post","link":"https:\/\/linguopeng.top\/?p=101","title":{"rendered":"WGCNA\u5b9e\u8df5"},"content":{"rendered":"\n

Langfelder P, Horvath <\/a>S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics, 2008, 9(559):1-13.<\/p>\n\n\n\n

#WGCNA \u5305\u901a\u8fc7 bioconductor \u5b89\u88c5\n#install.packages('BiocManager') \n#BiocManager::install('WGCNA')\nlibrary(WGCNA)\n#\u793a\u4f8b\u6570\u636e\uff0c\u57fa\u56e0\u8868\u8fbe\u503c\u77e9\u9635\ngene <- read.delim('gene_express.txt', row.names = 1, check.names = FALSE)\n#\u4f8b\u5982\u8fc7\u6ee4\u4f4e\u8868\u8fbe\u503c\u7684\u57fa\u56e0\uff0c\u53ea\u4fdd\u7559\u5e73\u5747\u8868\u8fbe\u503c\u5728 1 \u4ee5\u4e0a\u7684\ngene <- subset(gene, rowSums(gene)\/ncol(gene) >= 1)\n#\u8f6c\u7f6e\u77e9\u9635\uff0c\u884c\u662f\u57fa\u56e0\uff0c\u5217\u662f\u6837\u672c\ngene <- t(gene)\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
##power \u503c\u9009\u62e9\npowers <- 1:20\nsft <- pickSoftThreshold(gene, powerVector = powers, verbose = 5)\n#\u62df\u5408\u6307\u6570\u4e0e power \u503c\u6563\u70b9\u56fe\npar(mfrow = c(1, 2))\nplot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], type = 'n', \n    xlab = 'Soft Threshold (power)', ylab = 'Scale Free Topology Model Fit,signed R^2', \n    main = paste('Scale independence'))\ntext(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],\n    labels = powers, col = 'red');\nabline(h = 0.90, col = 'red')\n#\u5e73\u5747\u8fde\u901a\u6027\u4e0e power \u503c\u6563\u70b9\u56fe\nplot(sft$fitIndices[,1], sft$fitIndices[,5],\n    xlab = 'Soft Threshold (power)', ylab = 'Mean Connectivity', type = 'n',\n    main = paste('Mean connectivity'))\ntext(sft$fitIndices[,1], sft$fitIndices[,5], labels = powers, col = 'red')\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u4e0a\u4e00\u6b65\u4f30\u8ba1\u7684\u6700\u4f73 power \u503c\npowers <- sft$powerEstimate\n#\u83b7\u5f97 TOM \u77e9\u9635\nadjacency <- adjacency(gene, power = powers)\ntom_sim <- TOMsimilarity(adjacency)\nrownames(tom_sim) <- rownames(adjacency)\ncolnames(tom_sim) <- colnames(adjacency)\ntom_sim[1:6,1:6]\n#\u8f93\u51fa TOM \u77e9\u9635\n#write.table(tom_sim, 'TOMsimilarity.txt', sep = '\\t', col.names = NA, quote = FALSE)\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#TOM \u76f8\u5f02\u5ea6 = 1 \u2013 TOM \u76f8\u4f3c\u5ea6\ntom_dis  <- 1 - tom_sim\n#\u5c42\u6b21\u805a\u7c7b\u6811\uff0c\u4f7f\u7528\u4e2d\u503c\u7684\u975e\u6743\u91cd\u6210\u5bf9\u7ec4\u6cd5\u7684\u5e73\u5747\u805a\u5408\u805a\u7c7b\ngeneTree <- hclust(as.dist(tom_dis), method = 'average')\nplot(geneTree, xlab = '', sub = '', main = 'Gene clustering on TOM-based dissimilarity',\n    labels = FALSE, hang = 0.04)\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u4f7f\u7528\u52a8\u6001\u526a\u5207\u6811\u6316\u6398\u6a21\u5757\nminModuleSize <- 30  #\u6a21\u5757\u57fa\u56e0\u6570\u76ee\ndynamicMods <- cutreeDynamic(dendro = geneTree, distM = tom_dis,\n    deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize)\ntable(dynamicMods)\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u6a21\u5757\u989c\u8272\u6307\u4ee3\ndynamicColors <- labels2colors(dynamicMods)\ntable(dynamicColors)\nplotDendroAndColors(geneTree, dynamicColors, 'Dynamic Tree Cut',\n    dendroLabels = FALSE, addGuide = TRUE, hang = 0.03, guideHang = 0.05,\n    main = 'Gene dendrogram and module colors')\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u57fa\u56e0\u8868\u8fbe\u805a\u7c7b\u6811\u548c\u5171\u8868\u8fbe\u62d3\u6251\u70ed\u56fe\nplot_sim <- -(1-tom_sim)\n#plot_sim <- log(tom_sim)\ndiag(plot_sim) <- NA\nTOMplot(plot_sim, geneTree, dynamicColors, \n    main = 'Network heatmap plot, selected genes')\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u8ba1\u7b97\u57fa\u56e0\u8868\u8fbe\u77e9\u9635\u4e2d\u6a21\u5757\u7684\u7279\u5f81\u57fa\u56e0\uff08\u7b2c\u4e00\u4e3b\u6210\u5206\uff09\nMEList <- moduleEigengenes(gene, colors = dynamicColors)\nMEs <- MEList$eigengenes\nhead(MEs)[1:6]\n#\u8f93\u51fa\u6a21\u5757\u7279\u5f81\u57fa\u56e0\u77e9\u9635\n#write.table(MEs, 'moduleEigengenes.txt', sep = '\\t', col.names = NA, quote = FALSE)\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u901a\u8fc7\u6a21\u5757\u7279\u5f81\u57fa\u56e0\u8ba1\u7b97\u6a21\u5757\u95f4\u76f8\u5173\u6027\uff0c\u8868\u5f81\u6a21\u5757\u95f4\u76f8\u4f3c\u5ea6\nME_cor <- cor(MEs)\nME_cor[1:6,1:6]\n#\u7ed8\u5236\u805a\u7c7b\u6811\u89c2\u5bdf\nMETree <- hclust(as.dist(1-ME_cor), method = 'average')\nplot(METree, main = 'Clustering of module eigengenes', xlab = '', sub = '')\n#\u63a2\u7d22\u6027\u5206\u6790\uff0c\u89c2\u5bdf\u6a21\u5757\u95f4\u7684\u76f8\u4f3c\u6027\n#height \u503c\u53ef\u4ee3\u8868\u6a21\u5757\u95f4\u7684\u76f8\u5f02\u5ea6\uff0c\u5e76\u786e\u5b9a\u4e00\u4e2a\u5408\u9002\u7684\u9608\u503c\u4f5c\u4e3a\u526a\u5207\u9ad8\u5ea6\n#\u4ee5\u4fbf\u4e3a\u4f4e\u76f8\u5f02\u5ea6\uff08\u9ad8\u76f8\u4f3c\u5ea6\uff09\u7684\u6a21\u5757\u5408\u5e76\u63d0\u4f9b\u4f9d\u636e\nabline(h = 0.2, col = 'blue')\nabline(h = 0.25, col = 'red')\n#\u6a21\u5757\u7279\u5f81\u57fa\u56e0\u805a\u7c7b\u6811\u70ed\u56fe\nplotEigengeneNetworks(MEs, '', cex.lab = 0.8, xLabelsAngle= 90, \n    marDendro = c(0, 4, 1, 2), marHeatmap = c(3, 4, 1, 2))\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u76f8\u4f3c\u6a21\u5757\u5408\u5e76\uff0c\u4ee5 0.25 \u4f5c\u4e3a\u5408\u5e76\u9608\u503c\uff08\u526a\u5207\u9ad8\u5ea6\uff09\uff0c\u5728\u6b64\u9ad8\u5ea6\u4e0b\u7684\u6a21\u5757\u5c06\u5408\u5e76\n#\u8fd1\u4f3c\u7406\u89e3\u4e3a\u76f8\u5173\u7a0b\u5ea6\u9ad8\u4e8e 0.75 \u7684\u6a21\u5757\u5c06\u5408\u5e76\u5230\u4e00\u8d77\nmerge_module <- mergeCloseModules(gene, dynamicColors, cutHeight = 0.25, verbose = 3)\nmergedColors <- merge_module$colors\ntable(mergedColors)\n#\u57fa\u56e0\u8868\u8fbe\u548c\u6a21\u5757\u805a\u7c7b\u6811\nplotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c('Dynamic Tree Cut', 'Merged dynamic'),\n    dendroLabels = FALSE, addGuide = TRUE, hang = 0.03, guideHang = 0.05)\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u60a3\u8005\u7684\u4e34\u5e8a\u8868\u578b\u6570\u636e\ntrait <- read.delim('trait.txt', row.names = 1, check.names = FALSE)\n#\u4f7f\u7528\u4e0a\u4e00\u6b65\u65b0\u7ec4\u5408\u7684\u5171\u8868\u8fbe\u6a21\u5757\u7684\u7ed3\u679c\nmodule <- merge_module$newMEs\n#\u60a3\u8005\u57fa\u56e0\u5171\u8868\u8fbe\u6a21\u5757\u548c\u4e34\u5e8a\u8868\u578b\u7684\u76f8\u5173\u6027\u5206\u6790\nmoduleTraitCor <- cor(module, trait, use = 'p')\nmoduleTraitCor[1:6,1:6]  #\u76f8\u5173\u77e9\u9635\n#\u76f8\u5173\u7cfb\u6570\u7684 p \u503c\u77e9\u9635\nmoduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(module))\n#\u8f93\u51fa\u76f8\u5173\u7cfb\u6570\u77e9\u9635\u6216 p \u503c\u77e9\u9635\n#write.table(moduleTraitCor, 'moduleTraitCor.txt', sep = '\\t', col.names = NA, quote = FALSE)\n#write.table(moduleTraitPvalue, 'moduleTraitPvalue.txt', sep = '\\t', col.names = NA, quote = FALSE)\n#\u76f8\u5173\u56fe\u7ed8\u5236\ntextMatrix <- paste(signif(moduleTraitCor, 2), '\\n(', signif(moduleTraitPvalue, 1), ')', sep = '')\ndim(textMatrix) <- dim(moduleTraitCor)\npar(mar = c(5, 10, 3, 3))\nlabeledHeatmap(Matrix = moduleTraitCor, main = paste('Module-trait relationships'), \n    xLabels = names(trait), yLabels = names(module), ySymbols = names(module),\n    colorLabels = FALSE, colors = greenWhiteRed(50), cex.text = 0.7, zlim = c(-1,1), \n    textMatrix = textMatrix, setStdMargins = FALSE)\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u57fa\u56e0\u4e0e\u6a21\u5757\u7684\u5bf9\u5e94\u5173\u7cfb\u5217\u8868\ngene_module <- data.frame(gene_name = colnames(gene), module = mergedColors, stringsAsFactors = FALSE)\nhead(gene_module)\n#\u201cblack\u201d\u6a21\u5757\u5185\u7684\u57fa\u56e0\u540d\u79f0\ngene_module_select <- subset(gene_module, module == 'black')$gene_name\n#\u201cblack\u201d\u6a21\u5757\u5185\u57fa\u56e0\u5728\u5404\u6837\u672c\u4e2d\u7684\u8868\u8fbe\u503c\u77e9\u9635\uff08\u57fa\u56e0\u8868\u8fbe\u503c\u77e9\u9635\u7684\u4e00\u4e2a\u5b50\u96c6\uff09\ngene_select <- t(gene[,gene_module_select])\n#\u201cblack\u201d\u6a21\u5757\u5185\u57fa\u56e0\u7684\u5171\u8868\u8fbe\u76f8\u4f3c\u5ea6\uff08\u5728 TOM \u77e9\u9635\u4e2d\u63d0\u53d6\u5b50\u96c6\uff09\ntom_select <- tom_sim[gene_module_select,gene_module_select]\n#\u8f93\u51fa\n#write.table(gene_select, 'gene_select.txt', sep = '\\t', col.names = NA, quote = FALSE)\n#write.table(tom_select, 'tom_select.txt', sep = '\\t', col.names = NA, quote = FALSE)\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u9009\u62e9 black \u6a21\u5757\u5185\u7684\u57fa\u56e0\ngene_black <- gene[ ,gene_module_select]\n#\u57fa\u56e0\u7684\u6a21\u5757\u6210\u5458\u5ea6\uff08module membership\uff09\u8ba1\u7b97\n#\u5373\u5404\u57fa\u56e0\u8868\u8fbe\u503c\u4e0e\u76f8\u5e94\u6a21\u5757\u7279\u5f81\u57fa\u56e0\u7684\u76f8\u5173\u6027\uff0c\u5176\u8861\u91cf\u4e86\u57fa\u56e0\u5728\u5168\u5c40\u7f51\u7edc\u4e2d\u7684\u4f4d\u7f6e\ngeneModuleMembership <- signedKME(gene_black, module['MEblack'], outputColumnName = 'MM')\nMMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nrow(module)))\n#\u5404\u57fa\u56e0\u8868\u8fbe\u503c\u4e0e\u4e34\u5e8a\u8868\u578b\u7684\u76f8\u5173\u6027\u5206\u6790\ngeneTraitSignificance <- as.data.frame(cor(gene_black, trait['TNM'], use = 'p'))\nGSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nrow(trait)))\n#\u9009\u62e9\u663e\u8457\uff08p<0.01\uff09\u3001\u9ad8 black \u6a21\u5757\u6210\u5458\u5ea6\uff08MM>=0.8\uff09\uff0c\u4e0e TNM \u8868\u578b\u9ad8\u5ea6\u76f8\u5173\uff08r>=0.8\uff09\u7684\u57fa\u56e0\ngeneModuleMembership[geneModuleMembership<0.8 | MMPvalue>0.01] <- 0\ngeneTraitSignificance[geneTraitSignificance<0.8 | MMPvalue>0.01] <- 0\nselect <- cbind(geneModuleMembership, geneTraitSignificance)\nselect <- subset(select, geneModuleMembership>=0.8 & geneTraitSignificance>=0.8)\nhead(select)<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u5019\u9009\u57fa\u56e0\u4e0e TNM \u5206\u671f\u7684\u76f8\u5173\u6027\u6563\u70b9\u56fe\ndir.create('black_TNM_cor', recursive = TRUE)\nfor (i in rownames(select)) {\n    png(paste('black_TNM_cor\/', i, '-TNM.png', sep = ''), \n        width = 4, height = 4, res = 400, units = 'in', type = 'cairo')\n    plot(trait[ ,'TNM'], gene[ ,i], xlab = 'TNM', ylab = i, \n        main = paste('MM_black = ', select[i,'MMblack'], '\\ncor_TNM = ', select[i,'TNM']), cex = 0.8, pch = 20)\n    fit <- lm(gene[ ,i]~trait[ ,'TNM'])\n    lines(trait[ ,'TNM'], fitted(fit))\n    dev.off()\n}\n#\u5019\u9009\u57fa\u56e0\u95f4\u7684 TOM \u77e9\u9635\nplotNetworkHeatmap(gene,\n    plotGenes = rownames(select),\n    networkType = 'unsigned',\n    useTOM = TRUE,\n    power = powers,\n    main = 'unsigned correlations')\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
#\u8f93\u51fa\u5404\u6a21\u5757\u5b50\u7f51\u7edc\u7684\u8fb9\u5217\u8868\u548c\u8282\u70b9\u5217\u8868\uff0c\u540e\u7eed\u53ef\u5bfc\u5165 cytoscape \u7ed8\u5236\u7f51\u7edc\u56fe\n#\u5176\u4e2d\uff0c\u8fb9\u7684\u6743\u91cd\u4e3a TOM \u77e9\u9635\u4e2d\u7684\u503c\uff0c\u8bb0\u5f55\u4e86\u57fa\u56e0\u95f4\u5171\u8868\u8fbe\u76f8\u4f3c\u6027\ndir.create('cytoscape', recursive = TRUE)\nfor (mod in 1:nrow(table(mergedColors))) {\n    modules <- names(table(mergedColors))[mod]\n    probes <- colnames(gene)\n    inModule <- (mergedColors == modules)\n    modProbes <- probes[inModule]\n    modGenes <- modProbes\n    modtom_sim <- tom_sim[inModule, inModule] \n    dimnames(modtom_sim) <- list(modProbes, modProbes)\n    outEdge <- paste('cytoscape\/', modules, '.edge_list.txt',sep = '')\n    outNode <- paste('cytoscape\/', modules, '.node_list.txt', sep = '')\n    exportNetworkToCytoscape(modtom_sim,\n        edgeFile = outEdge,\n        nodeFile = outNode,\n        weighted = TRUE,\n        threshold = 0.3,  #\u8be5\u53c2\u6570\u53ef\u63a7\u5236\u8f93\u51fa\u7684\u8fb9\u6570\u91cf\uff0c\u8fc7\u6ee4\u4f4e\u6743\u91cd\u7684\u8fb9\n        nodeNames = modProbes,\n        altNodeNames = modGenes,\n        nodeAttr = mergedColors[inModule])\n}\n<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
<\/p>\n","protected":false},"excerpt":{"rendered":"
Langfelder P, Horvath S. WGCNA: an R package for weight […]<\/p>\n","protected":false},"author":1,"featured_media":115,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[3],"class_list":["post-101","post","type-post","status-publish","format-standard","has-post-thumbnail","hentry","category-uncategorized","tag-r"],"_links":{"self":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/101","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=101"}],"version-history":[{"count":3,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/101\/revisions"}],"predecessor-version":[{"id":123,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/101\/revisions\/123"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/media\/115"}],"wp:attachment":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=101"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=101"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=101"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}