bulkRNASeq_tutorial_jupyter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2371c328",
   "metadata": {},
   "source": [
    "# RNA-Seq Hands-on\n",
    "\n",
    "## Introduction\n",
    "\n",
    "### Intended learning outcomes\n",
    "\n",
    "This tutorial is meant to show you the basics of a bulk tissue RNA-Seq differential expression analysis. Today, you will learn...\n",
    "\n",
    "* terminology of bulk RNA-Seq bioinformatics analyses.\n",
    "* to name some important bulk RNA-Seq tools and describe their basic functionalities.\n",
    "* to name and explain the steps of a bulk RNA-Seq analysis workflow.\n",
    "* to explain the most important plot types and how to plot and visualize bulk RNA-Seq results with ggplot2.\n",
    "* to apply R code to your bulk RNA-Seq data to identify differentially expressed genes, plot results, and find enriched gene ontology terms.\n",
    "* to interpret the results you get from bulk RNA-Seq analyses.\n",
    "\n",
    "The nuances of bulk RNA-Seq data generation like experimental design and choice of library preparation methods are not part of this lecture and tutorial. If you'd like to look into that, we designed an OLAT-course which takes approx. 60 minutes on that topic. You can find it at https://lms.uni-kiel.de/url/Catalog/0/Search/0/Infos/5273846135 using biomedinf2022 as a password.\n",
    "\n",
    "### Roadmap\n",
    "\n",
    "This is the roadmap for today.\n",
    "\n",
    "1. Initialize R environment\n",
    "2. Inspect FastQC results\n",
    "3. Read files, data preprocessing\n",
    "4. Exploratory data analysis with PCA\n",
    "5. \n",
    "    a. Differential expression analysis, using DESeq2\n",
    "    b. Volcano plot\n",
    "    c. Heatmap\n",
    "6. Enrichment analysis of differentially expressed genes\n",
    "\n",
    "Most of the sections contain tasks. The tasks under the \"Tasks\" headings are basic tasks to test yourselves if you understood what just happened in the code. The \"Expert Questions\" tasks are a little more complicated and require some fiddling with code. We recommend to go through the script and the basic tasks first, and then continuing with advanced tasks you find interesting.\n",
    "\n",
    "Try to solve the tasks without directly looking up the solution to improve your learning. Try it below!\n",
    "\n",
    "#### Tasks"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5efec50",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "Example task: Which tissue was sampled to generate the data in this tutorial?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7cc8aea",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden"
   },
   "source": [
    "The data was generated from full peripheral blood samples."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9405e16",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "### Terminology\n",
    "\n",
    "Let's briefly review the terminology:\n",
    "\n",
    "* *bulk tissue* as opposed to *single-cell* means that a tissue sample consisting of many cells and possibly different cell types was analyzed. The tissue used here is full peripheral blood.\n",
    "* *differential expression* means that the relative abundance of many types of transcripts is compared between samples. Here we compare samples of healthy donors and patients infected with COVID-19 in different stages of the disease.\n",
    "\n",
    "### Workflow\n",
    "\n",
    "The bulk RNA-Seq workflow can be thought as a procedure comprising these steps:\n",
    "\n",
    "1. Library prep and Sequencing\n",
    "2. Quality control of sequencing reads\n",
    "3. Gapped alignment of sequencing reads to a reference genome\n",
    "4. Feature counting (counting reads with genes as bins)\n",
    "5. Differential expression analysis\n",
    "6. Gene ontology enrichment analysis\n",
    "\n",
    "This hands-on tutorial focuses on steps 5. and 6. We did the steps 1.-4. for you using the nf-core/rnaseq pipeline. These steps can be standardized and pipelining them saves us a lot of time."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9874c828",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "## Initialize R environment\n",
    "**The following applies when working on this course elsewhere and can be ignored when working on the JupyterHub:**\n",
    "\n",
    "Unzip the directory you received from the course assistants and change the R working directory into the unzipped folder named 'burnhot' using `setwd(\"path/to/burnhot\")`.\n",
    "\n",
    "Alternatively, download the repository from github in a shell `git clone https://github.com/ikmb/burnhot.git` or from the website of the repo <https://github.com/ikmb/burnhot>.\n",
    "\n",
    "Change the working directory in the R console into the burnhot directory `setwd(\"path/to/burnhot\")`.\n",
    "\n",
    "Before installing R libraries, you need to install cairolib on your system using bash ```apt-get install libcairo2-dev```.\n",
    "\n",
    "**Today you can execute the cells within this Jupyter notebook directly. Select for this the 'rnaseq_microbiome'-kernel. Where all the libraries needed are preinstalled.**\n",
    "\n",
    "Next, load necessary external libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45f6642e",
   "metadata": {
    "deletable": false,
    "editable": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(data.table)\n",
    "library(tidyverse)\n",
    "library(ggrepel)\n",
    "library(DESeq2)\n",
    "library(ComplexHeatmap)\n",
    "library(RColorBrewer)\n",
    "library(topGO)\n",
    "library(Rgraphviz)\n",
    "library(genefilter)\n",
    "library(geneplotter)\n",
    "library(org.Hs.eg.db)\n",
    "folder <- \"~/bioinformatics_materials/bulkRNASeq/results\"\n",
    "\n",
    "if (!dir.exists(folder)){\n",
    "  dir.create(folder)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e18a0ca",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "## Inspect FastQC results\n",
    "\n",
    "After running the pipeline, the tool *FastQC* outputs a report on the quality of the sequencing data. Using the tool *MultiQC*, the FastQC result and log output from other tools were compiled into one html document. \n",
    "\n",
    "#### Tasks\n",
    "\n",
    "Open the report <multiqc_report1.html> in your file-browser and focus on the section \"FastQC\". It contains some quality checks, in which each sample is graded as \"pass\", \"warn\" or \"fail\". What do you find? Think for a minute if you can come up with an explanation. If you can't, no problem, as this will become clear in the following. \n",
    "\n",
    "## Read files and preprocessing\n",
    "\n",
    "Before the tutorial, we have run the RNA-Seq pre-processing pipeline for you. The main result is a count table, with sample IDs (the patients) as columns, gene IDs as rows, and transcript counts as the values. From the clinics, we know who of those patients was diagnosed with COVID-19, and how severe their status was. This data is saved as tables in the project folder you downloaded.\n",
    "\n",
    "We now load the data and bring it into a format that we can run the analysis functions on. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28ce3c25",
   "metadata": {
    "deletable": false,
    "editable": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "##Read sample info file\n",
    "covid_info_file <- read.csv(\"~/bioinformatics_materials/bulkRNASeq/data/Sequenced_sample_metadata_020720 (1).csv\")\n",
    "covid_info_file$Pseudotime <- as.character(covid_info_file$Pseudotime)\n",
    "##Read raw count file\n",
    "count_data <- read.csv(\"~/bioinformatics_materials/bulkRNASeq/data/COVID_control_raw_counts_070720 (1).txt\", sep = '\\t')\n",
    "count_data <- count_data[, as.character(covid_info_file$RNAseqID)]\n",
    "##Format sample col data, the patient data\n",
    "col_data <- covid_info_file\n",
    "rownames(col_data) <- col_data$RNAseqID\n",
    "col_data$Gender <- as.factor(col_data$Gender)\n",
    "col_data$Disease.traj <- as.factor(col_data$Disease.traj)\n",
    "col_data$Pseudotime <- as.factor(as.character(col_data$Pseudotime))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4269ddd1",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "#### Tasks\n",
    "\n",
    "1. Now, use the function ```head()``` to get an impression what your data actually look like. What do you need to type to get a first glance at the data?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c885a1b",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden"
   },
   "source": [
    "We can have a brief look at the sample info table now and make sure that the samples in the count table have the same naming pattern as those in the sample info table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e67fe7f2",
   "metadata": {
    "deletable": false,
    "editable": false,
    "scrolled": true,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "head(col_data)\n",
    "head(count_data)\n",
    "#View(col_data) #Run to display the whole phenotype table. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d6126e0",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Run your code here:\n",
    "head(col_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1d0a887",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "2. Which function allows you a quick look at the end of a table? Run ```?head``` if you don't know the answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6bf1d35",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "tail(col_data)\n",
    "tail(count_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b339197d",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Run your code here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b65168d",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "It is always a good idea to check if the numbers of columns and rows are as you expect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "806e9abb",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "dim(col_data)\n",
    "dim(count_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfe5e4b6",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "3. Are the number of rows in the sample info table same as the number of columns in the count table? What does the number of rows in the count table signify?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bdc8390",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden"
   },
   "source": [
    "It is the number of genes that reads were assigned to during the feature counting step. It is the total number of genes and other features in our annotation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b53397f",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "## Exploratory data analysis with PCA and heatmaps\n",
    "\n",
    "We would like to get an overview of our data, to see if samples from the same stage of the disease are more similar to each other than to healthy controls.\n",
    "\n",
    "Principal Component Analysis (PCA) is a statistical method to reduce the dimensionality of a dataset. Here, we want to project our count data into a 2-dimensional space on our screen. The first principal component describes a vector in the real coordinate space which minimizes the average square distance of the data points to the line. Each subsequent component tries the same, but has to be orthogonal to all other vectors (<https://en.wikipedia.org/wiki/Principal_component_analysis>). Sometimes, principal components capture biological insights to our data, because they cluster samples with similar expression profiles.\n",
    "\n",
    "#### Task\n",
    "\n",
    "Before you run the section below, think about what you expect to see in the PCA plot for this dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7bc6983",
   "metadata": {
    "deletable": false,
    "editable": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "##Calculate the principle components for the count table\n",
    "pc <-prcomp(t(count_data))\n",
    "pc_df <- as.data.frame(pc$x)\n",
    "pc_df$Disease.traj <- factor(col_data$Disease.traj, \n",
    "                             levels = c(\"Healthy\",\"Complicated (incremental)\",\n",
    "                                        \"Critical\",\"Complicated with hyperinflammatory syndrom\",\n",
    "                                        \"Complicated (recovering)\",\"Mild (recovering)\",\n",
    "                                        \"recovery/asymptomatic\",\"Recovered\"), \n",
    "                             ordered = T)\n",
    "pc_df$Patient <- col_data$PatientID\n",
    "\n",
    "## check the fraction of the total variance that is captured by the principal components\n",
    "summary(pc)\n",
    "\n",
    "pseudotime_colors <- c(\"#A7A9AC\",\"#E84F8C\", \"#7031B9\", \"#A04E9E\", \"#E65826\", \"#F99B1C\", \"#FDC077\", \"#51BBFE\")\n",
    "P <- ggplot(data = pc_df, mapping = aes(x=PC1, y=PC2, color=Disease.traj, label=Patient)) + \n",
    "  geom_point(size=5) + geom_text_repel(hjust=0, nudge_x = 2, nudge_y = 1, cex=3)\n",
    "P <- P + scale_color_manual(values = pseudotime_colors)\n",
    "P <- P + theme_bw() + theme(axis.text = element_text(size = 20), axis.title = element_text(size = 22))\n",
    "##Set output size of the plot. \n",
    "options(repr.plot.width=10, repr.plot.height=8)\n",
    "P\n",
    "ggsave(P,file=paste0(folder,\"/PCA.png\"),height = 8,width = 10, units = \"in\", dpi = 300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3d6627a",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "#### Tasks\n",
    "\n",
    "Compare the PCA plot above with your expectation. Are the samples clustered by the disease status of the patient?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0366cce1",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden"
   },
   "source": [
    "The clustering pattern of COVID-cases and controls is not as obvious as we might have expected. The controls are somewhat concentrated, but not convincingly. This might be because we have not transformed the count data yet. The DESeq2 library can help us out."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9dc44513",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "### Data transformation\n",
    "\n",
    "The data transformation by DESeq2 consists of correcting for the variable number of sequenced molecules, for transcripts with too high variance that we consider a statistical artifact and putting the values on a log-scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6f83a70",
   "metadata": {
    "deletable": false,
    "editable": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "dds_counts <- DESeqDataSetFromMatrix(countData = count_data, colData = col_data, design = ~ Pseudotime)\n",
    "dds_counts <- dds_counts[ rowSums(counts(dds_counts)) > 1, ]\n",
    "dds_counts <- estimateSizeFactors(dds_counts)\n",
    "\n",
    "vst_counts <- vst(dds_counts, blind = TRUE)\n",
    "vst.norm.counts <- assay(vst_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84ebb95b",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "Let's try again with the PCA, using transformed data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e7445d7",
   "metadata": {
    "deletable": false,
    "editable": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "pc <-prcomp(t(vst.norm.counts))\n",
    "pc_df <- as.data.frame(pc$x)\n",
    "pc_df$Disease.traj <- factor(col_data$Disease.traj, levels = c(\"Healthy\",\"Complicated (incremental)\",\"Critical\",\"Complicated with hyperinflammatory syndrom\",\n",
    "                                                               \"Complicated (recovering)\",\"Mild (recovering)\",\"recovery/asymptomatic\",\"Recovered\"), ordered = T)\n",
    "pc_df$Patient <- col_data$PatientID\n",
    "pc_df$Gender <- col_data$Gender\n",
    "pc_df$Age <- col_data$Age\n",
    "\n",
    "summary(pc)\n",
    "pseudotime_colors <- c(\"#A7A9AC\",\"#E84F8C\", \"#7031B9\", \"#A04E9E\", \"#E65826\", \"#F99B1C\", \"#FDC077\", \"#51BBFE\")\n",
    "P <- ggplot(data = pc_df, mapping = aes(x=PC1, y=PC2, color=Disease.traj, label=Patient)) + \n",
    "  geom_point(size=5) + geom_text(hjust=0, nudge_x = 2, nudge_y = 1, cex=3) + xlab(\"PC1\") + \n",
    "  ylab(\"PC2\")\n",
    "P <- P + scale_color_manual(values = pseudotime_colors)\n",
    "P <- P + theme_bw() + theme(axis.text = element_text(size = 20), axis.title = element_text(size = 22))\n",
    "P\n",
    "ggsave(P,file=paste0(folder,\"/PCA_vst.png\"),height = 10,width = 10, units = \"in\", dpi = 300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af0e33f7",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "#### Tasks\n",
    "\n",
    "1. Have a look at how the result changed. Does it support our hypothesis about differences between patients and controls?\n",
    "\n",
    "2. Can you change the ggplot command above so that you can see the sex or the age of the patients from the color of the points? Hint: The ```scale_color_manual()``` might give funny results when used with Age as variable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2f0d3e1",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "P <- ggplot(data = pc_df, mapping = aes(x = PC1, y = PC2, color = Gender, label = Patient)) + \n",
    "  geom_point(size = 5) + \n",
    "  geom_text(hjust = 0, nudge_x = 2, nudge_y = 1, cex = 3) + \n",
    "  xlab(\"PC1\") + \n",
    "  ylab(\"PC2\")\n",
    "P\n",
    "P <- ggplot(data = pc_df, mapping = aes(x = PC1, y = PC2, color = Age, label = Patient)) + \n",
    "  geom_point(size = 5) + \n",
    "  geom_text(hjust = 0, nudge_x = 2, nudge_y = 1, cex = 3) + \n",
    "  xlab(\"PC1\") + \n",
    "  ylab(\"PC2\")\n",
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a1873b4",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Run your solution here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfdbc483",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "## Differential expression analysis, using DESeq2\n",
    "\n",
    "We have seen in the PCA that patients in a critical state of their COVID-19 disease cluster away from the controls. Now we want to know which transcripts are differentially expressed between healthy controls and patients with COVID-19 with different disease severity (pseudotime)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97c1536e",
   "metadata": {
    "deletable": false,
    "editable": false,
    "scrolled": true,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#Read the file that maps gene Ensembl ID to gene name so that we can understand the genes better later on\n",
    "ensg2gene <- read.table(\"~/bioinformatics_materials/bulkRNASeq/data/ensg2gene.txt\", header = T, sep=\",\")\n",
    "unique_ensg2gene <- subset(ensg2gene, duplicated(ensg2gene$ensembl_gene_id) == FALSE)\n",
    "rownames(unique_ensg2gene) <- unique_ensg2gene$ensembl_gene_id\n",
    "\n",
    "#We will recursively perform differential expression analysis between the controls and each of the pseudotime samples\n",
    "pseudotime <- c(\"1\", \"2\", \"3\", \"4\", \"5\", \"6\", \"7\")\n",
    "output_folder <- folder\n",
    "\n",
    "for (i in 1:length(pseudotime)) {\n",
    "  \n",
    "  #Filter info data for desired pseudotime\n",
    "  col_data <- subset(covid_info_file, covid_info_file$Pseudotime %in% c('0', pseudotime[i]))\n",
    "  rownames(col_data) <- col_data$RNAseqID\n",
    "  col_data$Pseudotime <- as.factor(as.character(col_data$Pseudotime))\n",
    "  \n",
    "  #Filter count data for desired pseudotime\n",
    "  selected_count_data <- count_data[, as.character(col_data$RNAseqID)]\n",
    "  \n",
    "  #Create deseq object\n",
    "  dds_counts <- DESeqDataSetFromMatrix(countData = selected_count_data, colData = col_data, design = ~ Pseudotime)\n",
    "  dds_counts <- dds_counts[ rowSums(counts(dds_counts)) > 1, ]\n",
    "  dds_counts <- estimateSizeFactors(dds_counts)\n",
    "  \n",
    "  #Run DESeq\n",
    "  dds <- DESeq(dds_counts, betaPrior = FALSE)\n",
    "  res <- results(dds, independentFiltering = TRUE, alpha = 0.05)\n",
    "  res_sorted <- res[order(res$padj), ]\n",
    "  \n",
    "  #Add gene symbols\n",
    "  res_sorted$gene <- unique_ensg2gene[rownames(res_sorted), ]$hgnc_symbol\n",
    "    \n",
    "  fwrite(as.data.table(res_sorted, keep.rownames=\"ensembl_gene_id\"), file = file.path(folder, paste(\"/DESeq2result_genenames_Control_vs_COVID_\", pseudotime[i], \".txt\", sep = '')), \n",
    "              sep=\"\\t\", quote=FALSE)  \n",
    "  #write.table(res_sorted, file = file.path(folder, paste(\"/DESeq2result_genenames_Control_vs_COVID_\", pseudotime[i], \".txt\", sep = '')), \n",
    "  #            sep=\"\\t\", quote=FALSE)\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd48c12f",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "You can load each table into the memory using the following commands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3640efe6",
   "metadata": {
    "deletable": false,
    "editable": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#change the 1 for a number from 1-7 to get the other results.\n",
    "res_sorted<-fread(paste0(folder,\"/DESeq2result_genenames_Control_vs_COVID_1.txt\"))\n",
    "res_sorted[1:20,]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75ede534",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "#### Tasks\n",
    "\n",
    "1. Try to find the gene with the strongest expression in the dataset (baseMean). The command ```res_sorted <- res[order(res$padj), ]``` has been used before to sort the table. Can you adapt the command so you get the strongest expressed genes? \n",
    "\n",
    "Hint: Type ```?base::order``` if you need information how to tweak the behavior of the function. Function ```head()``` is also helpful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58574efc",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "head(res_sorted[order(res_sorted$baseMean,decreasing = T)])\n",
    "# or\n",
    "tail(res_sorted[order(res_sorted$baseMean)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf5fd1b7",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# Run your solution here:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "861d94fe",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "2. What is the transcript coding for? Look it up for example here: <https://www.genecards.org>. \n",
    "\n",
    "3. Look at the level of expression of the top two genes. Explain your FastQC results from the beginning considering this observation.\n",
    "\n",
    "#### Advanced\n",
    "\n",
    "1. Explain the gene's effect on the resolution of the experiment. What could you do to improve the resolution?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12920b30",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "### Volcano plot\n",
    "\n",
    "A volcano plot is a visualization of the significance and signed fold-change of genes in a differential expression experiment. It somehow resembles a volcano eruption with the most significant genes as rocks tossed into the sky."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccf3a839",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot_data <- data.frame()\n",
    "\n",
    "#Read each deseq out file and filter for significant differentially expressed genes\n",
    "for (i in 1:length(pseudotime)) {\n",
    "  input_folder <- folder\n",
    "  result <- read.csv(file.path(folder, paste(\"DESeq2result_genenames_Control_vs_COVID_\", pseudotime[i], \".txt\", sep = '')), \n",
    "                     sep = '\\t')\n",
    "  result <- subset(result, result$padj < 0.95 & result$baseMean > 100 & abs(result$log2FoldChange) > 0.01)\n",
    "  result$pseudotime <- pseudotime[i]\n",
    "  result <- rownames_to_column(result, \"Gene_ID\")\n",
    "  plot_data <- rbind(plot_data, result)\n",
    "  \n",
    "}\n",
    "\n",
    "#Select interesting genes to label\n",
    "gene_to_label <- c(\"MAOA\", \"HTRA3\", \"ERG\", \"TMEM119\", \"INHBB\", \"PTX3\", \"FCGR1CP\", \"PDZK1IP1\", \"FHL2\", \"NRIR\", \"SLC4A1\", \"C1QB\", \"TESC\", \"MERTK\", \n",
    "                   \"IDO1\", \"P2RY14\", \"F13A1\", \"STAT1\", \"GATA2\", \"NFKBIA\", \"TREML1\", \"NCR1\", \"NLRP6\", \"NLRP3\", \"JAK3\", \"TXK\", \"IL5RA\", \"SATB1-AS1\", \n",
    "                   \"SHARPIN\", \"TNFRSF10C\", \"CXCR1\", \"ATF6\", \"IFI16\", \"AHR\", \"IL13RA1\", \"C5AR1\", \"CSF2RA\", \"FYN\", \"IGHG1\", \"JCHAIN\", \"TOP2A\", \n",
    "                   \"IFI27\", \"BIRC5\", \"IGHA1\", \"IL1R2\", \"DEFA4\", \"DEFA3\", \"LTF\", \"RORC\", \"IFIT1B\", \"AIM2\", \"IL24\", \"CD38\", \"TNFRSF13C\", \"ERBB2\", \n",
    "                   \"TCF7\", \"CIITA\", \"TGFBI\", \"SIGIRR\", \"IGHA2\", \"CDC6\", \"CHI3L1\", \"HLA-DPA1\", \"HLA-DMB\", \"CD4\", \"IL6ST\", \"CD274\", \"C3AR1\", \"ADAM17\", \n",
    "                   \"C1QA\", \"CARD11\", \"CD14\", \"BCL3\", \"TICAM1\")\n",
    "gene_to_label2 <- c(\"HLA-DMA\", \"HLA-DMB\", \"HLA-DOA\", \"HLA-DPB1\", \"HLA-DPB\", \"HLA-DQA1\", \"HLA-DQB1\", \"HLA-DR1\", \n",
    "                    \"HLA-DRB1\", \"HLA-DRB5\")\n",
    "\n",
    "plot_data$color <- ifelse(plot_data$log2FoldChange > 1 & plot_data$padj < 0.05, \"red\", \n",
    "                             ifelse(plot_data$log2FoldChange < -1 & plot_data$padj < 0.05, \"blue\", \n",
    "                                    ifelse(plot_data$padj < 0.05, \"black\", \"grey\")))\n",
    "\n",
    "volcano1 <- ggplot(data = plot_data[plot_data$pseudotime==\"1\",], mapping=aes(x= log2FoldChange, y=-log10(padj), color=color)) + \n",
    "  geom_point() +\n",
    "  geom_label_repel(data = subset(plot_data, plot_data$gene %in% gene_to_label & abs(plot_data$log2FoldChange) > 2 & pseudotime==\"1\"), \n",
    "                   aes(x=log2FoldChange, y=-log10(padj), label=gene), \n",
    "                   color=\"red\", nudge_x = -0.5, min.segment.length = 0) + \n",
    "  geom_label_repel(data = subset(plot_data, plot_data$gene %in% gene_to_label2& abs(plot_data$log2FoldChange) > 1 & pseudotime == \"1\"), \n",
    "                   aes(x=log2FoldChange, y= -log10(padj), label=gene), \n",
    "                   color=\"red\", nudge_x = -0.5, min.segment.length = 0)\n",
    "volcano1 <- volcano1 + geom_vline(xintercept = 0, lty=3) + geom_vline(xintercept = 1, lty=4) + geom_vline(xintercept = -1, lty=4) + \n",
    "  geom_hline(yintercept = 1.30103, lty=4) + geom_hline(yintercept = 3, lty=4)\n",
    "volcano1 <- volcano1 + scale_color_manual(values = c('#606060','#004C99','#C0C0C0','#990000'))\n",
    "volcano1 <- volcano1 + theme_bw() + theme(axis.text = element_text(size = 20), axis.title = element_text(size = 22))\n",
    "volcano1\n",
    "ggsave(filename=paste0(folder,\"/volcano.png\"),plot=volcano1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e57e8d4",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "#### Tasks\n",
    "\n",
    "1. Change the code so you see the volcano plot of pseudotime 2 compared to controls."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d109277",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "#change the filter in the following line, which selects the data: plot_data[plot_data$pseudotime==\"2\",]\n",
    "#this selects the data at pseudo-timepoint 2\n",
    "#volcano1 <- ggplot(data = plot_data[plot_data$pseudotime==\"2\",], mapping=aes(x= log2FoldChange, y=-log10(padj), color=color))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3240109b",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Insert your code here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6eb8bb51",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# change the pseudotime==\"1\" to pseudotime==\"2\" et voilà.  \n",
    "volcano1 <- ggplot(data = plot_data[plot_data$pseudotime==\"2\",], mapping=aes(x= log2FoldChange, y=-log10(padj), color=color)) +\n",
    "  geom_point() +\n",
    "  geom_label_repel(data = subset(plot_data, plot_data$gene %in% gene_to_label & abs(plot_data$log2FoldChange) > 2 & pseudotime==\"2\"), \n",
    "                   aes(x=log2FoldChange, y=-log10(padj), label=gene), \n",
    "                   color=\"red\", nudge_x = -0.5, min.segment.length = 0) + \n",
    "  geom_label_repel(data = subset(plot_data, plot_data$gene %in% gene_to_label2& abs(plot_data$log2FoldChange) > 1 & pseudotime == \"2\"), \n",
    "                   aes(x=log2FoldChange, y= -log10(padj), label=gene), \n",
    "                   color=\"red\", nudge_x = -0.5, min.segment.length = 0)\n",
    "volcano1 <- volcano1 + geom_vline(xintercept = 0, lty=3) + geom_vline(xintercept = 1, lty=4) + geom_vline(xintercept = -1, lty=4) + \n",
    "  geom_hline(yintercept = 1.30103, lty=4) + geom_hline(yintercept = 3, lty=4)\n",
    "volcano1 <- volcano1 + scale_color_manual(values = c('#606060','#004C99','#C0C0C0','#990000'))\n",
    "volcano1 <- volcano1 + theme_bw() + theme(axis.text = element_text(size = 20), axis.title = element_text(size = 22))\n",
    "volcano1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f07a8286",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "2. Change the code so you only see genes which have a negative log2foldchange."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9e60fec",
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "# in the line below, add a filtering condition log2FoldChange < 0 to the data and connect the two filters with the & -sign.\n",
    "#volcano1 <- ggplot(data = plot_data[plot_data$pseudotime==\"1\" & log2FoldChange < 0,], mapping=aes(x= log2FoldChange, y=-log10(padj), color=color))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caff2d67",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Insert your code here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5138f308",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# abs() gives the absolute value of the input. We remove the function and set the comparison to \"plot_data$log2FoldChange < 0\"\n",
    "# For the labeling of significant genes, we can filter by removing abs() and changing to \"plot_data$log2FoldChange < -2\"\n",
    "volcano1 <- ggplot(data = plot_data[plot_data$pseudotime==\"1\" & plot_data$log2FoldChange < 0,], mapping=aes(x= log2FoldChange, y=-log10(padj), color=color)) +\n",
    "  geom_point() +\n",
    "  geom_label_repel(data = subset(plot_data, plot_data$gene %in% gene_to_label & plot_data$log2FoldChange < -2 & pseudotime==\"1\"), \n",
    "                   aes(x=log2FoldChange, y=-log10(padj), label=gene), \n",
    "                   color=\"red\", nudge_x = -0.5, min.segment.length = 0) + \n",
    "  geom_label_repel(data = subset(plot_data, plot_data$gene %in% gene_to_label2 & plot_data$log2FoldChange < -1 & pseudotime == \"1\"), \n",
    "                   aes(x=log2FoldChange, y= -log10(padj), label=gene), \n",
    "                   color=\"red\", nudge_x = -0.5, min.segment.length = 0)\n",
    "volcano1 <- volcano1 + geom_vline(xintercept = 0, lty=3) + geom_vline(xintercept = 1, lty=4) + geom_vline(xintercept = -1, lty=4) + \n",
    "  geom_hline(yintercept = 1.30103, lty=4) + geom_hline(yintercept = 3, lty=4)\n",
    "volcano1 <- volcano1 + scale_color_manual(values = c('#606060','#004C99','#C0C0C0','#990000'))\n",
    "volcano1 <- volcano1 + theme_bw() + theme(axis.text = element_text(size = 20), axis.title = element_text(size = 22))\n",
    "volcano1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eda2686d",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "## Gene expression heatmaps\n",
    "\n",
    "Heatmaps are an effective way to visualize gene expression of multiple genes across different samples or groups. The ComplexHeatmap package <https://jokergoo.github.io/ComplexHeatmap-reference/book/> provides multiple ways to add annotations to the heatmap as well as clustering options. Normalized genes expression data should be used to create Heatmaps. Here we plot a heatmap for the 100 most differentially expressed genes grouped by different disease severities (pseudotimes). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1033ac96",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "top100 <- res_sorted[1:100,]\n",
    "top100 <- subset(top100, top100$gene != '')\n",
    "\n",
    "col_data <- covid_info_file\n",
    "rownames(col_data) <- col_data$RNAseqID\n",
    "col_data$Gender <- as.factor(col_data$Gender)\n",
    "col_data$Disease.traj <- as.factor(col_data$Disease.traj)\n",
    "col_data$Pseudotime <- as.factor(as.character(col_data$Pseudotime))\n",
    "col_data$Sample_name <- paste(col_data$Name, col_data$Pseudotime, sep = '_')\n",
    "\n",
    "#Normalize the expression counts\n",
    "dds_counts <- DESeqDataSetFromMatrix(countData = count_data, colData = col_data, design = ~ Pseudotime)\n",
    "dds_counts <- dds_counts[ rowSums(counts(dds_counts)) > 1, ]\n",
    "dds_counts <- estimateSizeFactors(dds_counts)\n",
    "normalized_counts <- counts(dds_counts, normalized = TRUE)\n",
    "\n",
    "top100_expression_data <- normalized_counts[as.numeric(rownames(top100)), ]\n",
    "colnames(top100_expression_data) <- col_data$Sample_name\n",
    "\n",
    "#Calculate mean expression level for each pseudotime\n",
    "vecPT <- str_split_fixed(colnames(top100_expression_data), '_', 4)[, 4]\n",
    "vecUniquePT <- unique(vecPT)\n",
    "\n",
    "matDataHeat <- do.call(cbind, lapply(vecUniquePT, function(pt) {\n",
    "  vecidxCols <- which(vecPT %in% pt)\n",
    "  if (length(vecidxCols) > 1) {\n",
    "    return(rowMeans(top100_expression_data[, vecidxCols], na.rm = TRUE))\n",
    "  } else {\n",
    "    return(top100_expression_data[, vecidxCols])\n",
    "  }\n",
    "}))\n",
    "colnames(matDataHeat) <- vecUniquePT\n",
    "rownames(matDataHeat) <- top100$gene\n",
    "\n",
    "matDataHeat <- matDataHeat[, c(\"0\", \"1\", \"2\", \"3\", \"4\", \"5\", \"6\", \"7\")]\n",
    "\n",
    "#Scale the data row-wise to highlight the genes expression differences across pseudotimes\n",
    "matDataHeat <- t(scale(t(matDataHeat)))\n",
    "names(pseudotime_colors) <- colnames(matDataHeat)\n",
    "\n",
    "#Annotate the columns by pseudotimes and plot heatmap\n",
    "col_ha = HeatmapAnnotation(df = data.frame(Pseudotime = colnames(matDataHeat)),\n",
    "                           col = list(Pseudotime = pseudotime_colors), \n",
    "                           gp = gpar(col = \"black\"))\n",
    "Heatmap(matDataHeat, col=rev(brewer.pal(11, \"RdYlBu\")), top_annotation = col_ha, cluster_columns = FALSE, row_names_gp = gpar(fontsize = 8), border = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdda2d32",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "#### Tasks\n",
    "\n",
    "In the code for the above graph, you find this line ```matDataHeat <- t(scale(t(matDataHeat)))```. It performs a transformation on the averaged count data to be plotted. \n",
    "1. Find out what it does, try out and explain why we need it for an informative figure output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65998500",
   "metadata": {
    "solution2": "hidden"
   },
   "outputs": [],
   "source": [
    "#You can use the ?-sign to find out what a function does.\n",
    "?scale"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68c35766",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden"
   },
   "source": [
    "The data is scaled on gene level. This means that the mean value is subtracted and the difference is divided by the standard deviation of the expression values. By commenting out the line, you see that the expression level then dominates the clustering of the genes, not similarity in expression pattern. Scaling is a useful transformation if features to be compared have widely different means."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa5538d0",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "#### Expert task\n",
    "\n",
    "1. Plot a heatmap with all individual samples grouped by different pseudotimes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90bb3b98",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#order the sample data by pseudotime\n",
    "col_data_copy <- copy(col_data)\n",
    "setorder(col_data_copy, Pseudotime)\n",
    "#reorder the columns of the expression data accordingly\n",
    "top100_expression_data_sorted <- top100_expression_data[,col_data_copy$Sample_name]\n",
    "#give gene names as rownames, overwriting the ENSG... naming\n",
    "rownames(top100_expression_data_sorted) <- top100$gene\n",
    "#scale the expression data on the gene level, to make them comparable across samples\n",
    "top100_expression_data_sorted <- t(scale(t(top100_expression_data_sorted)))\n",
    "#prepare the annotation\n",
    "name_colors <- col_data_copy$color_code\n",
    "names(name_colors) <- col_data_copy$Sample_name\n",
    "names(pseudotime_colors) <- colnames(matDataHeat)\n",
    "col_ha2 <- HeatmapAnnotation(df = data.frame(#Name = col_data$Sample_name,\n",
    "                           Pseudotime = col_data_copy$Pseudotime),\n",
    "                           col = list(Pseudotime = pseudotime_colors), \n",
    "                           gp = gpar(col = \"black\"))\n",
    "Heatmap(top100_expression_data_sorted, col=rev(brewer.pal(11, \"RdYlBu\")), top_annotation = col_ha2, cluster_columns = FALSE, row_names_gp = gpar(fontsize = 8), border = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bcc7bde",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "## Enrichment analysis of differentially expressed genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26127f45",
   "metadata": {
    "deletable": false,
    "editable": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "res_sorted <- read.table(paste0(folder,\"/DESeq2result_genenames_Control_vs_COVID_1.txt\"), header = TRUE, sep = '\\t')\n",
    "write(res_sorted$gene[1:100],file = paste0(folder, \"/top100genes.txt\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d4df188",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "You now have a list of significant differentially expressed genes in patients with complicated incremental disease. You can use <https://maayanlab.cloud/Enrichr/> to check if you can find an enriched term that could give you insight on the underlying processes. Tip: Check the \"Pathways\" and \"Ontologies\" tabs. Click on the red captions of the database name. You get to a view where you can hover over the bars to check the q-value of the enrichment. The q-value is a multiple testing corrected p-value and should be less than 0.05 to speak of a significant enrichment.\n",
    "\n",
    "#### Tasks\n",
    "\n",
    "Do you find a term which seems related to COVID? Which one?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e841886",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "### TopGO\n",
    "\n",
    "As more accurate approach, you can use a bioconductor package called topGO to perform the GO Term enrichment analysis. topGO package provides tools for testing GO terms while accounting for the topology of the GO graph. It can also account for the gene background (genes that were detected in the experiment at least once). Different test statistics and different methods for eliminating local similarities and dependencies between GO terms can be implemented and applied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b702874",
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "rownames(res_sorted) <- res_sorted$ensembl_gene_id\n",
    "#filter out genes that were never detected\n",
    "res_sorted <- subset(res_sorted, res_sorted$baseMean > 0)\n",
    "degs <- subset(res_sorted, res_sorted$padj < 0.05 & abs(res_sorted$log2FoldChange) > 1)\n",
    "top100 <- rownames(degs[1:100,])\n",
    "all_genes <- rownames(res_sorted)\n",
    "#top100 <- degs$ensembl_gene_id[1:100]\n",
    "#all_genes <- res_sorted$ensembl_gene_id\n",
    "\n",
    "inUniverse <- all_genes %in% all_genes\n",
    "inSelection <- all_genes %in% top100 \n",
    "alg <- factor( as.integer( inSelection[inUniverse] ) )\n",
    "names(alg) <- all_genes[inUniverse]\n",
    "\n",
    "## prepare data\n",
    "tgd <- new( \"topGOdata\", ontology=\"BP\", allGenes = alg, nodeSize=5,\n",
    "            annot=annFUN.org, mapping=\"org.Hs.eg.db\", ID = \"ensembl\" )\n",
    "\n",
    "## run tests\n",
    "resultTopGO.elim <- runTest(tgd, algorithm = \"elim\", statistic = \"Fisher\" )\n",
    "resultTopGO.classic <- runTest(tgd, algorithm = \"classic\", statistic = \"Fisher\" )\n",
    "\n",
    "## look at results\n",
    "if(length(nodes(graph(tgd))) < 200){\n",
    "  tab <- GenTable( tgd, Fisher.elim = resultTopGO.elim, \n",
    "                   Fisher.classic = resultTopGO.classic,\n",
    "                   orderBy = \"Fisher.elim\" , topNodes = length(nodes(graph(tgd))))\n",
    "}else{\n",
    "  tab <- GenTable( tgd, Fisher.elim = resultTopGO.elim, \n",
    "                   Fisher.classic = resultTopGO.classic,\n",
    "                   orderBy = \"Fisher.elim\" , topNodes = 200)}\n",
    "\n",
    "printGraph(tgd, resultTopGO.elim, firstSigNodes = 5, fn.prefix = paste0(folder,\"/\"), useInfo = \"all\", pdfSW = TRUE)\n",
    "printGraph(tgd, resultTopGO.classic, firstSigNodes = 15, fn.prefix = paste0(folder,\"/\"), useInfo = \"all\", pdfSW = TRUE)\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ea5381e",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "#### Tasks\n",
    "\n",
    "Do you see any biological processes that are related to the immune response or inflammation?\n",
    "1. Change the code to identify the GO terms that describe molecular function (MF) and/or cellular components (CC) instead of biological processes (BP). Hint: change 'ontology' while creating the topGOdata object.\n",
    "2. Change the code to identify GO terms that are enriched in genes that are downregulated in COVID-19."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4f1a6ef",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden",
    "solution2_first": true
   },
   "source": [
    "#### Expert questions\n",
    "\n",
    "1. Can you see any difference in the p-values generated by the Fisher classic algorithm and Fisher elim algorithm? Check the graphs for the difference in p-values at the different levels of the GO graph.\n",
    "2. Is it correct to calculate the enrichment of GO terms for our genes of interest against all expressed genes? Might it be better if we compare genes that have similar expression?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cd5ae85",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden"
   },
   "source": [
    "Run the code and find out!\n",
    "\n",
    "The function genefinder from the genefilter package will be used to find background genes that are similar in expression to our genes of interest. The function tries to identify 10 genes for each gene of interest that match its expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3517bcf6",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "overallBaseMean <- as.matrix(res_sorted[, c(\"baseMean\"), drop = F])\n",
    "\n",
    "sig_idx <- match(top100, rownames(overallBaseMean))\n",
    "\n",
    "backG <- c()\n",
    "\n",
    "for(i in sig_idx){\n",
    "  ind <- genefinder(overallBaseMean, i, 10, method = \"manhattan\")[[1]]$indices\n",
    "  backG <- c(backG, ind)\n",
    "  \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a26fbf0",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden"
   },
   "source": [
    "We then check whether the background has roughly the same distribution of average expression strength as the foreground by plotting the densities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59fb4b41",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "backG <- unique(backG)\n",
    "backG <- rownames(overallBaseMean)[backG]\n",
    "backG <- setdiff(backG, top100)\n",
    "\n",
    "multidensity( list( \n",
    "  all= log2(res_sorted[,\"baseMean\"]) ,\n",
    "  foreground =log2(res_sorted[top100, \"baseMean\"]), \n",
    "  background =log2(res_sorted[backG, \"baseMean\"])), \n",
    "  xlab=\"log2 mean normalized counts\", main = \"Matching for enrichment analysis\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "289e056a",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden"
   },
   "source": [
    "Run topGO with unbiased background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45ab6391",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "inUniverse <- all_genes %in% c(top100, backG)\n",
    "inSelection <- all_genes %in% top100 \n",
    "alg <- factor( as.integer( inSelection[inUniverse] ) )\n",
    "names(alg) <- all_genes[inUniverse]\n",
    "\n",
    "## prepare data\n",
    "tgd <- new( \"topGOdata\", ontology=\"BP\", allGenes = alg, nodeSize=5,\n",
    "            annot=annFUN.org, mapping=\"org.Hs.eg.db\", ID = \"ensembl\" )\n",
    "\n",
    "## run tests\n",
    "resultTopGO.elim <- runTest(tgd, algorithm = \"elim\", statistic = \"Fisher\" )\n",
    "resultTopGO.classic <- runTest(tgd, algorithm = \"classic\", statistic = \"Fisher\" )\n",
    "\n",
    "## look at results\n",
    "if(length(nodes(graph(tgd))) < 200){\n",
    "  tab <- GenTable( tgd, Fisher.elim = resultTopGO.elim, \n",
    "                   Fisher.classic = resultTopGO.classic,\n",
    "                   orderBy = \"Fisher.elim\" , topNodes = length(nodes(graph(tgd))))\n",
    "}else{\n",
    "  tab <- GenTable( tgd, Fisher.elim = resultTopGO.elim, \n",
    "                   Fisher.classic = resultTopGO.classic,\n",
    "                   orderBy = \"Fisher.elim\" , topNodes = 200)}\n",
    "\n",
    "printGraph(tgd, resultTopGO.elim, firstSigNodes = 5, fn.prefix = paste0(folder,\"/selected_background\"), useInfo = \"all\", pdfSW = TRUE)\n",
    "printGraph(tgd, resultTopGO.classic, firstSigNodes = 15, fn.prefix = paste0(folder,\"/selected_background\"), useInfo = \"all\", pdfSW = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c28537f5",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    },
    "solution2": "hidden"
   },
   "source": [
    "We next find the genes associated with the enriched GO terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69a2c806",
   "metadata": {
    "deletable": false,
    "editable": false,
    "solution2": "hidden",
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "for(i in 1:nrow(tab)){\n",
    "  go_id <- as.vector(tab[i,1])\n",
    "  genes_in_cat <- intersect(unlist(genesInTerm(tgd, go_id)), top100)\n",
    "  gene_sym_in_cat <- unique_ensg2gene[as.character(genes_in_cat), \"hgnc_symbol\"]\n",
    "  gene_sym_in_cat_str <- \"\"\n",
    "  \n",
    "  if(length(genes_in_cat) > 0){\n",
    "    for(j in 1:length(gene_sym_in_cat)){\n",
    "      gene_sym_in_cat_str <- paste(gene_sym_in_cat_str, gene_sym_in_cat[j], sep = ',')\n",
    "    }\n",
    "  }\n",
    "  if (gene_sym_in_cat_str == \"\"){\n",
    "    gene_sym_in_cat_str <- NA\n",
    "  }\n",
    "\n",
    "  tab$Genes[i] <- gene_sym_in_cat_str\n",
    "\n",
    "}\n",
    "tab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d97b1e1",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "## Misc\n",
    "\n",
    "This code was used to download the ensemblID to gene name matching table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e3cb663",
   "metadata": {
    "deletable": false,
    "editable": false,
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# library('biomaRt')\n",
    "# ensembl = useEnsembl(biomart = \"genes\")\n",
    "# ensembl = useDataset(\"hsapiens_gene_ensembl\",mart=ensembl)\n",
    "# #ensembl = useMart(\"ensembl\", dataset = \"hsapiens_gene_ensembl\")\n",
    "# genes <- rownames(count_data)\n",
    "# #df<-df[,-4]\n",
    "# G_list <- getBM(filters= \"ensembl_gene_id\", attributes= c(\"ensembl_gene_id\",\"hgnc_symbol\"),values=genes,mart= ensembl)\n",
    "# G_DT<-as.data.table(G_list)\n",
    "# data.table::fwrite(G_DT, col.names = T,file = \"ensg2gene.txt\")\n",
    "# merge(df,G_list,by.x=\"gene\",by.y=\"ensembl_peptide_id\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "rnaseq_microbiome",
   "language": "R",
   "name": "rnaseq_microbiome"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.1.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}