WGCNA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the current working directory\n",
    "getwd();\n",
    "# If necessary, change the path below to the directory where the data files are stored.\n",
    "# \".\" means current directory. On Windows use a forward slash / instead of the usual \\.\n",
    "workingDir = \"/home/jovyan/work/datasets/paper/gtexhk/\";\n",
    "setwd(workingDir);\n",
    "# Load WGCNA package\n",
    "library(WGCNA)\n",
    "allowWGCNAThreads()\n",
    "# The following setting is important, do not omit.\n",
    "options(stringsAsFactors = FALSE);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "add_legend <- function(...) {\n",
    "  opar <- par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0), \n",
    "    mar=c(0, 0, 0, 0), new=TRUE)\n",
    "  on.exit(par(opar))\n",
    "  plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')\n",
    "  legend(...)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datExpr <- t(read.csv(\"mainTable.csv\", check.names=FALSE, row.names = \"\"))\n",
    "genes <- dimnames(datExpr)[[2]]\n",
    "samples <- dimnames(datExpr)[[1]]\n",
    "datExpr <- as.data.frame(apply(datExpr,2, as.numeric)) #numeric\n",
    "dimnames(datExpr)[[1]] <- samples\n",
    "allTraits <- read.csv(\"files.dat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log2_1 <- function(x){\n",
    "    log2(x+1)\n",
    "}\n",
    "datExpr <- as.data.frame(apply(datExpr,2, log2_1))\n",
    "dimnames(datExpr)[[1]] <- samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datExpr[1:10,1:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(allTraits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels <- c('SMTS', 'SMTSD')\n",
    "datTraits <- allTraits[,labels]\n",
    "rownames(datTraits) <- allTraits$SAMPID #<--CHECK\n",
    "datTraits <- datTraits[rownames(datExpr),]\n",
    "head(datTraits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table( dimnames(datExpr)[[1]]==datTraits$SMTS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = datTraits$SMTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "meanExpressionBySample=apply( datExpr,1,mean, na.rm=T)\n",
    "NumberMissingBySample=apply( is.na(data.frame(datExpr)),1, sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gsg = goodSamplesGenes(datExpr, verbose = 3);\n",
    "gsg$allOK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sampleTree = hclust(dist(datExpr), method = \"average\");\n",
    "# Plot the sample tree: Open a graphic output window of size 12 by 9 inches\n",
    "# The user should change the dimensions if the window is too large or too small.\n",
    "sizeGrWindow(12,9)\n",
    "svg(\"samplehier.svg\")\n",
    "par(cex = 0.6);\n",
    "par(mar = c(0,4,2,0))\n",
    "plot(sampleTree, main = \"Sample clustering to detect outliers\", sub=\"\", xlab=\"\", cex.lab = 1.5,\n",
    "cex.axis = 1.5, cex.main = 2)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine cluster under the line\n",
    "clust = cutreeStatic(sampleTree, cutHeight = 1e5, minSize = 10)\n",
    "table(clust)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clust 1 contains the samples we want to keep.\n",
    "keepSamples = (clust==1)\n",
    "datExpr = datExpr[keepSamples, ]\n",
    "nGenes = ncol(datExpr)\n",
    "nSamples = nrow(datExpr)\n",
    "\n",
    "genes <- dimnames(datExpr)[[2]]\n",
    "samples <- dimnames(datExpr)[[1]]\n",
    "datExpr <- as.data.frame(apply(datExpr,2, as.numeric)) #numeric\n",
    "dimnames(datExpr)[[1]] <- samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head(datExpr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rownames(allTraits) <- allTraits$SAMPID\n",
    "datTraits <- allTraits[keepSamples,labels]\n",
    "rownames(datTraits) <- rownames(allTraits[keepSamples,labels])\n",
    "datTraits <- datTraits[rownames(datExpr),]\n",
    "head(datTraits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "collectGarbage();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-cluster samples\n",
    "sampleTree2 = hclust(dist(datExpr), method = \"average\")\n",
    "fontsize=.5\n",
    "# Convert traits to a color representation: white means low, red means high, grey means missing entry\n",
    "traitColors = labels2colors(datTraits);\n",
    "# Plot the sample dendrogram and the colors underneath.\n",
    "pdf(\"samplehiercut.pdf\")\n",
    "plotDendroAndColors(sampleTree2, traitColors,\n",
    "    groupLabels = labels,\n",
    "    main = \"Sample dendrogram and trait heatmap\",\n",
    "    dendroLabels=FALSE,\n",
    "    rowText=datTraits,\n",
    "    rowTextIgnore=labels[[2]]\n",
    "    )\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose a set of soft-thresholding powers\n",
    "powers = c(c(1:10), seq(from = 12, to=50, by=2))\n",
    "# Call the network topology analysis function\n",
    "sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)\n",
    "# Plot the results:\n",
    "sizeGrWindow(9, 5)\n",
    "pdf(\"scaletopology.pdf\")\n",
    "par(mfrow = c(1,2));\n",
    "cex1 = 0.9;\n",
    "# Scale-free topology fit index as a function of the soft-thresholding power\n",
    "plot(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\"));\n",
    "text(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\n",
    "abline(h=0.90,col=\"red\")\n",
    "# Mean connectivity as a function of the soft-thresholding power\n",
    "plot(sft$fitIndices[,1], sft$fitIndices[,5],\n",
    "xlab=\"Soft Threshold (power)\",ylab=\"Mean Connectivity\", type=\"n\",\n",
    "main = paste(\"Mean connectivity\"))\n",
    "text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col=\"red\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "net = blockwiseModules(datExpr, power = 30,\n",
    "TOMType = \"unsigned\", minModuleSize = 5,\n",
    "reassignThreshold = 0, mergeCutHeight = 0.2,\n",
    "numericLabels = TRUE, pamRespectsDendro = FALSE,\n",
    "saveTOMs = TRUE,\n",
    "saveTOMFileBase = \"wgcna\",\n",
    "verbose = 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# open a graphics window\n",
    "sizeGrWindow(12, 9)\n",
    "pdf(\"netdendograms.pdf\")\n",
    "# Convert labels to colors for plotting\n",
    "mergedColors = labels2colors(net$colors)\n",
    "# Plot the dendrogram and the module colors underneath\n",
    "plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],\n",
    "\"Module colors\",\n",
    "dendroLabels = FALSE, hang = 0.03,\n",
    "addGuide = TRUE, guideHang = 0.05)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moduleLabels = net$colors\n",
    "moduleColors = labels2colors(net$colors)\n",
    "MEs = net$MEs;\n",
    "geneTree = net$dendrograms[[1]];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datTraits$SMTS <- match(datTraits$SMTS, unique(datTraits$SMTS))\n",
    "datTraits$SMTSD <- match(datTraits$SMTSD, unique(datTraits$SMTSD))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define numbers of genes and samples\n",
    "nGenes = ncol(datExpr);\n",
    "nSamples = nrow(datExpr);\n",
    "# Recalculate MEs with color labels\n",
    "MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes\n",
    "MEs = orderMEs(MEs0)\n",
    "moduleTraitCor = cor(MEs, datTraits, use = \"p\");\n",
    "moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf(\"moduleTrait.pdf\", width=9, height=18)\n",
    "# Will display correlations and their p-values\n",
    "textMatrix = paste(signif(moduleTraitCor, 2), \"\\n(\",\n",
    "signif(moduleTraitPvalue, 1), \")\", sep = \"\");\n",
    "dim(textMatrix) = dim(moduleTraitCor)\n",
    "par(mar = c(6, 8.5, 3, 3));\n",
    "# Display the correlation values within a heatmap plot\n",
    "labeledHeatmap(Matrix = moduleTraitCor,\n",
    "xLabels = labels,\n",
    "yLabels = names(MEs),\n",
    "ySymbols = names(MEs),\n",
    "colorLabels = FALSE,\n",
    "colors = blueWhiteRed(50),\n",
    "textMatrix = textMatrix,\n",
    "setStdMargins = FALSE,\n",
    "cex.text = 0.5,\n",
    "zlim = c(-1,1),\n",
    "main = paste(\"Module-trait relationships\"))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variable weight containing the weight column of datTrait\n",
    "weight = as.data.frame(rownames(datTraits));\n",
    "names(weight) = \"weight\"\n",
    "# names (colors) of the modules\n",
    "modNames = substring(names(MEs), 3)\n",
    "geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = \"p\"));\n",
    "MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));\n",
    "names(geneModuleMembership) = paste(\"MM\", modNames, sep=\"\");\n",
    "names(MMPvalue) = paste(\"p.MM\", modNames, sep=\"\");\n",
    "geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = \"p\"));\n",
    "GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));\n",
    "names(geneTraitSignificance) = paste(\"GS.\", names(weight), sep=\"\");\n",
    "names(GSPvalue) = paste(\"p.GS.\", names(weight), sep=\"\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "module = \"gray\"\n",
    "column = match(module, modNames);\n",
    "moduleGenes = moduleColors==module;\n",
    "sizeGrWindow(7, 7);\n",
    "pdf(\"modulemembership.pdf\")\n",
    "par(mfrow = c(1,1));\n",
    "verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),\n",
    "    abs(geneTraitSignificance[moduleGenes, 1]),\n",
    "    xlab = paste(\"Module Membership in\", module, \"module\"),\n",
    "    ylab = \"Gene significance for primary_site\",\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",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the starting data frame\n",
    "geneInfo0 = data.frame(\n",
    "moduleColor = moduleColors,\n",
    "geneTraitSignificance,\n",
    "GSPvalue)\n",
    "# Order modules by their significance for weight\n",
    "modOrder = order(-abs(cor(MEs, weight, use = \"p\")));\n",
    "# Add module membership information in the chosen order\n",
    "for (mod in 1:ncol(geneModuleMembership))\n",
    "{\n",
    "oldNames = names(geneInfo0)\n",
    "geneInfo0 = data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]],\n",
    "MMPvalue[, modOrder[mod]]);\n",
    "names(geneInfo0) = c(oldNames, paste(\"MM.\", modNames[modOrder[mod]], sep=\"\"),\n",
    "paste(\"p.MM.\", modNames[modOrder[mod]], sep=\"\"))\n",
    "}\n",
    "# Order the genes in the geneInfo variable first by module color, then by geneTraitSignificance\n",
    "geneOrder = order(geneInfo0$moduleColor, -abs(geneInfo0$GS.weight));\n",
    "geneInfo = geneInfo0[geneOrder, ]\n",
    "write.csv(geneInfo, file = \"geneInfo.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-cluster samples\n",
    "sampleTree3 = hclust(dist(MEs0), method = \"average\")\n",
    "# Convert traits to a color representation: white means low, red means high, grey means missing entry\n",
    "traitColors = labels2colors(datTraits);\n",
    "# Plot the sample dendrogram and the colors underneath.\n",
    "pdf(\"nethier.pdf\", width=10)\n",
    "# Set the plot margin: bottom, left, top & right\n",
    "par(mar = c(5, 4, 1.4, 0.2) + 0.1,\n",
    "    xpd = NA) # allow content to go into outer margin \n",
    "plotDendroAndColors(sampleTree3, traitColors,\n",
    "groupLabels = labels,\n",
    "main = \"Network dendrogram and trait heatmap\",\n",
    "dendroLabels=FALSE,\n",
    "#rowText=allTraits[keepSamples,labels],\n",
    "#rowTextIgnore=labels[[2]],\n",
    "addGuide=TRUE)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rownames(allTraits)=allTraits$file_name\n",
    "unique(allTraits[rownames(datTraits),'Subtype_Selected'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svg(\"nethier_legend.svg\", width=10)\n",
    "add_legend(\"topleft\", legend = unique(allTraits[rownames(datTraits),'SMTS']), pch = 1, pt.cex = 1, cex = 1.5, bty = 'n', ncol=2,\n",
    "       inset = c(0.1, 0.), # place outside\n",
    "       title = \"SMTS\", \n",
    "       col = traitColors)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "collectGarbage()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clust = cutreeStatic(sampleTree3, cutHeight = 0.4, minSize = 10)\n",
    "table(clust)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Re-cluster samples after cut\n",
    "sampleTree4 = hclust(dist(MEs0[(clust==1),]), method = \"average\")\n",
    "# Convert traits to a color representation: white means low, red means high, grey means missing entry\n",
    "traitColors = labels2colors(datTraits[(clust==1),]);\n",
    "# Plot the sample dendrogram and the colors underneath.\n",
    "svg(\"nethiercut.svg\")\n",
    "plotDendroAndColors(sampleTree4, traitColors,\n",
    "groupLabels = labels,\n",
    "main = \"Network dendrogram and trait heatmap\",\n",
    "dendroLabels=FALSE)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate topological overlap anew: this could be done more efficiently by saving the TOM\n",
    "# calculated during module detection, but let us do it again here.\n",
    "dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);\n",
    "# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap\n",
    "plotTOM = dissTOM^7;\n",
    "# Set diagonal to NA for a nicer plot\n",
    "diag(plotTOM) = NA;\n",
    "# Call the plot function\n",
    "sizeGrWindow(9,9)\n",
    "svg(\"networkmap.svg\")\n",
    "TOMplot(plotTOM, geneTree, moduleColors, main = \"Network heatmap plot, all genes\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#selct genes to be faster\n",
    "nSelect = nGenes*0.7\n",
    "# For reproducibility, we set the random seed\n",
    "set.seed(10);\n",
    "select = sample(nGenes, size = nSelect);\n",
    "selectTOM = dissTOM[select, select];\n",
    "# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.\n",
    "selectTree = hclust(as.dist(selectTOM), method = \"average\")\n",
    "selectColors = moduleColors[select];\n",
    "# Open a graphical window\n",
    "sizeGrWindow(9,9)\n",
    "svg(\"networkmap.svg\")\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\n",
    "plotDiss = selectTOM^7;\n",
    "diag(plotDiss) = NA;\n",
    "TOMplot(plotDiss, selectTree, selectColors, main = \"Network heatmap plot, selected genes\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recalculate module eigengenes\n",
    "MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes\n",
    "# Isolate weight from the clinical traits\n",
    "weight = as.data.frame(datTraits$SMTS);\n",
    "names(weight) = \"SMTS\"\n",
    "# Add the weight to existing module eigengenes\n",
    "MET = orderMEs(cbind(MEs, weight))\n",
    "# Plot the relationships among the eigengenes and the trait\n",
    "sizeGrWindow(5,7.5);\n",
    "svg(\"eigenHeat.svg\")\n",
    "par(cex = 0.9)\n",
    "plotEigengeneNetworks(MET, \"\", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle = 90)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the dendrogram\n",
    "sizeGrWindow(6,6);\n",
    "svg(\"eigendendogram.svg\")\n",
    "par(cex = 1.0)\n",
    "plotEigengeneNetworks(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)\n",
    "par(cex = 1.0)\n",
    "plotEigengeneNetworks(MET, \"Eigengene adjacency heatmap\", marHeatmap = c(3,4,2,2),\n",
    "plotDendrograms = FALSE, xLabelsAngle = 90)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Values of *k* must be read from *clustersizes.txt file*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write.csv(cutree(sampleTree3, k=2), \"wgcna/wgcna_level_4_labels.csv\")\n",
    "write.csv(cutree(sampleTree3, k=2), \"wgcna/wgcna_level_3_labels.csv\")\n",
    "write.csv(cutree(sampleTree3, k=6), \"wgcna/wgcna_level_2_labels.csv\")\n",
    "write.csv(cutree(sampleTree3, k=13), \"wgcna/wgcna_level_1_labels.csv\")\n",
    "write.csv(cutree(sampleTree3, k=256), \"wgcna/wgcna_level_0_labels.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write.csv(geneModuleMembership, \"wgcna/wgcna_level_0_word-dist.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "write.csv(MEs, \"wgcna/wgcna_level_0_topic-dist.csv\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}