WGCNA_R.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "WGCNA.R",
      "version": "0.3.2",
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/Yixiangzhang1996/2019Bio.Script.Test/blob/master/WGCNA_R.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_-TU1k9tN2kV",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "library(WGCNA)\n",
        "\n",
        "exp0 <- read.delim(\"~/naeem/1.transcription_ana/data/gene_exp_tpm.txt\", row.names=1)\n",
        "\n",
        "#=====================================================================================\n",
        "#\n",
        "#  1.导入表达数据\n",
        "#\n",
        "#=====================================================================================\n",
        "\n",
        "  \n",
        "  # 过滤掉在所有样品中表达之和小于10 的基因\n",
        "   exp <- exp0[rowSums(exp0) > 10, ]\n",
        "  \n",
        "  # 转置\n",
        "  datExpr <- as.data.frame(t(exp0))\n",
        "  dim(datExpr)\n",
        "  rownames(datExpr)\n",
        "\n",
        "\n",
        "\n",
        "#2.寻找最佳 β 值\n",
        "#\n",
        "#=====================================================================================\n",
        "{\n",
        "  # 开启多线程模式\n",
        "  enableWGCNAThreads(nThreads = 20)\n",
        "  \n",
        "  # 通过对 power 的多次迭代，确定最佳 power\n",
        "  powers <- c(1:30)\n",
        "  sft <- pickSoftThreshold(datExpr, \n",
        "                           powerVector = powers, \n",
        "                           verbose = 5,\n",
        "                           networkType = \"signed\"\n",
        "  )\n",
        "  # 计算出来的最佳 β 存放在\n",
        "  sft$powerEstimate\n",
        "  \n",
        "  # 画图\n",
        "  sizeGrWindow(9, 5)\n",
        "  par(mfrow = c(1,2))\n",
        "  cex1 <- 0.9\n",
        "  #  R2 ~ 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.8,col=\"red\")\n",
        "  # Mean connectivity ~ 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",
        "}\n",
        "#=====================================================================================\n",
        "#\n",
        "#  3. 构建网络\n",
        "##  3.1. 计算相关系数\n",
        "##  3.2. 计算邻接矩阵\n",
        "##  3.3. 计算 TOM 矩阵\n",
        "##  3.4. 聚类并划分模块\n",
        "##  3.5. 合并相似模块\n",
        "#=====================================================================================\n",
        "{\n",
        "  net <- blockwiseModules(\n",
        "    # 0.输入数据\n",
        "    datExpr, \n",
        "    \n",
        "    # 1. 计算相关系数\n",
        "    corType = \"pearson\", # 相关系数算法，pearson|bicor\n",
        "    \n",
        "    # 2. 计算邻接矩阵\n",
        "    power = 12, # 前面得到的 soft power\n",
        "    networkType = \"unsigned\", # unsigned | signed | signed hybrid\n",
        "    \n",
        "    # 3. 计算 TOM 矩阵\n",
        "    TOMType = \"unsigned\", # none | unsigned | signed\n",
        "    saveTOMs = TRUE, # 是否保存\n",
        "    saveTOMFileBase = \"blockwiseTOM\", # 保存文件前缀\n",
        "    \n",
        "    # 4. 聚类并划分模块\n",
        "    deepSplit = 2, # 0|1|2|3|4, 值越大得到的模块就越多越小\n",
        "    minModuleSize = 30,\n",
        "    \n",
        "    # 5. 合并相似模块\n",
        "    ## 5.1 计算模块特征向量（module eigengenes， MEs），即第一个主成分（1st PC）\n",
        "    ## 5.2 计算 MEs 与 datTrait 之间的相关性\n",
        "    ## 5.3 对距离小于 mergeCutHeight （1-cor）的模块进行合并\n",
        "    mergeCutHeight = 0.25, \n",
        "    \n",
        "    # 其他参数\n",
        "    numericLabels = TRUE, # 以数字命名模块\n",
        "    nThreads = 0 # 0 则使用所有可用线程\n",
        "  )\n",
        "  # 查看每个模块包含基因数目\n",
        "  table(net$colors) \n",
        "}\n",
        "\n",
        "#  4. 可视化\n",
        "#\n",
        "#=====================================================================================\n",
        "{\n",
        "  # 模块名称修改为颜色\n",
        "  moduleColors <- labels2colors(net$colors)\n",
        "  # 同时绘制聚类图和模块颜色\n",
        "  # pdf(file = \"output/plotDendroAndColors.pdf\")\n",
        "  plotDendroAndColors(\n",
        "    net$dendrograms[[1]], \n",
        "    moduleColors[net$blockGenes[[1]]],\n",
        "    \"Module colors\",\n",
        "    dendroLabels = FALSE, \n",
        "    addGuide = TRUE)\n",
        "  # dev.off()\n",
        "}\n"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}