clusterprofiler.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyNqQAR0BSK3UIz3L9Mfk1MK",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "ir",
      "display_name": "R"
    },
    "language_info": {
      "name": "R"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/francescopatane96/RNA-seq_Pipeline/blob/main/clusterprofiler.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "if (!require(\"BiocManager\", quietly = TRUE))\n",
        "    install.packages(\"BiocManager\")\n",
        "\n",
        "BiocManager::install(\"clusterProfiler\")"
      ],
      "metadata": {
        "id": "MWmNQNhywzIi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "o37s4NfnwfNU"
      },
      "outputs": [],
      "source": [
        "library(clusterProfiler)\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "organism = \"org.Dm.eg.db\"\n",
        "BiocManager::install(organism, character.only = TRUE)\n",
        "library(organism, character.only = TRUE)"
      ],
      "metadata": {
        "id": "KUvgnW7ixYX_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "genelistup <- read.csv('your_file.csv', header=F, sep=' ')\n",
        "genelist_UP <- genelistup$V2                                    #select gene column\n",
        "\n",
        "\n",
        "genelist <- bitr(genelist_UP, fromType = \"SYMBOL\",          \n",
        "                toType=c(\"ENTREZID\", \"ENSEMBL\"),\n",
        "                OrgDb=\"org.Hs.eg.db\")\n",
        "\n",
        "gene_list=genelist$ENSEMBL\n",
        "\n",
        "gene_list_UP=as.character(gene_list)"
      ],
      "metadata": {
        "id": "YpWHl0rzw246"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "GO classification.\n",
        "\n",
        "- Functional Profile of a gene set at specific GO level. Given a vector of genes, this function will return the GO profile at a specific level"
      ],
      "metadata": {
        "id": "prIGzwjfxapK"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "ggo <- groupGO(gene=gene_list_UP, \n",
        "               OrgDb=\"org.Hs.eg.db\",\n",
        "               keyType='ENSEMBL',\n",
        "               ont=\"BP\", \n",
        "               level=3,\n",
        "               readable=T)\n",
        "head(ggo)"
      ],
      "metadata": {
        "id": "GlPG16GwxeBh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "barplot(ggo, drop=TRUE, showCategory=12)"
      ],
      "metadata": {
        "id": "tdI1m_nGxliJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "GO ORA (over-representation analysis)"
      ],
      "metadata": {
        "id": "95oDprkIxo1S"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "ego <- enrichGO(gene          = gene_list_UP,\n",
        "                OrgDb         = org.Hs.eg.db,\n",
        "                keyType       = \"ENSEMBL\",\n",
        "                ont           = \"BP\",\n",
        "                pAdjustMethod = \"BH\",\n",
        "                pvalueCutoff  = 0.1,\n",
        "                qvalueCutoff  = 0.1,\n",
        "                readable      = TRUE)\n",
        "head(ego)"
      ],
      "metadata": {
        "id": "SUwu7KAQxqzD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "library(DOSE)\n",
        "library(enrichplot)\n",
        "library(wordcloud)"
      ],
      "metadata": {
        "id": "VRQS2oR6xx4L"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "upsetplot(ego)"
      ],
      "metadata": {
        "id": "MF5qI6bGx0-6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "wcdf<-read.table(text=ego$GeneRatio, sep = \"/\")[1]\n",
        "wcdf$term<-ego[,2]\n",
        "wordcloud(words = wcdf$term, freq = wcdf$V1, scale=(c(4, .1)), colors=brewer.pal(8, \"Dark2\"), max.words = 25)"
      ],
      "metadata": {
        "id": "_HaOGzWGx3Pb"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "barplot(ego, \n",
        "        drop = TRUE, \n",
        "        showCategory = 10, \n",
        "        title = \"GO Biological Pathways\",\n",
        "        font.size = 8)"
      ],
      "metadata": {
        "id": "Z9HVKLPBx5-i"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "dotplot(ego)"
      ],
      "metadata": {
        "id": "L37l4uqFx6Vy"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "emapplot(ego)"
      ],
      "metadata": {
        "id": "QKFiETzUx9Gy"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "goplot(ego, showCategory = 10)"
      ],
      "metadata": {
        "id": "z_om6iwJx-3y"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plotGOgraph(ego)\n"
      ],
      "metadata": {
        "id": "40vD3jq6yAPc"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# categorySize can be either 'pvalue' or 'geneNum'\n",
        "cnetplot(ego, categorySize=\"pvalue\", foldChange=gene_list)"
      ],
      "metadata": {
        "id": "wm-daPjFyCbS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "GO gene set enrichment analysis\n",
        "\n",
        "- geneList input would be a ranked gene list, with 3 features: numeric vector,\n",
        "named vector and sorted vector.\n",
        "if import data from csv, you must have 2 columns,\n",
        "one for gene ID and another one for fold change"
      ],
      "metadata": {
        "id": "ILeWZG_gyDOx"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "d <- read.csv('path.csv', header=T, sep=',')"
      ],
      "metadata": {
        "id": "f76nUL_hyFRJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "## assume 2nd column is ID\n",
        "## 3rd column is FC\n",
        "\n",
        "## feature 1: numeric vector\n",
        "geneList = d[,2]\n",
        "\n",
        "## feature 2: named vector\n",
        "names(geneList) = as.character(d[,1])\n",
        "\n",
        "## feature 3: decreasing order\n",
        "geneList = sort(geneList, decreasing = TRUE)"
      ],
      "metadata": {
        "id": "qZjfrrqjyMWa"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "gse <- gseGO(geneList=geneList, \n",
        "             OrgDb='org.Hs.eg.db',\n",
        "             ont =\"ALL\", \n",
        "             keyType = \"ENSEMBL\", \n",
        "             nPerm = 10000, \n",
        "             minGSSize = 3, \n",
        "             maxGSSize = 800, \n",
        "             pvalueCutoff = 0.1, \n",
        "             verbose = TRUE, \n",
        "             pAdjustMethod = \"none\")\n",
        "\n",
        "head(gse)"
      ],
      "metadata": {
        "id": "reLMqUDmySER"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#Dotplot\n",
        "\n",
        "dotplot(gse, showCategory=15, split=\".sign\") "
      ],
      "metadata": {
        "id": "JX6bIcRNyVLD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#The goplot() function can accept the output of enrichGO and visualize the enriched \n",
        "#GO induced graph.\n",
        "goplot(ego)"
      ],
      "metadata": {
        "id": "v-_NbrHdyXZy"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Use the `Gene Set` param for the index in the title, and as the value for geneSetId\n",
        "gseaplot(gse, by = \"all\", title = gse$Description[1], geneSetID = 1)"
      ],
      "metadata": {
        "id": "CFAgAKPyyZGD"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "KEGG ORA analysis"
      ],
      "metadata": {
        "id": "WNEtgcz-ybfx"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "search_kegg_organism('sap', by='kegg_code')\n",
        "search_kegg_organism('sapiens', by='scientific_name')\n",
        "#hsa\n",
        "\n",
        "#Input ID type can be kegg, ncbi-geneid, ncbi-proteinid or uniprot\n",
        "\n",
        "prova <- bitr(gene_list, fromType = \"ENSEMBL\",\n",
        "                 toType=c(\"UNIPROT\"),\n",
        "                 OrgDb=\"org.Hs.eg.db\")"
      ],
      "metadata": {
        "id": "vLvFVRyWyiWV"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "kegg_type <- bitr_kegg(prova$UNIPROT, fromType='uniprot', toType='kegg', organism='hsa')\n",
        "\n",
        "kk <- enrichKEGG(gene         = kegg_type$kegg,\n",
        "                 organism     = 'hsa',\n",
        "                 pvalueCutoff = 0.05)\n",
        "head(kk)\n"
      ],
      "metadata": {
        "id": "dvjWZxo7yjLd"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "KEGG GSEA"
      ],
      "metadata": {
        "id": "fQ78judbymW6"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#geneList must be a decreasing sorted vector\n",
        "kk2 <- gseKEGG(geneList     = kegg_type$kegg,\n",
        "               organism     = 'hsa',\n",
        "               nPerm        = 1000,\n",
        "               minGSSize    = 120,\n",
        "               pvalueCutoff = 0.05,\n",
        "               verbose      = FALSE)\n",
        "head(kk2)"
      ],
      "metadata": {
        "id": "xhuF21qmyoRh"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "gseaplot(kk2, geneSetID = \"hsa04145\")"
      ],
      "metadata": {
        "id": "f5gbf1sayqeq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "KEGG module ORA"
      ],
      "metadata": {
        "id": "tQfkihAMytfq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "mkk <- enrichMKEGG(gene = kegg_type$kegg,\n",
        "                   organism = 'hsa',\n",
        "                   pvalueCutoff = 1,\n",
        "                   qvalueCutoff = 1)\n",
        "head(mkk)           "
      ],
      "metadata": {
        "id": "K4uI0_5iyvgT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "KEGG module GSEA"
      ],
      "metadata": {
        "id": "y0V6P8HvyxOc"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "mkk2 <- gseMKEGG(geneList = geneList,\n",
        "                 organism = 'hsa',\n",
        "                 pvalueCutoff = 1)\n",
        "head(mkk2)"
      ],
      "metadata": {
        "id": "eFOA_6ieyzBi"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "visualize enriched KEGG pathway"
      ],
      "metadata": {
        "id": "2SkapzFKy1bR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "browseKEGG(kk, 'hsa05205')"
      ],
      "metadata": {
        "id": "025k9nYky0xq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "library(\"pathview\")\n",
        " pathview(gene.data=kegg_type$kegg,\n",
        "                     pathway.id = \"hsa05205\",\n",
        "                     species    = \"hsa\",\n",
        "                     limit      = list(gene=max(abs(geneList)), cpd=1))"
      ],
      "metadata": {
        "id": "2Wer-YNry6Gx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Produce the native KEGG plot (PNG)\n",
        " dme <- pathview(gene.data=kegg_type$kegg, pathway.id=\"hsa05205\", species = 'hsa')\n",
        " "
      ],
      "metadata": {
        "id": "kJNjz5Fny8Pq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Produce a different plot (PDF) \n",
        " dme <- pathview(gene.data=kegg_type$kegg, pathway.id=\"hsa05205\", species = 'hsa', kegg.native = F)"
      ],
      "metadata": {
        "id": "n0WCWl94y-F5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "barplot(kk, \n",
        "         showCategory = 10, \n",
        "         title = \"Enriched Pathways\",\n",
        "         font.size = 8)"
      ],
      "metadata": {
        "id": "4AKfsCmJzAPx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "dotplot(kk, \n",
        "         showCategory = 10, \n",
        "         title = \"Enriched Pathways\",\n",
        "         font.size = 8)"
      ],
      "metadata": {
        "id": "0s_0lGIzzB0x"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        " # categorySize can be either 'pvalue' or 'geneNum'\n",
        " cnetplot(kk, categorySize=\"pvalue\", foldChange=gene_list)"
      ],
      "metadata": {
        "id": "ae270hoPzE6E"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "WikiPathways analysis"
      ],
      "metadata": {
        "id": "cUi2Bo5ozGmx"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#ORA\n",
        " enrichWP(genelist$ENTREZID, organism = \"Homo sapiens\") "
      ],
      "metadata": {
        "id": "j5KqjnxyzHEx"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#GSEA (decreasing sorted vector)\n",
        " gseWP(gene_list, organism = \"Homo sapiens\")"
      ],
      "metadata": {
        "id": "U0vQsvSbzJVK"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "REACTOME"
      ],
      "metadata": {
        "id": "DvtzZ555zLKZ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "library(ReactomePA)"
      ],
      "metadata": {
        "id": "iRBn2snszNFR"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#ORA\n",
        " x <- enrichPathway(gene=genelist$ENTREZID, pvalueCutoff = 0.05, readable=TRUE)\n",
        " head(x)"
      ],
      "metadata": {
        "id": "OYzPXkRKzPHJ"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#GSEA\n",
        " y <- gsePathway(geneList, \n",
        "                 pvalueCutoff = 0.2,\n",
        "                 pAdjustMethod = \"BH\", \n",
        "                 verbose = FALSE)\n",
        " head(y)"
      ],
      "metadata": {
        "id": "HUXtx607zQay"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#pathway visualization\n",
        " \n",
        " \n",
        " viewPathway(\"E2F mediated regulation of DNA replication\", \n",
        "             readable = TRUE, \n",
        "             foldChange = d$V2)"
      ],
      "metadata": {
        "id": "IgXWHZFOzSIR"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "disease enrichment analysis (DO)"
      ],
      "metadata": {
        "id": "hAF_iiaBzT8p"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "library(DOSE)"
      ],
      "metadata": {
        "id": "_-eFCk4GzYhF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#ORA\n",
        " x <- enrichDO(gene          = genelist$ENTREZID,\n",
        "               ont           = \"DO\",\n",
        "               pvalueCutoff  = 0.05,\n",
        "               pAdjustMethod = \"BH\",\n",
        "               universe      = names(genelist$ENTREZID),\n",
        "               minGSSize     = 5,\n",
        "               maxGSSize     = 500,\n",
        "               qvalueCutoff  = 0.05,\n",
        "               readable      = FALSE)\n",
        " head(x)"
      ],
      "metadata": {
        "id": "sM9AoiuCzaBK"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#GSEA\n",
        " gene2 <- names(genelist$ENTREZID)\n",
        " ncg <- enrichNCG(gene2) \n",
        " head(ncg)"
      ],
      "metadata": {
        "id": "IHXWd9pwzdNp"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}