Tutorial_3_DEG_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RNA-Seq Analysis Training Demo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial will show you the following workflow that is used to generate the list of differentially expressed genes. Running the code in this tutorial should take approximately 15 minutes.\n",
    "\n",
    "This overflow centers roughly around the useage of the deseq2 package. Documentation for this package can be found here: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html\n",
    "\n",
    "![RNA-Seq workflow](images/count-workflow.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 1: Library Installation\n",
    "\n",
    "First install the R packages that will be used. <strong>This installation may take around 15 minutes.</strong>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###install the libraries first, if required\n",
    "install.packages(\"BiocManager\")\n",
    "\n",
    "if (!require(\"ComplexHeatmap\"))\n",
    "    BiocManager::install(\"ComplexHeatmap\")\n",
    "\n",
    "if (!require(\"DESeq2\"))\n",
    "    BiocManager::install(\"DESeq2\")\n",
    "\n",
    "if (!require(\"EnhancedVolcano\"))\n",
    "    BiocManager::install('EnhancedVolcano')\n",
    "\n",
    "if (!require(\"dplyr\"))\n",
    "    install.packages(\"dplyr\")\n",
    "\n",
    "if (!require(\"pheatmap\"))\n",
    "    install.packages(\"pheatmap\")\n",
    "\n",
    "if (!require(\"ggrepel\"))\n",
    "    install.packages(\"ggrepel\")\n",
    "\n",
    "if (!require(\"ggfortify\"))\n",
    "    install.packages(\"ggfortify\")\n",
    "\n",
    "\n",
    "###load the libraries\n",
    "library(DESeq2)\n",
    "library(dplyr)\n",
    "library(ComplexHeatmap)\n",
    "library(EnhancedVolcano)\n",
    "library(ggplot2)\n",
    "library(ggrepel)\n",
    "library(ggfortify)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 2: Reading in Data\n",
    "\n",
    "Next read-in the gene count data that will be analyzed.\n",
    "\n",
    "This read-in step often involves reformatting or adjusting the data.\n",
    "\n",
    "To make things easier, we created a matrix, or a table, of just our read-count data -- excluding the non-numerical metadata columns. This makes it easy for us to feed this table into the deseq2 tool later on.\n",
    "\n",
    "Additionally, later on we will use normalization involving logarithms. So, to prevent errors, we will replace all readcount values of '0' in our data, with '1'. This will change the data only slightly, and will prevent these zeroes from causing an 'undefined' or 'N/A' math error in the logarithm normalization step.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "system(\"gsutil cp gs://nigms-sandbox/me-inbre-rnaseq-pipelinev2/data/salmon/SRP300216_merged_quants.txt --output data/gene_counts/SRP300216_genecounts.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#download genecounts and annotation file from cloud storage bucket\n",
    "system(\"mkdir -p data/gene_counts\", intern=TRUE)\n",
    "system(\"mkdir -p data/reference\", intern=TRUE)\n",
    "system(\"gsutil cp gs://nigms-sandbox/me-inbre-rnaseq-pipelinev2/data/salmon/SRP300216_merged_quants.txt data/gene_counts/SRP300216_genecounts.txt\")\n",
    "#try downloading with curl\n",
    "system(\"curl https://storage.googleapis.com/nigms-sandbox/me-inbre-rnaseq-pipelinev2/data/reference/GCF_001632805.1_ASM163280v1_feature_table.txt --output data/reference/GCF_001632805.1_ASM163280v1_feature_table.txt\", intern=TRUE)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read-in the raw gene count file to a dataframe variable we named 'read_counts'\n",
    "read_counts <- read.table('data/gene_counts/SRP300216_genecounts.txt',head=TRUE)\n",
    "\n",
    "#replace all count entries of '0' in the raw gene count file with '1'.\n",
    "read_counts[read_counts==0] <- 1\n",
    "\n",
    "#assign the numerical read counts to a matrix variable we named 'rnaseqMatrix'\n",
    "rnaseqMatrix <- round(read_counts[,c(2:13)])\n",
    "\n",
    "#label the rownames of this matrix with the rowname column from the gene count dataframe.\n",
    "rownames(rnaseqMatrix) <- read_counts[,1]\n",
    "\n",
    "head(rnaseqMatrix)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 4: Specifying Experimental Design\n",
    "\n",
    "Next specify the experimental design. \n",
    "\n",
    "The deseq2 tool will later use this design to group samples together, and output information about the statistical differences in gene expression between these specified groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the sample experimental design, in this case 6 wildtype and 6 bacteriophage infected samples\n",
    "samples_treatment <- data.frame(matrix(c(rep(\"WT\",6),rep(\"BPs_lysogen\",6)),ncol=1))\n",
    "samples_ID <- data.frame(matrix(colnames(rnaseqMatrix),ncol=1))\n",
    "samples <- cbind(samples_ID,samples_treatment)\n",
    "names(samples) <- c(\"ID\",\"Treatment\")\n",
    "rownames(samples) <- samples[,1]\n",
    "\n",
    "print(\"An example of how a deseq2 experimental design table might look.\")\n",
    "samples\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 5: Creating Deseq2 Object\n",
    "\n",
    "Now use the treatment design matrix in combination with the readcount matrix to create a deseq2 object. \n",
    "\n",
    "Once created, this is also a good opportunity to filter out lowly expresseed genes, and to inspect the pre-normalized data using pairwise comparison plots.\n",
    "\n",
    "Finally, the deseq2 analysis can be run on the deseq2 object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create the DEseq2DataSet object\n",
    "deseq2Data <- DESeqDataSetFromMatrix(countData = rnaseqMatrix, colData = samples, design = ~ Treatment)\n",
    "\n",
    "# Perform pre-filtering of the data\n",
    "deseq2Data <- deseq2Data[rowSums(counts(deseq2Data)) > 10, ]\n",
    "\n",
    "# Inspect pre-normalized data \n",
    "pairs(log2((counts(deseq2Data)+1)))\n",
    "boxplot(log2((counts(deseq2Data)+1)))\n",
    "\n",
    "# Run pipeline for differential expression steps\n",
    "deseq2Data <- DESeq(deseq2Data)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 6: Normalize Data\n",
    "\n",
    "Now normalize the data using the 'regular normalization' function from deseq2.\n",
    "\n",
    "Normalization of gene readcounts across samples will help to minimize sample variation noise and ensure more accurate comparisons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#perform rlog normalization on deseq2 object\n",
    "rld <- rlog(deseq2Data, blind=FALSE)\n",
    "\n",
    "#a handy variable so we can easily reference the rlog genecount matrix in the future, as deseq objects contain many tables of data.\n",
    "rlogcounts <- assay(rld)\n",
    "rownames(rlogcounts) <- rownames(deseq2Data)\n",
    "\n",
    "head(rlogcounts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalization of data can is often verified by inspecting by various plots, for instance histograms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist(rlogcounts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Boxplots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boxplot(rlogcounts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or pairwise comparison plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs(rlogcounts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the plots of the pre-normalized and normalized data to see the effect of normalization.\n",
    "\n",
    "Notice how the samples in the box plot, and in-group samples in the pair plots display similar distributions.\n",
    "\n",
    "Looking at normalized plots can be an easy way to look out for any obvious inconsistencies or errors in the sample data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 7: PCA Plot\n",
    "\n",
    "We can now begin to produce plots to further analyze the differences between groups.\n",
    "\n",
    "PCA plots can be an effective way to visualize variation within and between sample sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_counts <- prcomp(t(rlogcounts))\n",
    "autoplot(pca_counts, data = samples, colour=\"Treatment\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plots can be output as above, or also saved as pdf or image files, as below. Try to navigate and click to the created file, after you run the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf(\"data/gene_counts/PCA_Plot.pdf\")\n",
    "autoplot(pca_counts, data = samples, colour=\"Treatment\")\n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 8: Contrast Comparisons\n",
    "\n",
    "Contrast comparisons between two treatment groups can be performed using the deseq2 tool to identify differentially up and down regulated genes.\n",
    "\n",
    "Using R, these results can be summarized, reordered, and/or trimmed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Output the results of comparing WT and BPs_lysogen group\n",
    "res_WT_BPs <- results(deseq2Data, contrast=c(\"Treatment\",\"WT\",\"BPs_lysogen\"))\n",
    "summary(res_WT_BPs)\n",
    "print(\"Number of genes under a p value of 0.05\")\n",
    "sum(res_WT_BPs$pvalue < 0.05, na.rm=TRUE)\n",
    "print(\"Number of genes under a p adjusted value of 0.05\")\n",
    "sum(res_WT_BPs$padj < 0.05, na.rm=TRUE)\n",
    "\n",
    "\n",
    "#deseq2 optionally output some padj values as 'NA' when genecounts meet certain criteria\n",
    "#that criteria can be found here: https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pvaluesNA\n",
    "#for practical uses, one could replace the NA values with '1', essentially NA values not significant.\n",
    "res_WT_BPs$padj[is.na(res_WT_BPs$padj)] <- 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 9: Annotating and Exporting Results\n",
    "\n",
    "These differentially expressed genes can be annotated and exported.\n",
    "\n",
    "Combining gene counts with gene annotation isn't always simple, and will vary heavily depending on the annotation file itself.\n",
    "\n",
    "For example, when merging genecounts by 'gene id' with feature tables from NCBI's assembly database, its very common to have duplicate rows, as a single 'gene id' may match to multiple features. \n",
    "\n",
    "Depending on how you want to use your annotated results, you may decide to further reformat your table in whatever suits your further analysis best, for instance by downloading and using a spreadsheet editor.\n",
    "\n",
    "Because the transcript table for m chelonae we have here is so simple, and contains essentially only one 'CDS' or '..RNA' feature for every gene identifier we are using, it is easy for us to reformat it here in the way we want just using R code. However, it is important to note that, again, this is very often not the case, and reformatting will vary depending on both the annotation file itself, and what that annotation file will be used for."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read-in the annotation file (see extended tutorial for more info on how to get annotation files) and store it into a variable.\n",
    "annotation <- read.table(\"data/reference/GCF_001632805.1_ASM163280v1_feature_table.txt\", header=TRUE,sep='\\t', quote=\"\\\"\",comment.char = \"\")\n",
    "\n",
    "#put the previous deseq results into an easier to manipulate data.frame format\n",
    "results <- data.frame(res_WT_BPs)\n",
    "\n",
    "#merge the annotation with the deseq2 contrast comparison table\n",
    "results.annot <- merge(results,annotation,by.x=0,by.y=17,all.x=TRUE)\n",
    "\n",
    "\n",
    "\n",
    "#Open the annotation file above to see the format. If you further download and examine in a spreadsheet editor,\n",
    "#you can see there is one CDS or ..RNA identifier for each 'gene' feature\n",
    "#this makes reformatting in R simple by just removing every 'gene' feature, and leaving the more detailed 'cds' or '...RNA' feature.\n",
    "#this is very often not always the case, but simple and easy here.\n",
    "results.annot <- results.annot[annotation$X..feature!='gene',]\n",
    "\n",
    "#labeling the columns properly after the merge\n",
    "colnames(results.annot)[1] <- 'locus_tag'\n",
    "\n",
    "#write out the file\n",
    "write.table(results.annot,\"data/gene_counts/WT_BPs_DESeq2_annot.txt\",row.names=FALSE,sep=\"\\t\")\n",
    "\n",
    "head(results.annot)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STEP 10: More Plots and Clustering Data\n",
    "\n",
    "Finally, although several other different kind of analysis can be done, it is common to plot the results of differential gene expression analysis in forms such as volcano, MA, and heatmaps. These plots help one to further examine the output list of differentially expressed genes in different ways, and potentially identify or further investigate candidate gene(s)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#volcano and ma plots can be used to visualize significant differences in gene expression\n",
    "results_volc_df <- results.annot\n",
    "EnhancedVolcano(results_volc_df,\n",
    "    lab = results_volc_df[,1],\n",
    "    x = 'log2FoldChange',\n",
    "    y = 'padj')\n",
    "\n",
    "#ma plot\n",
    "plotMA(res_WT_BPs)\n",
    "\n",
    "\n",
    "#heatmaps can used to look at clustering and expression of various gene lists.\n",
    "filtered_res <- results.annot %>% filter(padj < 0.05)\n",
    "filtered_res_counts <- subset(rlogcounts, rownames(rlogcounts) %in% filtered_res[,1])\n",
    "filtered_res_counts = t(scale(t(filtered_res_counts)))\n",
    "Heatmap(filtered_res_counts[1:50,], name = \"Z-Score\", row_names_gp = gpar(fontsize = 8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <a name=\"workflow\">Additional Workflows</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that you have ran some analysis on read count data, feel feel free to explore and revisit the previous workflows for creating readcounts, such as the standard or extended tutorials, or the snakemake tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "[Workflow One:](Tutorial_1.ipynb) A short introduction to downloading and mapping sequences to a transcriptome using Trimmomatic and Salmon. Here is a link to the YouTube video demonstrating the tutorial: <https://youtu.be/ChGfBR4do_Y>.\n",
    "\n",
    "[Workflow One (Extended):](Tutorial_1B_Extended.ipynb) An extended version of workflow one. Once you have got your feet wet, you can retry workflow one with this extended version that covers the entire dataset, and includes elaboration such as using SRA tools for sequence downloading, and examples of running batches of fastq files through the pipeline. This workflow may take around an hour to run.\n",
    "\n",
    "[Workflow One (Using Snakemake):](Tutorial_2_Snakemake.ipynb) Using snakemake to run workflow one.\n",
    "\n",
    "[Workflow Two (DEG Analysis):](Tutorial_3_DEG_Analysis.ipynb) Using Deseq2 and R to conduct clustering and differential gene expression analysis.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![RNA-Seq workflow](images/RNA-Seq_Notebook_Homepage.png)"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 4
}