Phase 2/Bio_Phase_2.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wcsSGbbglOyC"
      },
      "source": [
        "# Phase 2\n",
        "Analysing AML microarray."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-6SspO77lk6d"
      },
      "source": [
        "## 1. Prepairing Data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "piMJAAcYrKSM"
      },
      "source": [
        "### Installing Packages"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pmcs6qBWY-UH"
      },
      "outputs": [],
      "source": [
        "!pip uninstall rpy2 -y\n",
        "!pip install rpy2==3.5.1\n",
        "%load_ext rpy2.ipython"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4uk61UNmZTg9"
      },
      "outputs": [],
      "source": [
        "%%R\n",
        "if (!require(\"BiocManager\", quietly = TRUE))\n",
        "    install.packages(\"BiocManager\", repos='http://cran.us.r-project.org')\n",
        "BiocManager::install()\n",
        "BiocManager::install(\"GEOquery\")\n",
        "\n",
        "install.packages(\"limma\")\n",
        "install.packages(\"umap\")\n",
        "install.packages(\"maptools\")\n",
        "install.packages(\"corrplot\")\n",
        "install.packages(\"pheatmap\")\n",
        "install.packages(\"reshape2\")\n",
        "install.packages(\"mvtnorm\")\n",
        "install.packages(\"Rtsne\")\n",
        "install.packages(\"M3C\")\n",
        "install.packages(\"MASS\")\n",
        "install.packages(\"ggpubr\")\n",
        "install.packages(c(\"vegan\", \"ecodist\", \"labdsv\", \"ape\", \"ade4\", \"smacof\"))\n",
        "install.packages(\"fastICA\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uA58YioPrMud"
      },
      "source": [
        "### Libraries"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QEwShNlirP4h"
      },
      "outputs": [],
      "source": [
        "%%R\n",
        "\n",
        "library(Biobase)\n",
        "library(GEOquery)\n",
        "library(limma)\n",
        "library(umap)\n",
        "library(maptools)\n",
        "library(tidyr)\n",
        "library(dplyr)\n",
        "library(ggplot2)\n",
        "library(devtools)\n",
        "install_github(\"vqv/ggbiplot\")\n",
        "library(ggbiplot)\n",
        "library(data.table)\n",
        "library(corrplot)\n",
        "library(pheatmap)\n",
        "library(reshape2)\n",
        "library(plyr)\n",
        "library(MASS)\n",
        "library(mvtnorm)\n",
        "library(Rtsne)\n",
        "library(\"ggpubr\")\n",
        "library(vegan)\n",
        "library(ecodist)\n",
        "library(labdsv)\n",
        "library(ape)\n",
        "library(ade4)\n",
        "library(fastICA)\n",
        "library(smacof)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "av8SnzRLuu23"
      },
      "source": [
        "### Loading Data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "XLvROrm-l3fg",
        "outputId": "41ba455f-a896-4f43-c4af-2f21e38a9b98"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "WARNING:rpy2.rinterface_lib.callbacks:R[write to console]: Found 1 file(s)\n",
            "\n",
            "WARNING:rpy2.rinterface_lib.callbacks:R[write to console]: GSE48558_series_matrix.txt.gz\n",
            "\n"
          ]
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "$GSE48558_series_matrix.txt.gz\n",
            "ExpressionSet (storageMode: lockedEnvironment)\n",
            "assayData: 32321 features, 170 samples \n",
            "  element names: exprs \n",
            "protocolData: none\n",
            "phenoData\n",
            "  sampleNames: GSM1180750 GSM1180751 ... GSM1180919 (170 total)\n",
            "  varLabels: title geo_accession ... phenotype:ch1 (32 total)\n",
            "  varMetadata: labelDescription\n",
            "featureData\n",
            "  featureNames: 7892501 7892502 ... 8180418 (32321 total)\n",
            "  fvarLabels: ID Gene title ... GO:Component ID (21 total)\n",
            "  fvarMetadata: Column Description labelDescription\n",
            "experimentData: use 'experimentData(object)'\n",
            "  pubMedIds: 23836560 \n",
            "Annotation: GPL6244 \n",
            "\n"
          ]
        }
      ],
      "source": [
        "%%R\n",
        "\n",
        "gds <- getGEO(\"GSE48558\", GSEMatrix = TRUE, AnnotGPL=TRUE)\n",
        "show(gds)\n",
        "gds48558 <- gds[[1]]\n",
        "\n",
        "fvarLabels(gds48558) <- make.names(fvarLabels(gds48558))\n",
        "\n",
        "# group membership for all samples\n",
        "gsms <- paste0(\"1111111111111XXXXXXXXXXXXXXXXXXXXXXXXXXX0XXX0XXXXX\",\n",
        "               \"XXXXXXXXXXXXXXXXXX0X0XXX0X0000X0XX00XX00X0X0X0X0X0\",\n",
        "               \"XXX0XXX0XXXXXXXXXXXXXXXXXXXXXXXXXXXXX0000000110111\",\n",
        "               \"00000000000000000000\")\n",
        "sml <- strsplit(gsms, split=\"\")[[1]]\n",
        "\n",
        "sel <- which(sml != \"X\")\n",
        "sml <- sml[sel]\n",
        "gds48558 <- gds48558[ ,sel]\n",
        "\n",
        "gs <- factor(sml)\n",
        "groups <- make.names(c(\"Healthy\", \"AML\"))\n",
        "levels(gs) <- groups\n",
        "gds48558$group <- gs"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "n_SQ7gRwsmcN",
        "outputId": "4bebc786-0360-4e4f-c73b-169d1a56951c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Features  Samples \n",
            "   32321       67 \n"
          ]
        }
      ],
      "source": [
        "%%R\n",
        "\n",
        "expr <- exprs(gds48558)\n",
        "features <- fData(gds48558)\n",
        "phenotypes <- pData(gds48558)\n",
        "dim(gds48558)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bgbQ4IMSEEs8"
      },
      "source": [
        "### Normalizing the Data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "v2v7PkQNsHJx",
        "outputId": "8b78b342-7d11-408f-d5b6-dba435fcb299"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[1] 1.611473\n",
            "[1] 13.76154\n"
          ]
        }
      ],
      "source": [
        "%%R\n",
        "print(min(expr))\n",
        "print(max(expr))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Mn_3XzaDmn7N",
        "outputId": "4d430955-0933-4d92-bf39-f49d13845879"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[1] 32321    67\n"
          ]
        }
      ],
      "source": [
        "%%R\n",
        "\n",
        "# log2 transform\n",
        "qx <- as.numeric(quantile(expr, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))\n",
        "LogC <- (qx[5] > 100) || (qx[6]-qx[1] > 50 && qx[2] > 0)\n",
        "if (LogC) { \n",
        "    expr[which(expr <= 0)] <- NaN\n",
        "    expr <- log2(expr) \n",
        "}\n",
        "\n",
        "dim(expr)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qs9WDBOpe4VY"
      },
      "source": [
        "## Differential Expression Analysis"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "7Tpn8uzFbAsO",
        "outputId": "d97c5f30-6106-4c23-f049-393604230786"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "        Gene.symbol Gene.ID    adj.P.Val     logFC\n",
            "8016932         MPO    4353 3.617813e-19 -5.563501\n",
            "7970737        FLT3    2322 4.835716e-19 -5.250065\n",
            "7989647    KIAA0101    9768 6.308160e-19 -4.559135\n",
            "7982663       BUB1B     701 1.664043e-18 -2.756554\n",
            "8083422      SUCNR1   56670 1.938573e-18 -2.996816\n",
            "7926259       MCM10   55388 3.712137e-18 -2.318848\n"
          ]
        }
      ],
      "source": [
        "%%R\n",
        "\n",
        "design <- model.matrix(~group + 0, gds48558)\n",
        "colnames(design) <- levels(gs)\n",
        "\n",
        "fit <- lmFit(gds48558, design)  # fit linear model\n",
        "\n",
        "# set up contrasts of interest and recalculate model coefficients\n",
        "cts <- paste(groups[1], groups[2], sep=\"-\")\n",
        "cont.matrix <- makeContrasts(contrasts=cts, levels=design)\n",
        "fit2 <- contrasts.fit(fit, cont.matrix)\n",
        "\n",
        "# compute statistics and table of top significant genes\n",
        "fit2 <- eBayes(fit2, 0.01)\n",
        "tT <- topTable(fit2, adjust=\"fdr\", sort.by=\"B\", number=Inf)\n",
        "\n",
        "tT <- subset(tT, select=c(\"Gene.symbol\", \"Gene.ID\",\"adj.P.Val\",\"logFC\"))\n",
        "write.table(tT, file='AML_Healthy.tsv', row.names=F, sep=\"\\t\", quote=F)\n",
        "\n",
        "head(tT)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-hnXI43CBnvI",
        "outputId": "4d972c12-fb9b-4a40-d60a-295abdc3dbeb"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[1] \"STK38\" \"CBX7\"  \"PLCL2\" \"PECR\"  \"\"      \"HLA-F\"\n"
          ]
        }
      ],
      "source": [
        "%%R\n",
        "\n",
        "result <- read.table(file = 'AML_Healthy.tsv', sep = '\\t', header = TRUE)\n",
        "aml.up <- subset(result, logFC > 1 & adj.P.Val < 0.05)\n",
        "# aml.up.gene <- unique(aml.up$Gene.symbol)\n",
        "# aml.up.gene <- sub(\"///.*\", \"\", aml.up.gene)\n",
        "aml.up.gene <- unique(as.character(strsplit2(aml.up$Gene.symbol, \"///\")))\n",
        "write.table(aml.up.gene, file='AML_Healthy_Up.tsv', row.names=F, col.names=F, quote=F)\n",
        "\n",
        "aml.down <- subset(result, logFC < -1 & adj.P.Val < 0.05)\n",
        "aml.down.gene <- unique(as.character(strsplit2(aml.down$Gene.symbol, \"///\")))\n",
        "write.table(aml.down.gene, file='AML_Healthy_Down.tsv', row.names=F, col.names=F, quote=F)\n",
        "\n",
        "head(aml.up.gene)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "bgbQ4IMSEEs8"
      ],
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}