code/WGCNA_proteomics/wgcna.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b7a4b5f3-90f9-4270-b55c-e421fea75e5a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "── \u001b[1mAttaching packages\u001b[22m ─────────────────────────────────────── tidyverse 1.3.1 ──\n",
      "\n",
      "\u001b[32m✔\u001b[39m \u001b[34mggplot2\u001b[39m 3.3.5     \u001b[32m✔\u001b[39m \u001b[34mpurrr  \u001b[39m 0.3.4\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtibble \u001b[39m 3.1.5     \u001b[32m✔\u001b[39m \u001b[34mdplyr  \u001b[39m 1.0.7\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtidyr  \u001b[39m 1.1.4     \u001b[32m✔\u001b[39m \u001b[34mstringr\u001b[39m 1.4.0\n",
      "\u001b[32m✔\u001b[39m \u001b[34mreadr  \u001b[39m 2.0.2     \u001b[32m✔\u001b[39m \u001b[34mforcats\u001b[39m 0.5.1\n",
      "\n",
      "── \u001b[1mConflicts\u001b[22m ────────────────────────────────────────── tidyverse_conflicts() ──\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mfilter()\u001b[39m masks \u001b[34mstats\u001b[39m::filter()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mlag()\u001b[39m    masks \u001b[34mstats\u001b[39m::lag()\n",
      "\n",
      "\n",
      "Attaching package: ‘jsonlite’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:purrr’:\n",
      "\n",
      "    flatten\n",
      "\n",
      "\n",
      "Loading required package: dynamicTreeCut\n",
      "\n",
      "Loading required package: fastcluster\n",
      "\n",
      "\n",
      "Attaching package: ‘fastcluster’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:stats’:\n",
      "\n",
      "    hclust\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘WGCNA’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:stats’:\n",
      "\n",
      "    cor\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘flashClust’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:fastcluster’:\n",
      "\n",
      "    hclust\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:stats’:\n",
      "\n",
      "    hclust\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘matrixStats’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:dplyr’:\n",
      "\n",
      "    count\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "'/data/projects/robin/cell_free/proteomics_WGCNA'"
      ],
      "text/latex": [
       "'/data/projects/robin/cell\\_free/proteomics\\_WGCNA'"
      ],
      "text/markdown": [
       "'/data/projects/robin/cell_free/proteomics_WGCNA'"
      ],
      "text/plain": [
       "[1] \"/data/projects/robin/cell_free/proteomics_WGCNA\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(tidyverse)\n",
    "library(jsonlite)\n",
    "library(WGCNA)\n",
    "library(flashClust)\n",
    "library(matrixStats)\n",
    "\n",
    "# Set number of threads\n",
    "nThreads = 30\n",
    "\n",
    "# Set directlory for storing experiment results\n",
    "setwd(\"./\")\n",
    "\n",
    "getwd()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3e74d71f-d12b-4627-be97-a53e3c92fd22",
   "metadata": {},
   "outputs": [],
   "source": [
    "base_dir = \"/data/projects/robin/cell_free/proteomics_WGCNA/\"\n",
    "data = \"processed/lfq.csv\"\n",
    "metadata = \"processed/metadata.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b4e019e6-f8a7-468a-99d7-81934ae9dbdc",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj = paste(base_dir, data, sep=\"\")\n",
    "vst_mat = read.csv(obj, row.names=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5cd075e3-7063-48da-be6e-66376841e9fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "obj_metadata = paste(base_dir, metadata, sep=\"\")\n",
    "anno = read.csv(obj_metadata, row.names=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "11667657-596d-48ca-8f83-43717937903f",
   "metadata": {},
   "outputs": [],
   "source": [
    "vst_mat_corrected <- empiricalBayesLM(t(vst_mat), anno$batch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d4f3caa2-d5d4-4a0d-9258-6dea5c92ba73",
   "metadata": {},
   "outputs": [],
   "source": [
    "vst_mat <- t(vst_mat_corrected[[1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be001075-b0f7-4db3-b6fc-d5be1923513d",
   "metadata": {},
   "outputs": [],
   "source": [
    "anno$group <- anno[, \"Neuronal.subclass..0.low..1.high.\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b396b2e4-f8ca-46a3-99b1-656911b2c85b",
   "metadata": {},
   "outputs": [],
   "source": [
    "conditions<-factor(sapply(rownames(anno),function(id){paste(anno[id,c(\"group\")],collapse = '_')}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c9a5867f-04fb-4827-9e95-0d53cd5777f0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".dl-inline {width: auto; margin:0; padding: 0}\n",
       ".dl-inline>dt, .dl-inline>dd {float: none; width: auto; display: inline-block}\n",
       ".dl-inline>dt::after {content: \":\\0020\"; padding-right: .5ex}\n",
       ".dl-inline>dt:not(:first-of-type) {padding-left: .5ex}\n",
       "</style><dl class=dl-inline><dt>F1</dt><dd>1</dd><dt>F2</dt><dd>0</dd><dt>F3</dt><dd>1</dd><dt>F4</dt><dd>0</dd><dt>F5</dt><dd>0</dd><dt>F6</dt><dd>0</dd><dt>F7</dt><dd>0</dd><dt>F8</dt><dd>0</dd><dt>F9</dt><dd>0</dd><dt>F10</dt><dd>0</dd><dt>F11</dt><dd>0</dd><dt>F13</dt><dd>1</dd><dt>F14</dt><dd>0</dd><dt>F15</dt><dd>1</dd><dt>F16</dt><dd>1</dd><dt>F17</dt><dd>1</dd><dt>F18</dt><dd>1</dd><dt>F19</dt><dd>0</dd><dt>F20</dt><dd>0</dd><dt>F21</dt><dd>0</dd><dt>F23</dt><dd>0</dd><dt>F24</dt><dd>0</dd><dt>F26</dt><dd>0</dd><dt>F28</dt><dd>0</dd><dt>F29</dt><dd>0</dd><dt>F30</dt><dd>1</dd><dt>F36</dt><dd>1</dd><dt>F38</dt><dd>1</dd></dl>\n",
       "\n",
       "<details>\n",
       "\t<summary style=display:list-item;cursor:pointer>\n",
       "\t\t<strong>Levels</strong>:\n",
       "\t</summary>\n",
       "\t<style>\n",
       "\t.list-inline {list-style: none; margin:0; padding: 0}\n",
       "\t.list-inline>li {display: inline-block}\n",
       "\t.list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "\t</style>\n",
       "\t<ol class=list-inline><li>'0'</li><li>'1'</li></ol>\n",
       "</details>"
      ],
      "text/latex": [
       "\\begin{description*}\n",
       "\\item[F1] 1\n",
       "\\item[F2] 0\n",
       "\\item[F3] 1\n",
       "\\item[F4] 0\n",
       "\\item[F5] 0\n",
       "\\item[F6] 0\n",
       "\\item[F7] 0\n",
       "\\item[F8] 0\n",
       "\\item[F9] 0\n",
       "\\item[F10] 0\n",
       "\\item[F11] 0\n",
       "\\item[F13] 1\n",
       "\\item[F14] 0\n",
       "\\item[F15] 1\n",
       "\\item[F16] 1\n",
       "\\item[F17] 1\n",
       "\\item[F18] 1\n",
       "\\item[F19] 0\n",
       "\\item[F20] 0\n",
       "\\item[F21] 0\n",
       "\\item[F23] 0\n",
       "\\item[F24] 0\n",
       "\\item[F26] 0\n",
       "\\item[F28] 0\n",
       "\\item[F29] 0\n",
       "\\item[F30] 1\n",
       "\\item[F36] 1\n",
       "\\item[F38] 1\n",
       "\\end{description*}\n",
       "\n",
       "\\emph{Levels}: \\begin{enumerate*}\n",
       "\\item '0'\n",
       "\\item '1'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "F1\n",
       ":   1F2\n",
       ":   0F3\n",
       ":   1F4\n",
       ":   0F5\n",
       ":   0F6\n",
       ":   0F7\n",
       ":   0F8\n",
       ":   0F9\n",
       ":   0F10\n",
       ":   0F11\n",
       ":   0F13\n",
       ":   1F14\n",
       ":   0F15\n",
       ":   1F16\n",
       ":   1F17\n",
       ":   1F18\n",
       ":   1F19\n",
       ":   0F20\n",
       ":   0F21\n",
       ":   0F23\n",
       ":   0F24\n",
       ":   0F26\n",
       ":   0F28\n",
       ":   0F29\n",
       ":   0F30\n",
       ":   1F36\n",
       ":   1F38\n",
       ":   1\n",
       "\n",
       "\n",
       "**Levels**: 1. '0'\n",
       "2. '1'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       " F1  F2  F3  F4  F5  F6  F7  F8  F9 F10 F11 F13 F14 F15 F16 F17 F18 F19 F20 F21 \n",
       "  1   0   1   0   0   0   0   0   0   0   0   1   0   1   1   1   1   0   0   0 \n",
       "F23 F24 F26 F28 F29 F30 F36 F38 \n",
       "  0   0   0   0   0   1   1   1 \n",
       "Levels: 0 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "03247797-0780-46bf-ad43-887ec5721cce",
   "metadata": {},
   "outputs": [],
   "source": [
    "trait_df<-data.frame(t(sapply(conditions,function(condition){table(condition)})),row.names=rownames(anno))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e08388e9-2a82-491e-84a4-68943841392d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# This creates an object called \"datExpr\" that contains the normalized counts file output from DESeq2\n",
    "datExpr = vst_mat\n",
    "datExpr = t(datExpr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "74af1c1d-56dc-4be8-b97d-84b99edb767e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       ".list-inline {list-style: none; margin:0; padding: 0}\n",
       ".list-inline>li {display: inline-block}\n",
       ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
       "</style>\n",
       "<ol class=list-inline><li>'X0'</li><li>'X1'</li></ol>\n"
      ],
      "text/latex": [
       "\\begin{enumerate*}\n",
       "\\item 'X0'\n",
       "\\item 'X1'\n",
       "\\end{enumerate*}\n"
      ],
      "text/markdown": [
       "1. 'X0'\n",
       "2. 'X1'\n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "[1] \"X0\" \"X1\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "colnames(trait_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "5df9e23e-71af-445c-9f3e-895e42db8bce",
   "metadata": {},
   "outputs": [],
   "source": [
    "colnames(trait_df) <- c(\"low\", \"high\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "9b5e53ca-523c-45f7-a72e-fb1dbc4c5464",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "TRUE \n",
       "  28 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "datTraits = trait_df\n",
    "table(rownames(datTraits)==rownames(datExpr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e798e0f4-14f5-4507-b14d-4647cc51bb0f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "6d651721-43b0-4453-8269-892f4a8475db",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>png:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{png:} 2"
      ],
      "text/markdown": [
       "**png:** 2"
      ],
      "text/plain": [
       "png \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Cluster samples by expression ----------------------------------------------------------------\n",
    "\n",
    "A = adjacency(t(datExpr),type=\"signed\") # this calculates the whole network connectivity\n",
    "k = as.numeric(apply(A,2,sum))-1 # standardized connectivity\n",
    "Z.k = scale(k)\n",
    "thresholdZ.k = -2.5 # often -2.5\n",
    "outlierColor = ifelse(Z.k<thresholdZ.k,\"red\",\"black\")\n",
    "sampleTree = flashClust(as.dist(1-A), method = \"average\")\n",
    "# Convert traits to a color representation where red indicates high values\n",
    "traitColors = data.frame(numbers2colors((datTraits),signed=FALSE))\n",
    "dimnames(traitColors)[[2]] = paste(names(datTraits))\n",
    "datColors = data.frame(outlier = outlierColor,traitColors)\n",
    "\n",
    "png(filename=\"sampleDendorgramAndTraitHeatmap.png\")\n",
    "plotDendroAndColors(sampleTree,groupLabels=names(datColors),\n",
    "                    colors=datColors,main=\"Sample Dendrogram and Trait Heatmap\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c93a3013-6267-4230-8635-76574f8f6459",
   "metadata": {},
   "source": [
    "**remove outliers**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9b609765-47e5-422b-9dac-00621cd8ab3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "outliers <- c(\"F17\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "e66f000c-f0db-4358-a0da-28ffb38642d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Flagging genes and samples with too many missing values...\n",
      "  ..step 1\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "TRUE"
      ],
      "text/latex": [
       "TRUE"
      ],
      "text/markdown": [
       "TRUE"
      ],
      "text/plain": [
       "[1] TRUE"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gsg = goodSamplesGenes(datExpr, verbose = 3);\n",
    "gsg$allOK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "54a45fb6-7fda-4cc2-ac47-c2adfd9036ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "datTraits <- datTraits[!(row.names(datTraits) %in% outliers),]\n",
    "datExpr <- datExpr[!(row.names(datExpr) %in% outliers),]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "cd9ad641-41d4-481c-8f81-5f4396119a17",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "TRUE \n",
       "  27 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table(rownames(datTraits)==rownames(datExpr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1e110c5a-4f36-496c-a387-2b5de2826e5f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pickSoftThreshold: will use block size 4716.\n",
      " pickSoftThreshold: calculating connectivity for given powers...\n",
      "   ..working on genes 1 through 4716 of 4716\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“executing %dopar% sequentially: no parallel backend registered”\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Power SFT.R.sq    slope truncated.R.sq mean.k. median.k. max.k.\n",
      "1     1.0 5.06e-01  5.16000          0.892  2670.0  2.67e+03   3210\n",
      "2     2.0 1.66e-01  1.37000          0.590  1610.0  1.58e+03   2310\n",
      "3     3.0 5.10e-07  0.00144          0.385  1020.0  9.69e+02   1740\n",
      "4     4.0 1.55e-01 -0.65400          0.231   674.0  6.14e+02   1370\n",
      "5     5.0 3.73e-01 -1.09000          0.262   466.0  3.99e+02   1120\n",
      "6     6.0 4.72e-01 -1.37000          0.323   334.0  2.66e+02    940\n",
      "7     7.0 5.72e-01 -1.54000          0.450   249.0  1.82e+02    813\n",
      "8     8.0 5.73e-01 -1.60000          0.460   191.0  1.26e+02    721\n",
      "9     9.0 6.42e-01 -1.57000          0.570   151.0  8.87e+01    652\n",
      "10   10.0 6.07e-01 -1.59000          0.561   123.0  6.32e+01    599\n",
      "11   10.0 6.07e-01 -1.59000          0.561   123.0  6.32e+01    599\n",
      "12   10.5 5.83e-01 -1.59000          0.552   112.0  5.38e+01    577\n",
      "13   11.0 5.52e-01 -1.59000          0.535   103.0  4.57e+01    558\n",
      "14   11.5 5.28e-01 -1.57000          0.537    94.7  3.92e+01    541\n",
      "15   12.0 5.92e-01 -1.51000          0.604    87.7  3.36e+01    526\n",
      "16   12.5 5.81e-01 -1.52000          0.599    81.8  2.89e+01    513\n",
      "17   13.0 5.56e-01 -1.47000          0.604    76.6  2.49e+01    501\n",
      "18   13.5 5.40e-01 -1.46000          0.595    72.0  2.15e+01    490\n",
      "19   14.0 5.20e-01 -1.42000          0.602    68.1  1.86e+01    480\n",
      "20   14.5 4.91e-01 -1.39000          0.593    64.6  1.61e+01    472\n",
      "21   15.0 4.63e-01 -1.36000          0.590    61.5  1.40e+01    464\n",
      "22   15.5 5.17e-01 -1.32000          0.644    58.8  1.23e+01    457\n",
      "23   16.0 5.21e-01 -1.30000          0.643    56.4  1.07e+01    450\n",
      "24   16.5 5.00e-01 -1.28000          0.633    54.2  9.39e+00    444\n",
      "25   17.0 4.86e-01 -1.25000          0.636    52.3  8.26e+00    439\n",
      "26   17.5 4.85e-01 -1.21000          0.650    50.6  7.29e+00    434\n",
      "27   18.0 4.62e-01 -1.19000          0.644    49.0  6.48e+00    429\n",
      "28   18.5 4.39e-01 -1.17000          0.644    47.6  5.74e+00    425\n",
      "29   19.0 4.15e-01 -1.17000          0.594    46.4  5.08e+00    421\n",
      "30   19.5 4.00e-01 -1.13000          0.605    45.2  4.53e+00    418\n",
      "31   20.0 3.88e-01 -1.11000          0.611    44.2  4.01e+00    414\n",
      "32   20.5 3.55e-01 -1.11000          0.583    43.3  3.58e+00    411\n",
      "33   21.0 3.96e-01 -1.11000          0.611    42.4  3.20e+00    409\n",
      "34   21.5 4.75e-01 -1.08000          0.609    41.6  2.86e+00    406\n",
      "35   22.0 4.70e-01 -1.08000          0.617    40.9  2.57e+00    404\n",
      "36   22.5 4.59e-01 -1.03000          0.639    40.2  2.30e+00    401\n",
      "37   23.0 4.69e-01 -1.02000          0.637    39.6  2.06e+00    399\n",
      "38   23.5 4.40e-01 -0.99700          0.634    39.1  1.86e+00    397\n",
      "39   24.0 4.26e-01 -0.98100          0.618    38.5  1.67e+00    395\n",
      "40   24.5 4.02e-01 -0.96200          0.620    38.1  1.51e+00    393\n",
      "41   25.0 3.98e-01 -0.94600          0.618    37.6  1.37e+00    392\n",
      "42   25.5 3.87e-01 -0.92000          0.621    37.2  1.24e+00    390\n",
      "43   26.0 3.83e-01 -0.90400          0.619    36.8  1.12e+00    389\n",
      "44   26.5 3.82e-01 -0.89400          0.621    36.4  1.01e+00    387\n",
      "45   27.0 3.69e-01 -0.87600          0.622    36.1  9.15e-01    386\n",
      "46   27.5 3.36e-01 -0.87700          0.600    35.8  8.29e-01    385\n",
      "47   28.0 3.26e-01 -0.86300          0.605    35.5  7.53e-01    384\n",
      "48   28.5 3.65e-01 -0.86800          0.649    35.2  6.84e-01    382\n",
      "49   29.0 3.65e-01 -0.86000          0.640    35.0  6.23e-01    381\n",
      "50   29.5 3.89e-01 -0.85500          0.644    34.7  5.69e-01    380\n",
      "51   30.0 3.85e-01 -0.83800          0.645    34.5  5.17e-01    379\n",
      "52   30.5 3.84e-01 -0.82700          0.645    34.2  4.70e-01    378\n",
      "53   31.0 3.84e-01 -0.83500          0.636    34.0  4.30e-01    377\n",
      "54   31.5 3.75e-01 -0.82400          0.642    33.8  3.93e-01    377\n",
      "55   32.0 3.74e-01 -0.81700          0.628    33.6  3.59e-01    376\n",
      "56   32.5 3.72e-01 -0.81700          0.585    33.5  3.29e-01    375\n",
      "57   33.0 3.58e-01 -0.80900          0.580    33.3  2.98e-01    374\n",
      "58   33.5 3.57e-01 -0.80000          0.584    33.1  2.74e-01    373\n",
      "59   34.0 3.54e-01 -0.80200          0.549    33.0  2.52e-01    373\n",
      "60   34.5 3.47e-01 -0.78200          0.568    32.8  2.32e-01    372\n",
      "61   35.0 3.32e-01 -0.75700          0.599    32.7  2.13e-01    371\n",
      "62   35.5 3.30e-01 -0.74600          0.606    32.5  1.96e-01    371\n",
      "63   36.0 3.28e-01 -0.73500          0.603    32.4  1.80e-01    370\n",
      "64   36.5 3.26e-01 -0.72700          0.603    32.3  1.66e-01    369\n",
      "65   37.0 3.49e-01 -0.73300          0.604    32.2  1.54e-01    369\n",
      "66   37.5 3.96e-01 -0.76100          0.623    32.0  1.42e-01    368\n",
      "67   38.0 4.26e-01 -0.75100          0.626    31.9  1.31e-01    368\n",
      "68   38.5 4.48e-01 -0.74600          0.609    31.8  1.20e-01    367\n",
      "69   39.0 4.35e-01 -0.74500          0.580    31.7  1.11e-01    367\n",
      "70   39.5 4.73e-01 -0.75400          0.605    31.6  1.02e-01    366\n",
      "71   40.0 4.71e-01 -0.74600          0.620    31.5  9.48e-02    366\n",
      "72   40.5 4.67e-01 -0.74500          0.612    31.4  8.77e-02    365\n",
      "73   41.0 4.66e-01 -0.73800          0.613    31.3  8.14e-02    364\n",
      "74   41.5 4.65e-01 -0.72800          0.617    31.2  7.56e-02    364\n",
      "75   42.0 4.65e-01 -0.72800          0.601    31.1  7.02e-02    364\n",
      "76   42.5 4.62e-01 -0.72600          0.587    31.0  6.53e-02    363\n",
      "77   43.0 4.60e-01 -0.72000          0.589    31.0  6.05e-02    363\n",
      "78   43.5 4.59e-01 -0.71500          0.590    30.9  5.60e-02    362\n",
      "79   44.0 4.57e-01 -0.71000          0.590    30.8  5.18e-02    362\n",
      "80   44.5 4.41e-01 -0.71400          0.578    30.7  4.81e-02    361\n",
      "81   45.0 4.38e-01 -0.70900          0.578    30.6  4.44e-02    361\n",
      "82   45.5 4.36e-01 -0.70200          0.580    30.6  4.16e-02    361\n",
      "83   46.0 4.33e-01 -0.69600          0.583    30.5  3.86e-02    360\n",
      "84   46.5 4.31e-01 -0.68800          0.586    30.4  3.59e-02    360\n",
      "85   47.0 4.26e-01 -0.68500          0.570    30.4  3.35e-02    359\n",
      "86   47.5 4.19e-01 -0.68700          0.542    30.3  3.12e-02    359\n",
      "87   48.0 4.17e-01 -0.68500          0.542    30.2  2.90e-02    359\n",
      "88   48.5 4.16e-01 -0.67800          0.541    30.1  2.70e-02    358\n",
      "89   49.0 4.20e-01 -0.67700          0.520    30.1  2.52e-02    358\n",
      "90   49.5 4.34e-01 -0.65700          0.564    30.0  2.36e-02    357\n",
      "91   50.0 4.32e-01 -0.65300          0.575    30.0  2.22e-02    357\n",
      "92   50.5 4.77e-01 -0.66400          0.607    29.9  2.06e-02    357\n",
      "93   51.0 4.74e-01 -0.66000          0.605    29.8  1.93e-02    356\n",
      "94   51.5 4.72e-01 -0.65500          0.608    29.8  1.81e-02    356\n",
      "95   52.0 4.71e-01 -0.64900          0.617    29.7  1.69e-02    356\n",
      "96   52.5 4.68e-01 -0.64500          0.610    29.7  1.57e-02    355\n",
      "97   53.0 4.67e-01 -0.64200          0.612    29.6  1.47e-02    355\n",
      "98   53.5 4.46e-01 -0.64600          0.605    29.6  1.36e-02    355\n",
      "99   54.0 4.45e-01 -0.64400          0.605    29.5  1.27e-02    354\n",
      "100  54.5 4.44e-01 -0.64200          0.603    29.4  1.19e-02    354\n",
      "101  55.0 4.42e-01 -0.64000          0.601    29.4  1.11e-02    354\n",
      "102  55.5 4.42e-01 -0.63500          0.602    29.3  1.03e-02    353\n",
      "103  56.0 4.43e-01 -0.63300          0.592    29.3  9.64e-03    353\n",
      "104  56.5 4.15e-01 -0.64300          0.580    29.2  8.98e-03    353\n",
      "105  57.0 4.13e-01 -0.63800          0.582    29.2  8.42e-03    353\n",
      "106  57.5 4.11e-01 -0.63600          0.585    29.1  7.86e-03    352\n",
      "107  58.0 4.10e-01 -0.63200          0.586    29.1  7.38e-03    352\n",
      "108  58.5 4.10e-01 -0.62800          0.587    29.0  6.90e-03    352\n",
      "109  59.0 4.09e-01 -0.62500          0.589    29.0  6.41e-03    351\n",
      "110  59.5 4.09e-01 -0.62300          0.588    29.0  5.99e-03    351\n",
      "111  60.0 4.07e-01 -0.62000          0.592    28.9  5.62e-03    351\n",
      "112  60.5 4.06e-01 -0.61800          0.595    28.9  5.27e-03    350\n",
      "113  61.0 4.04e-01 -0.61500          0.597    28.8  4.92e-03    350\n",
      "114  61.5 4.01e-01 -0.61000          0.601    28.8  4.63e-03    350\n",
      "115  62.0 4.58e-01 -0.62400          0.627    28.7  4.36e-03    350\n",
      "116  62.5 4.57e-01 -0.62100          0.625    28.7  4.08e-03    349\n",
      "117  63.0 4.54e-01 -0.62300          0.610    28.6  3.84e-03    349\n",
      "118  63.5 4.51e-01 -0.62200          0.615    28.6  3.63e-03    349\n",
      "119  64.0 4.45e-01 -0.62500          0.591    28.6  3.40e-03    348\n",
      "120  64.5 4.44e-01 -0.62300          0.588    28.5  3.20e-03    348\n",
      "121  65.0 4.43e-01 -0.62100          0.591    28.5  3.01e-03    348\n",
      "122  65.5 4.42e-01 -0.61900          0.592    28.4  2.83e-03    348\n",
      "123  66.0 4.41e-01 -0.61600          0.593    28.4  2.65e-03    347\n",
      "124  66.5 4.39e-01 -0.61500          0.590    28.4  2.48e-03    347\n",
      "125  67.0 4.37e-01 -0.61100          0.591    28.3  2.32e-03    347\n",
      "126  67.5 4.71e-01 -0.61600          0.583    28.3  2.18e-03    347\n",
      "127  68.0 4.69e-01 -0.61200          0.584    28.2  2.05e-03    346\n",
      "128  68.5 4.79e-01 -0.60700          0.596    28.2  1.92e-03    346\n",
      "129  69.0 4.87e-01 -0.59900          0.612    28.2  1.80e-03    346\n",
      "130  69.5 4.87e-01 -0.59600          0.614    28.1  1.69e-03    346\n",
      "131  70.0 4.85e-01 -0.59200          0.617    28.1  1.60e-03    345\n",
      "132  70.5 4.88e-01 -0.58700          0.624    28.0  1.51e-03    345\n",
      "133  71.0 4.88e-01 -0.58500          0.624    28.0  1.42e-03    345\n",
      "134  71.5 4.87e-01 -0.58300          0.625    28.0  1.33e-03    345\n",
      "135  72.0 4.87e-01 -0.58000          0.628    27.9  1.25e-03    344\n",
      "136  72.5 4.63e-01 -0.59600          0.621    27.9  1.17e-03    344\n",
      "137  73.0 4.63e-01 -0.59400          0.622    27.9  1.09e-03    344\n",
      "138  73.5 4.62e-01 -0.59000          0.624    27.8  1.03e-03    344\n",
      "139  74.0 4.56e-01 -0.59000          0.613    27.8  9.64e-04    344\n",
      "140  74.5 4.55e-01 -0.58800          0.614    27.8  9.10e-04    343\n",
      "141  75.0 4.54e-01 -0.58500          0.614    27.7  8.51e-04    343\n",
      "142  75.5 4.61e-01 -0.58200          0.613    27.7  8.02e-04    343\n",
      "143  76.0 4.60e-01 -0.58000          0.614    27.7  7.52e-04    343\n",
      "144  76.5 4.60e-01 -0.57800          0.614    27.6  7.09e-04    342\n",
      "145  77.0 4.60e-01 -0.57400          0.617    27.6  6.71e-04    342\n",
      "146  77.5 4.59e-01 -0.57200          0.620    27.6  6.33e-04    342\n",
      "147  78.0 4.56e-01 -0.56900          0.623    27.5  5.93e-04    342\n",
      "148  78.5 4.54e-01 -0.56800          0.625    27.5  5.61e-04    341\n",
      "149  79.0 4.55e-01 -0.56500          0.624    27.5  5.30e-04    341\n",
      "150  79.5 4.55e-01 -0.56300          0.621    27.4  5.01e-04    341\n",
      "151  80.0 4.54e-01 -0.56100          0.622    27.4  4.72e-04    341\n",
      "152  80.5 4.53e-01 -0.56000          0.622    27.4  4.47e-04    341\n",
      "153  81.0 4.52e-01 -0.55800          0.623    27.3  4.22e-04    340\n",
      "154  81.5 4.52e-01 -0.55600          0.620    27.3  3.97e-04    340\n",
      "155  82.0 4.86e-01 -0.56300          0.632    27.3  3.75e-04    340\n",
      "156  82.5 5.38e-01 -0.56800          0.644    27.2  3.55e-04    340\n",
      "157  83.0 5.43e-01 -0.56600          0.643    27.2  3.35e-04    340\n",
      "158  83.5 5.53e-01 -0.56400          0.646    27.2  3.15e-04    339\n",
      "159  84.0 5.46e-01 -0.56500          0.632    27.1  2.98e-04    339\n",
      "160  84.5 5.52e-01 -0.56300          0.631    27.1  2.81e-04    339\n",
      "161  85.0 5.53e-01 -0.56200          0.628    27.1  2.64e-04    339\n",
      "162  85.5 5.35e-01 -0.56600          0.602    27.1  2.48e-04    339\n",
      "163  86.0 5.37e-01 -0.56400          0.602    27.0  2.33e-04    338\n",
      "164  86.5 5.34e-01 -0.56300          0.602    27.0  2.20e-04    338\n",
      "165  87.0 5.35e-01 -0.56200          0.601    27.0  2.08e-04    338\n",
      "166  87.5 5.31e-01 -0.56200          0.600    26.9  1.96e-04    338\n",
      "167  88.0 5.32e-01 -0.56000          0.599    26.9  1.85e-04    337\n",
      "168  88.5 5.31e-01 -0.55900          0.600    26.9  1.75e-04    337\n",
      "169  89.0 5.33e-01 -0.55700          0.601    26.8  1.66e-04    337\n",
      "170  89.5 5.29e-01 -0.57400          0.606    26.8  1.56e-04    337\n",
      "171  90.0 4.99e-01 -0.58100          0.555    26.8  1.47e-04    337\n",
      "172  90.5 4.99e-01 -0.57900          0.556    26.8  1.38e-04    336\n",
      "173  91.0 4.96e-01 -0.57900          0.556    26.7  1.30e-04    336\n",
      "174  91.5 4.96e-01 -0.57700          0.557    26.7  1.23e-04    336\n",
      "175  92.0 5.36e-01 -0.57200          0.596    26.7  1.16e-04    336\n",
      "176  92.5 5.36e-01 -0.57000          0.597    26.7  1.10e-04    336\n",
      "177  93.0 5.36e-01 -0.56800          0.598    26.6  1.04e-04    336\n",
      "178  93.5 5.36e-01 -0.56800          0.596    26.6  9.82e-05    335\n",
      "179  94.0 5.65e-01 -0.56300          0.629    26.6  9.28e-05    335\n",
      "180  94.5 5.66e-01 -0.56200          0.628    26.5  8.78e-05    335\n",
      "181  95.0 5.65e-01 -0.56000          0.629    26.5  8.30e-05    335\n",
      "182  95.5 5.65e-01 -0.55900          0.630    26.5  7.85e-05    335\n",
      "183  96.0 5.66e-01 -0.55700          0.631    26.5  7.45e-05    334\n",
      "184  96.5 5.65e-01 -0.57400          0.633    26.4  7.02e-05    334\n",
      "185  97.0 5.69e-01 -0.57100          0.638    26.4  6.63e-05    334\n",
      "186  97.5 5.69e-01 -0.57000          0.639    26.4  6.26e-05    334\n",
      "187  98.0 5.69e-01 -0.56900          0.640    26.4  5.91e-05    334\n",
      "188  98.5 5.69e-01 -0.56800          0.638    26.3  5.57e-05    333\n",
      "189  99.0 5.74e-01 -0.56600          0.642    26.3  5.25e-05    333\n",
      "190  99.5 5.73e-01 -0.56500          0.641    26.3  4.96e-05    333\n",
      "191 100.0 5.64e-01 -0.56600          0.633    26.3  4.70e-05    333\n"
     ]
    }
   ],
   "source": [
    "# Choose a soft threshold power\n",
    "powers = c(c(1:10), seq(from =10, to=100, by=0.5)) #choosing a set of soft-thresholding powers\n",
    "sft = pickSoftThreshold(datExpr, dataIsExpr = TRUE, powerVector=powers, \n",
    "                        verbose =5, networkType=\"signed\",\n",
    "                        corOptions = list(use = 'p')) #call network topology analysis function\n",
    "\n",
    "sizeGrWindow(9,5)\n",
    "par(mfrow= c(1,2))\n",
    "cex1=0.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "434d1566-656a-4a2d-8f1d-90b0d7c04b60",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>png:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{png:} 2"
      ],
      "text/markdown": [
       "**png:** 2"
      ],
      "text/plain": [
       "png \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<strong>png:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{png:} 2"
      ],
      "text/markdown": [
       "**png:** 2"
      ],
      "text/plain": [
       "png \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "png(\"scaleIndependence.png\")\n",
    "plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab= \"Soft Threshold (power)\", \n",
    "     ylab=\"Scale Free Topology Model Fit, signed R^2\", type= \"n\", main= paste(\"Scale independence\"))\n",
    "text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers, cex=cex1, col=\"red\")\n",
    "abline(h=0.80, col=\"red\")\n",
    "dev.off()\n",
    "\n",
    "png(\"meanConnectivity.png\")\n",
    "plot(sft$fitIndices[,1], sft$fitIndices[,5], 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": 22,
   "id": "1a110c36-aa9a-4b2b-8607-df1e6485fe6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "softPower <- 9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6098bdb5-82ce-47eb-833e-821abbe52cf6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Allowing parallel execution with up to 50 working processes.\n"
     ]
    }
   ],
   "source": [
    "nThreads <- 50\n",
    "enableWGCNAThreads(nThreads=nThreads)\n",
    "adjacency = adjacency(datExpr, power = softPower, type = \"signed\") #specify network type\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "b1fcec27-24b3-48af-9e93-1f835b2756fe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..connectivity..\n",
      "..matrix multiplication (system BLAS)..\n",
      "..normalization..\n",
      "..done.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>png:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{png:} 2"
      ],
      "text/markdown": [
       "**png:** 2"
      ],
      "text/plain": [
       "png \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Construct Networks- USE A SUPERCOMPUTER IRL -----------------------------\n",
    "  \n",
    "#translate the adjacency into topological overlap matrix and calculate the corresponding dissimilarity:\n",
    "TOM = TOMsimilarity(adjacency, TOMType=\"signed\") # specify network type\n",
    "dissTOM = 1-TOM\n",
    "  \n",
    "# Generate Modules --------------------------------------------------------\n",
    "  \n",
    "# Generate a clustered gene tree\n",
    "geneTree = flashClust(as.dist(dissTOM), method=\"average\")\n",
    "\n",
    "png(\"geneClusteringTOM.png\")\n",
    "plot(geneTree, xlab=\"\", sub=\"\", main= \"Gene Clustering on TOM-based dissimilarity\", labels= FALSE, hang=0.04)\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "4e9a86ce-8255-4776-bfb6-c547e4ecb031",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] ‘4.2.0’"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "packageVersion(\"clusterProfiler\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "15761e31-4a2c-45a9-96fa-d1f2cf0f87a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " ..cutHeight not given, setting it to 0.983  ===>  99% of the (truncated) height range in dendro.\n",
      " ..done.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>png:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{png:} 2"
      ],
      "text/markdown": [
       "**png:** 2"
      ],
      "text/plain": [
       "png \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#This sets the minimum number of genes to cluster into a module\n",
    "minModuleSize = 10\n",
    "x = 4 \n",
    "dynamicMods = cutreeDynamic(dendro = geneTree, distM = as.matrix(dissTOM), \n",
    "                            method=\"hybrid\", pamStage = F, deepSplit = x, \n",
    "                            minClusterSize = minModuleSize)\n",
    "  \n",
    "dynamicColors= labels2colors(dynamicMods)\n",
    "MEList= moduleEigengenes(datExpr, colors= dynamicColors)#,softPower = 14)\n",
    "MEs= MEList$eigengenes\n",
    "MEDiss= 1-cor(MEs)\n",
    "METree= flashClust(as.dist(MEDiss), method= \"average\")\n",
    "  \n",
    "save(dynamicMods, MEList, MEs, MEDiss, METree, file= \"Network_allSamples_signed_RLDfiltered.RData\")\n",
    "  \n",
    "#plots tree showing how the eigengenes cluster together\n",
    "png(paste(\"eigenGenesClustering\", \".png\", sep=\"\"))\n",
    "plot(METree, main= \"Clustering of module eigengenes\", xlab= \"\", sub= \"\")\n",
    "dev.off()\n",
    "\n",
    "#softPower = 10 # see scale independence file \"softthreshold_SOD1.png\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "1daef8da-b1a1-43a3-b949-389d7f6ad394",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " mergeCloseModules: Merging modules whose distance is less than 0.4\n",
      "   multiSetMEs: Calculating module MEs.\n",
      "     Working on set 1 ...\n",
      "     moduleEigengenes: Calculating 30 module eigengenes in given set.\n",
      "   multiSetMEs: Calculating module MEs.\n",
      "     Working on set 1 ...\n",
      "     moduleEigengenes: Calculating 20 module eigengenes in given set.\n",
      "   multiSetMEs: Calculating module MEs.\n",
      "     Working on set 1 ...\n",
      "     moduleEigengenes: Calculating 19 module eigengenes in given set.\n",
      "   Calculating new MEs...\n",
      "   multiSetMEs: Calculating module MEs.\n",
      "     Working on set 1 ...\n",
      "     moduleEigengenes: Calculating 19 module eigengenes in given set.\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<strong>png:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{png:} 2"
      ],
      "text/markdown": [
       "**png:** 2"
      ],
      "text/plain": [
       "png \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#### -----------------------------------------------------------------------\n",
    "#set a threhold for merging modules. In this example we are not merging so MEDissThres=0.0\n",
    "MEDissThres = 0.4\n",
    "merge = mergeCloseModules(datExpr, dynamicColors, cutHeight= MEDissThres, verbose =3)\n",
    "mergedColors = merge$colors\n",
    "mergedMEs = merge$newMEs\n",
    "\n",
    "#plot dendrogram with module colors below it\n",
    "png(\"geneTree.png\")\n",
    "plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c(\"Dynamic Tree Cut\", \"Merged dynamic\"), dendroLabels= FALSE, hang=0.03, addGuide= TRUE, guideHang=0.05)\n",
    "dev.off()\n",
    "moduleColors = mergedColors\n",
    "colorOrder = c(\"grey\", standardColors(50))\n",
    "moduleLabels = match(moduleColors, colorOrder)-1\n",
    "MEs = mergedMEs\n",
    "\n",
    "save(MEs, moduleLabels, moduleColors, geneTree, file= \"Network_allSamples_signed_nomerge_RLDfiltered.RData\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "5a04e3a1-3416-46e1-b8bf-831ea0e46736",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define number of genes and samples\n",
    "nGenes = ncol(datExpr)\n",
    "nSamples = nrow(datExpr)\n",
    "\n",
    "#Recalculate MEs with color labels\n",
    "MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes\n",
    "MEs = orderMEs(MEs0)\n",
    "moduleTraitCor = cor(MEs, datTraits, use=\"p\")\n",
    "\n",
    "write.csv(moduleTraitCor, \"moduleTraitCor.csv\")\n",
    "\n",
    "moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)\n",
    "\n",
    "write.csv(moduleTraitPvalue, \"moduleTraitPvalue.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "ca9efad1-ddc8-4df3-b0c5-beb1cc4a1371",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<strong>png:</strong> 2"
      ],
      "text/latex": [
       "\\textbf{png:} 2"
      ],
      "text/markdown": [
       "**png:** 2"
      ],
      "text/plain": [
       "png \n",
       "  2 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#display the corelation values with a heatmap plot\n",
    "par(mar=  c(6, 9, 3, 3))\n",
    "png(\"moduleTraitRelationshipHeatmap_corrAll.png\", width=600, height=600)\n",
    "labeledHeatmap(Matrix= t(moduleTraitCor), \n",
    "               yLabels= names(datTraits), \n",
    "               xLabels= names(MEs),\n",
    "               xSymbols= names(MEs),\n",
    "               colorLabels= FALSE, \n",
    "               colors= blueWhiteRed(50), \n",
    "               #textMatrix= textMatrix[module_selection,], \n",
    "               setStdMargins= TRUE, \n",
    "               cex.text= 0.5, \n",
    "               zlim= c(-1,1), \n",
    "               main= paste(\"Module-trait relationships\"))\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "57f56867-e796-4fd2-af07-49043d9b511d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "'getOption(\"repos\")' replaces Bioconductor standard repositories, see\n",
      "'?repositories' for details\n",
      "\n",
      "replacement repositories:\n",
      "    CRAN: https://packagemanager.rstudio.com/all/__linux__/focal/2021-10-18+Y3JhbjoyMDIxLTEwLTE1LDI6NDUyNjIxNTsxMzE0RjdBNg\n",
      "\n",
      "\n",
      "Bioconductor version 3.14 (BiocManager 1.30.16), R 4.1.1 (2021-08-10)\n",
      "\n",
      "Warning message:\n",
      "“package(s) not installed when version(s) same as current; use `force = TRUE` to\n",
      "  re-install: 'pathview'”\n",
      "Installation paths not writeable, unable to update packages\n",
      "  path: /usr/local/lib/R/library\n",
      "  packages:\n",
      "    lattice, mgcv, nlme, survival\n",
      "\n",
      "Old packages: 'basilisk', 'basilisk.utils', 'BiocFileCache', 'BiocParallel',\n",
      "  'biomaRt', 'ChIPseeker', 'clusterProfiler', 'cytolib', 'dir.expiry',\n",
      "  'enrichplot', 'ensembldb', 'GeneTonic', 'GenomeInfoDb', 'GenomicFeatures',\n",
      "  'limma', 'MassSpecWavelet', 'quantiseqr', 'rhdf5', 'S4Vectors', 'singscore'\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#get selected genes\n",
    "selected_Colors<-sapply(names(MEs),function(ME){sub(\"ME\",\"\",ME)})\n",
    "#names(datExpr)[moduleColors==\"brown\"]\n",
    "module_genes<-sapply(selected_Colors,function(module_color){rownames(t(datExpr)[mergedColors==module_color,])})\n",
    "\n",
    "WGCNA_folder<-\"./\"\n",
    "saveRDS(module_genes,file=file.path(WGCNA_folder,'module_genes.rds'))\n",
    "write(toJSON(module_genes),file=file.path(WGCNA_folder,'module_genes.json'))\n",
    "\n",
    "## enrichment of each module\n",
    "if (!requireNamespace(\"BiocManager\", quietly = TRUE))\n",
    "  install.packages(\"BiocManager\")\n",
    "\n",
    "BiocManager::install(\"pathview\")\n",
    "\n",
    "suppressPackageStartupMessages({\n",
    "  library(clusterProfiler)\n",
    "  library(enrichplot)\n",
    "  library(DOSE)\n",
    "  library(org.Mm.eg.db)\n",
    "  library(org.Hs.eg.db)\n",
    "  library(ggplot2)\n",
    "  library(pathview)\n",
    "  library(tidyverse)\n",
    "})\n",
    "\n",
    "make_go_enr = function(genelist, \n",
    "                       fname, \n",
    "                       folder, \n",
    "                       fromType=\"SYMBOL\",\n",
    "                       #OrgDb=org.Mm.eg.db,\n",
    "                       OrgDb=org.Hs.eg.db,\n",
    "                       prefixes=c(\"BP\",\"MF\",\"CC\",\"All\"),\n",
    "                       image_width=30,\n",
    "                      image_height=30){\n",
    "  \n",
    "  gene.df <- bitr((genelist), fromType = fromType,\n",
    "                  toType = c(\"ENTREZID\"),\n",
    "                  OrgDb = OrgDb)\n",
    "  for (prefix in prefixes) {\n",
    "    ego <- enrichGO(gene          = gene.df$ENTREZID,\n",
    "                    OrgDb         = OrgDb,\n",
    "                    ont           = prefix,\n",
    "                    pAdjustMethod = \"BH\",\n",
    "                    pvalueCutoff  = 0.05,\n",
    "                    qvalueCutoff  = 0.05,\n",
    "                    readable = T)\n",
    "    \n",
    "    dir.create(file.path(folder, fname))\n",
    "    \n",
    "    if (!is.null(ego) && nrow(ego)!=0){\n",
    "      write.csv(ego, file.path(folder, fname, paste0(\"go_enr_\", prefix,\"_\",fname,\".csv\")))\n",
    "\n",
    "      p1 = dotplot(ego, showCategory=20) + ggtitle(paste0(\"GO enrichment \", prefix))\n",
    "      ggsave(p1, filename = file.path(folder,fname, paste0(\"dotplot_go_\", prefix,\"_\",fname,\".pdf\")), width = image_width, \n",
    "             height=image_height, units=\"cm\")\n",
    "\n",
    "      p2 = barplot(ego, showCategory=20) + ggtitle(paste0(\"GO enrichment \", prefix))\n",
    "      ggsave(p2, filename = file.path(folder, fname, paste0(\"barplot_go_\", prefix,\"_\",fname,\".pdf\")), width = image_width, \n",
    "             height=image_height, units=\"cm\")\n",
    "    }\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "303940e1-e42b-4c29-a994-6c0c879daa5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "make_go_enr1 = function(genelist, \n",
    "                       fname, \n",
    "                       folder, \n",
    "                       fromType=\"SYMBOL\",\n",
    "                       #OrgDb=org.Mm.eg.db,\n",
    "                       OrgDb=org.Hs.eg.db,\n",
    "                       prefixes=c(\"BP\",\"MF\",\"CC\",\"All\"),\n",
    "                       image_width=20,\n",
    "                        image_height=20){\n",
    "  \n",
    "  gene.df <- bitr((genelist), fromType = fromType,\n",
    "                  toType = c(\"ENTREZID\"),\n",
    "                  OrgDb = OrgDb)\n",
    "  for (prefix in prefixes) {\n",
    "    ego <- enrichGO(gene          = gene.df$ENTREZID,\n",
    "                    OrgDb         = OrgDb,\n",
    "                    ont           = prefix,\n",
    "                    pAdjustMethod = \"BH\",\n",
    "                    pvalueCutoff  = 0.05,\n",
    "                    qvalueCutoff  = 0.05,\n",
    "                    readable = T)\n",
    "    \n",
    "    dir.create(file.path(folder, fname))\n",
    "    \n",
    "    if (!is.null(ego) && nrow(ego)!=0){\n",
    "      write.csv(ego, file.path(folder, fname, paste0(\"go_enr_\", prefix,\"_\",fname,\".csv\")))\n",
    "\n",
    "      p1 = dotplot(ego, showCategory=10, font.size=20) + ggtitle(paste0(\"GO enrichment \", prefix))\n",
    "      ggsave(p1, filename = file.path(folder,fname, paste0(\"dotplot_go_\", prefix,\"_\",fname,\".pdf\")), width = image_width, \n",
    "             height=image_height, units=\"cm\")\n",
    "\n",
    "      p2 = barplot(ego, showCategory=10, font.size=20) + ggtitle(paste0(\"GO enrichment \", prefix))\n",
    "      ggsave(p2, filename = file.path(folder, fname, paste0(\"barplot_go_\", prefix,\"_\",fname,\".pdf\")), width = image_width, \n",
    "             height=image_height, units=\"cm\")\n",
    "    }\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "340113a8-c8ca-44b8-9763-93d62edf7b2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "modules_ <- c(\"MEblue\", \"MEbrown\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "fe0af972-20db-49a8-bceb-345a152be8de",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in dir.create(result_folder):\n",
      "“'results_MEblue_MEbrown' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“1.67% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results_MEblue_MEbrown/MEblue' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results_MEblue_MEbrown/MEblue' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results_MEblue_MEbrown/MEblue' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results_MEblue_MEbrown/MEblue' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“6.02% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results_MEblue_MEbrown/MEbrown' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results_MEblue_MEbrown/MEbrown' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results_MEblue_MEbrown/MEbrown' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results_MEblue_MEbrown/MEbrown' already exists”\n"
     ]
    }
   ],
   "source": [
    "result_folder<-\"results_MEblue_MEbrown\"\n",
    "overview_df<- data.frame()\n",
    "\n",
    "dir.create(result_folder)\n",
    "\n",
    "for (module in modules_){ \n",
    "\n",
    "  make_go_enr1(module_genes[[module]],module,file.path(WGCNA_folder,result_folder))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "ebbaeeb1-6f62-4428-871a-8f9c07f97790",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in dir.create(result_folder):\n",
      "“'results' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkorange' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkorange' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkorange' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkturquoise' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkturquoise' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkturquoise' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“2.22% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkred' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkred' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkred' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“26.67% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MElightcyan' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MElightcyan' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MElightcyan' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“4.55% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEroyalblue' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEroyalblue' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEroyalblue' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“9.82% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEyellow' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEyellow' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEyellow' already exists”\n",
      "'select()' returned 1:many mapping between keys and columns\n",
      "\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkgreen' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkgreen' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkgreen' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEgrey60' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEgrey60' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEgrey60' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEsaddlebrown' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEsaddlebrown' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEsaddlebrown' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEwhite' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEwhite' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEwhite' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“12.07% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEmagenta' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEmagenta' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEmagenta' already exists”\n",
      "'select()' returned 1:many mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“5.75% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkgrey' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkgrey' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEdarkgrey' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“5.41% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEturquoise' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEturquoise' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEturquoise' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“6.02% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEbrown' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEbrown' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEbrown' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“1.82% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEpurple' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEpurple' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEpurple' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“1.67% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEblue' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEblue' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEblue' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“8% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEblack' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEblack' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEblack' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“3.03% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEcyan' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEcyan' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEcyan' already exists”\n",
      "'select()' returned 1:1 mapping between keys and columns\n",
      "\n",
      "Warning message in bitr((genelist), fromType = fromType, toType = c(\"ENTREZID\"), :\n",
      "“4.7% of input gene IDs are fail to map...”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEgrey' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEgrey' already exists”\n",
      "Warning message in dir.create(file.path(folder, fname)):\n",
      "“'.//results/MEgrey' already exists”\n"
     ]
    }
   ],
   "source": [
    "result_folder<-\"results\"\n",
    "overview_df<- data.frame()\n",
    "\n",
    "dir.create(result_folder)\n",
    "\n",
    "for (module in names(module_genes)){ \n",
    "\n",
    "  make_go_enr(module_genes[[module]],module,file.path(WGCNA_folder,result_folder))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42edadb3-1b48-4e3a-92fc-f586f7aa013b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# modules_ <- c(\"MEblue\", \"MEbrown\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "becb8063-0fe3-493d-857d-81cecd3d33f8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# result_folder<-\"results\"\n",
    "# overview_df<- data.frame()\n",
    "# \n",
    "# dir.create(result_folder)\n",
    "# \n",
    "# for (module in modules_){ # names(module_genes)\n",
    "# \n",
    "#   make_go_enr(module_genes[[module]],module,file.path(WGCNA_folder,result_folder))\n",
    "# }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e9c3997-6ee1-489d-995a-5313b177cf90",
   "metadata": {},
   "outputs": [],
   "source": [
    "# modules_ <- c(\"MEroyalblue\", \"MEmagenta\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c86af21c-baff-46a8-a46e-5463319988c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# result_folder<-\"results\"\n",
    "# overview_df<- data.frame()\n",
    "\n",
    "# dir.create(result_folder)\n",
    "\n",
    "# for (module in modules_){ # names(module_genes)\n",
    "\n",
    "#   make_go_enr(module_genes[[module]],module,file.path(WGCNA_folder,result_folder))\n",
    "# }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db3ce386-7522-4016-8ae4-2c4b976cc5a1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": "4.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}