DESeq2_DEG_from_counts.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyP938f74jqgGDEYzLbh+4El",
      "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/DESeq2_DEG_from_counts.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(\"vsn\")"
      ],
      "metadata": {
        "id": "UvIbt0oTzyLP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "install.packages(\"pheatmap\")"
      ],
      "metadata": {
        "id": "RIXIgM3d0BC-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "if (!require(\"BiocManager\", quietly = TRUE))\n",
        "    install.packages(\"BiocManager\")\n",
        "\n",
        "BiocManager::install(\"DESeq2\")"
      ],
      "metadata": {
        "id": "COzTQMSq0CTv"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#--------------------------htseq-count_input---------------------------"
      ],
      "metadata": {
        "id": "W5g3oOhuwTie"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "dir <- 'directory_path'"
      ],
      "metadata": {
        "id": "kdfxADMjxLen"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "sampleFiles <- grep(\".count\", list.files(dir),value=TRUE, fixed=TRUE)\n",
        "coldata <- read.csv('descriptor.csv', header=TRUE, sep='\\t')"
      ],
      "metadata": {
        "id": "3s2EewyDwXtG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "rownames(coldata) <- coldata$ID"
      ],
      "metadata": {
        "id": "gfidd-D1wm1W"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "coldata$PT_ID <- as.factor(coldata$PT_ID)\n",
        "levels(coldata$PT_ID)=c('PT_1', 'PT_2', 'PT_3', 'PT_4', 'PT_5')"
      ],
      "metadata": {
        "id": "vJPGuErXw09z"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#add condition column($)\n",
        "coldata$condition <- factor(coldata$MOL_MARKER)"
      ],
      "metadata": {
        "id": "dwGVR0aRw9CX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "library(\"DESeq2\")"
      ],
      "metadata": {
        "id": "jj3YnRBwzs0W"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "dds <- DESeqDataSetFromHTSeqCount(sampleTable = coldata,\n",
        "\n",
        "                                  directory = dir,\n",
        "\n",
        "                                  design = ~ condition)"
      ],
      "metadata": {
        "id": "A8wFU1IaxDDH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "dds"
      ],
      "metadata": {
        "id": "H-uitHKwxTPN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#-----------------------------pre-filtering-----------------------------"
      ],
      "metadata": {
        "id": "RiPilhDXxWht"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "keep <- rowSums(counts(dds)) >= 10          #the choosen value depends on experimental design\n",
        "dds <- dds[keep,]\n",
        "dds"
      ],
      "metadata": {
        "id": "50o2c1ufxXcW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#-------------------------------DEG-------------------------------------"
      ],
      "metadata": {
        "id": "tAzkjr7Axq_O"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "dds <- DESeq(dds)\n",
        "res <- results(dds)\n",
        "res"
      ],
      "metadata": {
        "id": "OYMbYnVQxroG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "write.csv(res,'DEG.csv', row.names=TRUE)"
      ],
      "metadata": {
        "id": "LahZPVZRxvbI"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#------------------------------Analysis------------------------------------"
      ],
      "metadata": {
        "id": "aznUW554x5aP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "resOrdered <- res[order(res$pvalue),]\n",
        "summary(res)"
      ],
      "metadata": {
        "id": "nFRFMLq1x5_e"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "sum(res$padj < 0.1, na.rm=TRUE)"
      ],
      "metadata": {
        "id": "ZLm0-LhSx9xo"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "res05 <- results(dds, alpha=0.05)\n",
        "summary(res05)"
      ],
      "metadata": {
        "id": "wPQ9oHHxyA_n"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "sum(res05$padj < 0.05, na.rm=TRUE)"
      ],
      "metadata": {
        "id": "3a-IJhX_yDje"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#--------------------------------MA plots---------------------------------"
      ],
      "metadata": {
        "id": "fgOoQkBOyHVH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plotMA(res, ylim=c(-10,10))"
      ],
      "metadata": {
        "id": "uvZ65JbpyJp2"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plotMA(resLFC, ylim=c(-6,6))"
      ],
      "metadata": {
        "id": "AzZAtSD_yMD_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "idx <- identify(res$baseMean, res$log2FoldChange)\n",
        "rownames(res)[idx]"
      ],
      "metadata": {
        "id": "EtOiHXVJyP3_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#---------------------------plot counts----------------------------"
      ],
      "metadata": {
        "id": "tOApXEmoyS9-"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "library(\"ggplot2\")"
      ],
      "metadata": {
        "id": "A5ETSruyznUX"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plotCounts(dds, gene=which.min(res$padj), intgroup=\"condition\")\n",
        "\n",
        "d <- plotCounts(dds, gene=which.min(res$padj), intgroup=\"condition\", \n",
        "                returnData=TRUE)\n",
        "                \n",
        "ggplot(d, aes(x=condition, y=count)) + \n",
        "  geom_point(position=position_jitter(w=0.1,h=0)) + \n",
        "  scale_y_log10(breaks=c(25,100,400))"
      ],
      "metadata": {
        "id": "HVM1anaiyUZu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#----------------------------------report--------------------------------"
      ],
      "metadata": {
        "id": "A25a1vSyyY-m"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "mcols(res)$description"
      ],
      "metadata": {
        "id": "GVxx1Ibnyaff"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "write.csv(as.data.frame(resOrdered), \n",
        "          \"report_ordered.csv\")"
      ],
      "metadata": {
        "id": "lWI4J-j6ydZX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "resSig <- subset(resOrdered, padj <= 0.05)"
      ],
      "metadata": {
        "id": "temV7VwkyhUG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "write.csv(resSig, \"report_p005.csv\")"
      ],
      "metadata": {
        "id": "w1jOTNp8yknH"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#-------------------------Data transformation-------------------------"
      ],
      "metadata": {
        "id": "praNzoHwyt9A"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "vsd <- vst(dds, blind=FALSE)\n",
        "rld <- rlog(dds, blind=FALSE)\n",
        "head(assay(vsd), 3)"
      ],
      "metadata": {
        "id": "JBkCI9GvyvyW"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "write.csv(assay(vsd), \"matrix_vsdlog.csv\" )"
      ],
      "metadata": {
        "id": "k4NLvcVCyxdm"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "library(\"vsn\")"
      ],
      "metadata": {
        "id": "6kfC8EJVzi2H"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# this gives log2(n + 1)\n",
        "ntd <- normTransform(dds)\n",
        "\n",
        "meanSdPlot(assay(ntd))\n",
        "\n",
        "meanSdPlot(assay(vsd))\n",
        "\n",
        "meanSdPlot(assay(rld))"
      ],
      "metadata": {
        "id": "8KHcSpV8y3pY"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#---------quality assessment, clustering and visualization-------------"
      ],
      "metadata": {
        "id": "K_aHZ_2qy8SV"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "library(\"pheatmap\")"
      ],
      "metadata": {
        "id": "mLAr2dgGy81u"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#vsd\n",
        "select <- order(rowMeans(counts(dds,normalized=TRUE)),\n",
        "                decreasing=TRUE)[1:20]\n",
        "df <- as.data.frame(colData(dds)[,c(\"condition\", \"SAMPLE_ORIGIN\")])\n",
        "pheatmap(assay(vsd), cluster_rows=FALSE, show_rownames=FALSE,\n",
        "         cluster_cols=FALSE, annotation_col=df)"
      ],
      "metadata": {
        "id": "IBzSMa3Oy-zN"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#clustering\n",
        "pheatmap(assay(vsd), cluster_rows=F, show_rownames=F,\n",
        "         cluster_cols=T, annotation_col=df)"
      ],
      "metadata": {
        "id": "qUIoC5lHzJPg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#sampletosample distances, distance matrix\n",
        "sampleDists <- dist(t(assay(vsd)))"
      ],
      "metadata": {
        "id": "vbI15idUzLE8"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#PCA\n",
        "plotPCA(vsd, intgroup=c(\"condition\"))"
      ],
      "metadata": {
        "id": "67QaIlShzRVf"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "plotPCA(vsd, intgroup=c(\"condition\", \"SAMPLE_ORIGIN\"))"
      ],
      "metadata": {
        "id": "cJngNByVzSsg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "pcaData <- plotPCA(vsd, intgroup=c(\"condition\", \"SAMPLE_ORIGIN\"), returnData=TRUE)\n",
        "percentVar <- round(100 * attr(pcaData, \"percentVar\"))\n",
        "ggplot(pcaData, aes(PC1, PC2, color=condition, shape=SAMPLE_ORIGIN)) +\n",
        "    geom_point(size=2) +\n",
        "    xlab(paste0(\"PC1: \",percentVar[1],\"% variance\")) +\n",
        "    ylab(paste0(\"PC2: \",percentVar[2],\"% variance\")) + \n",
        "    coord_fixed()"
      ],
      "metadata": {
        "id": "6upmmwT8zW0Y"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "pcaData <- plotPCA(vsd, intgroup=c(\"condition\", \"SAMPLE_ORIGIN\"), returnData=TRUE)\n",
        "percentVar <- round(100 * attr(pcaData, \"percentVar\"))\n",
        "ggplot(pcaData, aes(PC1, PC2, color=condition, shape=PT_ID)) +\n",
        "  geom_point(size=2) +\n",
        "  xlab(paste0(\"PC1: \",percentVar[1],\"% variance\")) +\n",
        "  ylab(paste0(\"PC2: \",percentVar[2],\"% variance\")) + \n",
        "  coord_fixed()"
      ],
      "metadata": {
        "id": "wKZ5IDU3zXMY"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#dispersionestimates\n",
        "plotDispEsts(dds)"
      ],
      "metadata": {
        "id": "TO9joBc8zbme"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}