{"id":329,"date":"2018-04-17T01:33:03","date_gmt":"2018-04-16T17:33:03","guid":{"rendered":"http:\/\/www.wuchangsong.com\/?p=329"},"modified":"2021-05-16T22:55:50","modified_gmt":"2021-05-16T14:55:50","slug":"wgcna%ef%bc%88my-project%ef%bc%89","status":"publish","type":"post","link":"http:\/\/www.wuchangsong.com\/?p=329","title":{"rendered":"WGCNA\uff08my project\uff09"},"content":{"rendered":"<pre>setwd(\"D:\/Ma_transcription\/WGCNA\/\")\ndatabase &lt;- read.table(file = \"COS_deseq_counts_normalized.txt\", sep = \"\\t\", header = T, row.names = 1, stringsAsFactors = F)\nlibrary(reshape2)\nlibrary('WGCNA')\nenableWGCNAThreads()#\u6253\u5f00\u591a\u7ebf\u7a0b\nWGCNA_matrix = t(database[order(apply(database,1,mad), decreasing = T)[1:10000],])\nsubname=sapply(colnames(database),function(x) strsplit(x,\"_\")[[1]][1])\ndatTraits = data.frame(gsm=names(database),\nsubtype=subname)\nrownames(datTraits)=datTraits[,1]\nhead(datTraits)\n# Choose a set of soft-thresholding powers\npowers = c(c(1:10), seq(from = 12, to=20, by=2))\ndatExpr &lt;- WGCNA_matrix\n# Call the network topology analysis function\nsft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)\n# Plot the results:\npar(mfrow = c(1,2));\ncex1 = 0.9;\n# Scale-free topology fit index as a function of the soft-thresholding power\nplot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],\n     xlab=\"Soft Threshold (power)\",ylab=\"Scale Free Topology Model Fit,signed R^2\",type=\"n\",\n     main = paste(\"Scale independence\"));\ntext(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],\n     labels=powers,cex=cex1,col=\"red\");\n# this line corresponds to using an R^2 cut-off of h\nabline(h=0.90,col=\"green\")\n# Mean connectivity as a function of the soft-thresholding power\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, cex=cex1,col=\"red\")\nsft$powerEstimate\nnet = blockwiseModules(datExpr, power = sft$powerEstimate,\n                       maxBlockSize = 5000,TOMType = \"unsigned\", \n                       minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,\n                       numericLabels = TRUE, pamRespectsDendro = FALSE,\n                       saveTOMs = TRUE,\n                       saveTOMFileBase = \"AS-green-FPKM-TOM\",\n                       verbose = 3)\ntable(net$colors)\n\n# Convert labels to colors for plotting\nmergedColors = labels2colors(net$colors)\ntable(mergedColors)\n# Plot the dendrogram and the module colors underneath\nplotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],\n                    \"Module colors\",\n                    dendroLabels = FALSE, hang = 0.03,\n                    addGuide = TRUE, guideHang = 0.05)\n## assign all of the gene to their corresponding module \n## hclust for the genes.\n#\u660e\u786e\u6837\u672c\u6570\u548c\u57fa\u56e0\u6570\nnGenes = ncol(datExpr)\nnSamples = nrow(datExpr)\n#\u9996\u5148\u9488\u5bf9\u6837\u672c\u505a\u4e2a\u7cfb\u7edf\u805a\u7c7b\u6811\ndatExpr_tree&lt;-hclust(dist(datExpr), method = \"average\")\npar(mar = c(0,5,2,0))\nplot(datExpr_tree, main = \"Sample clustering\", sub=\"\", xlab=\"\", cex.lab = 2, \n     cex.axis = 1, cex.main = 1,cex.lab=1)\n## \u5982\u679c\u8fd9\u4e2a\u65f6\u5019\u6837\u672c\u662f\u6709\u6027\u72b6\uff0c\u6216\u8005\u4e34\u5e8a\u8868\u578b\u7684\uff0c\u53ef\u4ee5\u52a0\u8fdb\u53bb\u770b\u770b\u662f\u5426\u805a\u7c7b\u5408\u7406\n#\u9488\u5bf9\u524d\u9762\u6784\u9020\u7684\u6837\u54c1\u77e9\u9635\u6dfb\u52a0\u5bf9\u5e94\u989c\u8272\nsample_colors &lt;- numbers2colors(as.numeric(factor(datTraits$subtype)), \n                                colors = c(\"grey\",\"blue\",\"red\",\"green\"),signed = FALSE)\n## \u8fd9\u4e2a\u7ed9\u6837\u54c1\u6dfb\u52a0\u5bf9\u5e94\u989c\u8272\u7684\u4ee3\u7801\u9700\u8981\u81ea\u884c\u4fee\u6539\u4ee5\u9002\u5e94\u81ea\u5df1\u7684\u6570\u636e\u5206\u6790\u9879\u76ee\u3002\n#  sample_colors &lt;- numbers2colors( datTraits ,signed = FALSE)\n## \u5982\u679c\u6837\u54c1\u6709\u591a\u79cd\u5206\u7c7b\u60c5\u51b5\uff0c\u800c\u4e14 datTraits \u91cc\u9762\u90fd\u662f\u5206\u7c7b\u4fe1\u606f\uff0c\u90a3\u4e48\u53ef\u4ee5\u76f4\u63a5\u7528\u4e0a\u9762\u4ee3\u7801\uff0c\n##\u5f53\u7136\uff0c\u8fd9\u6837\u7ed9\u7684\u989c\u8272\u4e0d\u660e\u663e\uff0c\u610f\u4e49\u4e0d\u5927\u3002\n#\u6784\u902010\u4e2a\u6837\u54c1\u7684\u7cfb\u7edf\u805a\u7c7b\u6811\u53ca\u6027\u72b6\u70ed\u56fe\npar(mar = c(1,4,3,1),cex=0.8)\nplotDendroAndColors(datExpr_tree, sample_colors,\n                    groupLabels = colnames(sample),\n                    cex.dendroLabels = 0.8,\n                    marAll = c(1, 4, 3, 1),\n                    cex.rowText = 0.01,\n                    main = \"Sample dendrogram and trait heatmap\")\n\ndesign=model.matrix(~0+ datTraits$subtype)\ncolnames(design)=levels(datTraits$subtype)\nmoduleColors &lt;- labels2colors(net$colors)\n# Recalculate MEs with color labels\nMEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes\nMEs = orderMEs(MEs0); ##\u4e0d\u540c\u989c\u8272\u7684\u6a21\u5757\u7684ME\u503c\u77e9\u9635(\u6837\u672cvs\u6a21\u5757)\nmoduleTraitCor = cor(MEs, design , use = \"p\");\nmoduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)\n\nsizeGrWindow(10,6)\n# Will display correlations and their p-values\ntextMatrix = paste(signif(moduleTraitCor, 2), \"\\n(\",\n                   signif(moduleTraitPvalue, 1), \")\", sep = \"\");\ndim(textMatrix) = dim(moduleTraitCor)\npar(mar = c(6, 8.5, 3, 3));\n# Display the correlation values within a heatmap plot\nlabeledHeatmap(Matrix = moduleTraitCor,\n               xLabels = colnames(design),\n               yLabels = names(MEs),\n               ySymbols = names(MEs),\n               colorLabels = FALSE,\n               colors = greenWhiteRed(50),\n               textMatrix = textMatrix,\n               setStdMargins = FALSE,\n               cex.text = 0.5,\n               zlim = c(-1,1),\n               main = paste(\"Module-trait relationships\"))                 \nhead(design)\n# names (colors) of the modules\nmodNames = substring(names(MEs), 3)\ngeneModuleMembership = as.data.frame(cor(datExpr, MEs, use = \"p\"));\n## \u7b97\u51fa\u6bcf\u4e2a\u6a21\u5757\u8ddf\u57fa\u56e0\u7684\u76ae\u5c14\u68ee\u76f8\u5173\u7cfb\u6570\u77e9\u9635\n## MEs\u662f\u6bcf\u4e2a\u6a21\u5757\u5728\u6bcf\u4e2a\u6837\u672c\u91cc\u9762\u7684\u503c\n## datExpr\u662f\u6bcf\u4e2a\u57fa\u56e0\u5728\u6bcf\u4e2a\u6837\u672c\u7684\u8868\u8fbe\u91cf\nMMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));\nnames(geneModuleMembership) = paste(\"MM\", modNames, sep=\"\");\nnames(MMPvalue) = paste(\"p.MM\", modNames, sep=\"\");\n## \u53ea\u6709\u8fde\u7eed\u578b\u6027\u72b6\u624d\u80fd\u53ea\u6709\u8ba1\u7b97\n## \u8fd9\u91cc\u628a\u662f\u5426\u5c5e\u4e8e CB \u8868\u578b\u8fd9\u4e2a\u53d8\u91cf\u75280,1\u8fdb\u884c\u6570\u503c\u5316\u3002\nCB = as.data.frame(design[,2]);\nnames(CB) = \"CB\"\ngeneTraitSignificance = as.data.frame(cor(datExpr, CB, use = \"p\"));\nGSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));\nnames(geneTraitSignificance) = paste(\"GS.\", names(CB), sep=\"\");\nnames(GSPvalue) = paste(\"p.GS.\", names(CB), sep=\"\")\nmodule = \"turquoise\"\ncolumn = match(module, modNames);\nmoduleGenes = moduleColors==module;\nsizeGrWindow(7, 7);\npar(mfrow = c(1,1));\nverboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),\n                   abs(geneTraitSignificance[moduleGenes, 1]),\n                   xlab = paste(\"Module Membership in\", module, \"module\"),\n                   ylab = \"Gene significance for CB\",\n                   main = paste(\"Module membership vs. gene significance\\n\"),\n                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)\n\n#\u9996\u5148\u9488\u5bf9\u6240\u6709\u57fa\u56e0\u753b\u70ed\u56fe\nnGenes = ncol(datExpr)\nnSamples = nrow(datExpr)\ngeneTree = net$dendrograms[[1]]; \n#\u751f\u6210\u5168\u57fa\u56e0\u4e0d\u76f8\u4f3cTOM\u77e9\u9635\ndissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6); \nplotTOM = dissTOM^7; \n#\u5982\u679c\u8868\u8fbe\u77e9\u9635\u6709\u5f88\u591a0\u503c\uff0c\u4f1a\u5bfc\u81f4TOM\u77e9\u9635\u4e3aNaN\uff0c\u6b64\u65f6\u5e94\u8fc7\u6ee4\u6570\u636e\ndiag(plotTOM) = NA; \n#TOMplot(plotTOM, geneTree, moduleColors, main = \"Network heatmap plot, all genes\")\n#pdf(\"11.pdf\")\n#TOMplot(plotTOM, net$dendrograms, moduleColors, main = \"Network heatmap plot, all genes\")\n#dev.off()\n#\u7136\u540e\u968f\u673a\u9009\u53d6\u90e8\u5206\u57fa\u56e0\u4f5c\u56fe\nnSelect = 400\n# For reproducibility, we set the random seed\nset.seed(10);\nselect = sample(nGenes, size = nSelect);\nselectTOM = dissTOM[select, select];\n# There\u2019s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.\nselectTree = hclust(as.dist(selectTOM), method = \"average\")\nselectColors = moduleColors[select];\n# Open a graphical window\nsizeGrWindow(9,9)\n# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing\n# the color palette; setting the diagonal to NA also improves the clarity of the plot\nplotDiss = selectTOM^7;\ndiag(plotDiss) = NA;\nTOMplot(plotDiss, selectTree, selectColors, main = \"Network heatmap plot, selected genes\")\n\n#\u6700\u540e\u753b\u6a21\u5757\u548c\u6027\u72b6\u7684\u5173\u7cfb\n# Recalculate module eigengenes\nMEs = moduleEigengenes(datExpr, moduleColors)$eigengenes\n## \u53ea\u6709\u8fde\u7eed\u578b\u6027\u72b6\u624d\u80fd\u53ea\u6709\u8ba1\u7b97\n## \u8fd9\u91cc\u628a\u662f\u5426\u5c5e\u4e8e Luminal \u8868\u578b\u8fd9\u4e2a\u53d8\u91cf\u75280,1\u8fdb\u884c\u6570\u503c\u5316\u3002\nCB = as.data.frame(design[,2]);\nnames(CB) = \"CB\"\n# Add the weight to existing module eigengenes\nMET = orderMEs(cbind(MEs, CB))\n# Plot the relationships among the eigengenes and the trait\nsizeGrWindow(5,7.5);\npar(cex = 0.9)\nplotEigengeneNetworks(MET, \"\", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8,                         xLabelsAngle= 90)\n# Plot the dendrogram\nsizeGrWindow(6,6);\npar(cex = 1.0)\n## \u6a21\u5757\u7684\u805a\u7c7b\u56fe\nplotEigengeneNetworks(MET, \"Eigengene dendrogram\", marDendro = c(0,4,2,0),\n                      plotHeatmaps = FALSE)\n# Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)\npar(cex = 1.0)\n## \u6027\u72b6\u4e0e\u6a21\u5757\u70ed\u56fe\nplotEigengeneNetworks(MET, \"Eigengene adjacency heatmap\", marHeatmap = c(3,4,2,2),\n                      plotDendrograms = FALSE, xLabelsAngle = 90)\n# Recalculate topological overlap\nTOM = TOMsimilarityFromExpr(datExpr, power = 6); \n# Select module\nmodule = \"blue\";\n# Select module probes\nprobes = colnames(datExpr) ## \u6211\u4eec\u4f8b\u5b50\u91cc\u9762\u7684probe\u5c31\u662f\u57fa\u56e0\u540d\ninModule = (moduleColors==module);\nmodProbes = probes[inModule]; \n## \u4e5f\u662f\u63d0\u53d6\u6307\u5b9a\u6a21\u5757\u7684\u57fa\u56e0\u540d\n# Select the corresponding Topological Overlap\nmodTOM = TOM[inModule, inModule];\ndimnames(modTOM) = list(modProbes, modProbes)\nnTop = 200;\nIMConn = softConnectivity(datExpr[, modProbes]);\ntop = (rank(-IMConn) &lt;= nTop)\nvis = exportNetworkToCytoscape(modTOM[top, top], edgeFile = paste(\"CytoscapeInput-edges-top200\", paste(module, collapse=\"-\"), \".txt\", sep=\"\"), nodeFile = paste(\"CytoscapeInput-nodes-top200\", paste(module, collapse=\"-\"), \".txt\", sep=\"\"), weighted = TRUE, threshold = 0);\n<\/pre>\n<p>\u53c2\u8003\uff1a<a href=\"https:\/\/mp.weixin.qq.com\/s?__biz=MzAxMDkxODM1Ng==&amp;mid=2247485962&amp;idx=1&amp;sn=2b662ed3e7dbdaa36b565ff0f73b40f5&amp;chksm=9b484ab1ac3fc3a7f67cabf1339e0d8cbd365adbbf4b746afa9910b95728885e8aa6f93630f2&amp;mpshare=1&amp;scene=23&amp;srcid=0217zXnGAiJMTZ1EJpEDqG9g#rd\">lncRNA\u5b9e\u6218\u9879\u76ee-\u7b2c\u516d\u6b65-WGCNA\u76f8\u5173\u6027\u5206\u6790<\/a><\/p>\n<p>&nbsp;<\/p>\n","protected":false},"excerpt":{"rendered":"<p>setwd(&#8220;D:\/Ma_transcription\/WGCNA\/&#8221;) database &lt;- read.table(file = &#8220;COS_deseq_counts_normalized.txt&#8221;, sep = &#8220;\\t&#8221;, header = T, row.names = 1, stringsAsFactors = F) library(reshape2) library(&#8216;WGCNA&#8217;) enableWGCNAThreads()#\u6253\u5f00\u591a\u7ebf\u7a0b WGCNA_matrix = t(database[order(apply(database,1,mad), decreasing = T)[1:10000],]) subname=sapply(colnames(database),function(x) strsplit(x,&#8221;_&#8221;)[[1]][1]) datTraits = data.frame(gsm=names(database), subtype=subname) rownames(datTraits)=datTraits[,1] head(datTraits) # Choose a set of soft-thresholding powers powers = c(c(1:10), seq(from = 12, to=20, by=2)) datExpr &lt;- WGCNA_matrix # Call [&hellip;]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[4,10],"tags":[],"_links":{"self":[{"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=\/wp\/v2\/posts\/329"}],"collection":[{"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=329"}],"version-history":[{"count":13,"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=\/wp\/v2\/posts\/329\/revisions"}],"predecessor-version":[{"id":768,"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=\/wp\/v2\/posts\/329\/revisions\/768"}],"wp:attachment":[{"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=329"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=329"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/www.wuchangsong.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=329"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}