workshop-example.txt

---
title: "workshop-example"
author: "HBC Training Team"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Project details
For this analysis, we are using the full count matrix from an RNA-Seq dataset which is part of a larger study described in [Kenny PJ et al, Cell Rep 2014]().

The RNA-Seq was performed on HEK293F cells that were either transfected with a MOV10 transgene, or siRNA to knock down Mov10 expression, or non-specific (irrelevant) siRNA. This resulted in 3 conditions Mov10 oe (over expression), Mov10 kd (knock down) and Irrelevant kd, respectively. The number of replicates is as shown below.

Using these data, we will evaluate transcriptional patterns associated with perturbation of MOV10 expression. Please note that the irrelevant siRNA will be treated as our control condition.

## Setup

### Load Libraries

```{r load-libraries, message=FALSE}

### Bioconductor and CRAN libraries used
library(ggplot2)
library(RColorBrewer)
library(DESeq2)
library(pheatmap)
library(grid)
library(gridExtra)
library(RColorBrewer)
library(knitr)
```

### Load Data

```{r load_data, message=FALSE}
data <- read.table("data/Mov10_full_counts.txt", header=T, row.names=1) 
meta <- read.table("meta/Mov10_full_meta.txt", header=T, row.names=1)
kable(meta, format="markdown", row.names=F)
```

## DESeq2 Analysis

```{r DE_analysis, message=FALSE}
dds <- DESeqDataSetFromMatrix(countData = data, colData = meta, design = ~ sampletype)
dds <- DESeq(dds)
```

## Quality Control

Both sets of QC plots show the samples from the same group clustering together, this means that the data look good to move forward with the differential gene expression analysis.

```{r qc, message=FALSE}
### Transform counts for data visualization
rld <- rlog(dds, blind=TRUE)

### Plot PCA 
plotPCA(rld, intgroup="sampletype")

### Extract the rlog matrix and compute correlation
rld_mat <- assay(rld) 
rld_cor <- cor(rld_mat) 

### Plot heatmap
pheatmap(rld_cor)
```


## Extracting Results

```{r extract_results, message=FALSE}

## Extract results
res_tableOE <- results(dds, contrast = c("sampletype", "MOV10_overexpression", "control"))
res_tableKD <- results(dds, contrast = c("sampletype", "MOV10_knockdown", "control"))

## Set thresholds
padj.cutoff <- 0.05
lfc.cutoff <- 0.58

# Significant results for mov10OE relative to control
sigOE <- subset(data.frame(res_tableOE), padj < padj.cutoff & abs(log2FoldChange) > lfc.cutoff)

# Significant results for mov10KD relative to control
sigKD <- subset(data.frame(res_tableKD), padj < padj.cutoff & abs(log2FoldChange) > lfc.cutoff)

```

## Summarizing and Visualizing Results

At an adjusted p-value cutoff of `r padj.cutoff` and an absolute log2FoldChange cutoff of `r lfc.cutoff` we find the follwing number of significant genes:

* Mov10 overepression vs. Control: **`r nrow(sigOE)`** genes
* Mov10 knockdown vs. Control: **`r nrow(sigKD)`** genes

## Volcano plots

```{r volcano_plot, message=FALSE}

# Identify significant OE genes
res_tableOE$threshold <- res_tableOE$padj < padj.cutoff & abs(res_tableOE$log2FoldChange) > lfc.cutoff

# Convert the DESeq results object to a data frame
resOE_df <- data.frame(res_tableOE)

# Volcano plot
p1 <- ggplot(resOE_df) +
  geom_point(aes(x=log2FoldChange, y=-log10(padj), colour=threshold)) +
  ggtitle("Mov10 overexpression") +
  xlab("log2 fold change") + 
  ylab("-log10 adjusted p-value") +
  theme(legend.position = "none",
        plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = rel(1.25))) 

# Identify significant KD genes
res_tableKD$threshold <- res_tableKD$padj < padj.cutoff & abs(res_tableKD$log2FoldChange) > lfc.cutoff

# Convert the DESeq results object to a data frame
resKD_df <- data.frame(res_tableKD)

# Volcano plot
p2 <- ggplot(resKD_df) +
  geom_point(aes(x=log2FoldChange, y=-log10(padj), colour=threshold)) +
  ggtitle("Mov10 knockdown") +
  xlab("log2 fold change") + 
  ylab("-log10 adjusted p-value") +
  theme(legend.position = "none",
        plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = rel(1.25))) 

grid.arrange(p1, p2, nrow=1)
```

## Heatmap of significant genes

```{r heatmap, message=FALSE}

## normalized counts 
normalized_counts <- counts(dds, normalized=T)

### Extract normalized expression for significant genes
norm_OEsig <- normalized_counts[rownames(sigOE), which(meta$sampletype != "MOV10_knockdown")]
norm_KDsig <- normalized_counts[rownames(sigKD), which(meta$sampletype != "MOV10_overexpression")]

### Annotate our heatmap 
annotation_oe <- meta[which(meta$sampletype != "MOV10_knockdown"), 'sampletype', drop=F]
annotation_kd <- meta[which(meta$sampletype != "MOV10_overexpression"), 'sampletype', drop=F]

### Set a color palette
heat.colors <- brewer.pal(6, "YlOrRd")

### Run pheatmap
pheatmap(norm_OEsig, color = heat.colors, cluster_rows = T, show_rownames=F,
annotation= annotation_oe, border_color=NA, fontsize = 10, scale="row",
     fontsize_row = 10, height=20)

pheatmap(norm_KDsig, color = heat.colors, cluster_rows = T, show_rownames=F,
annotation= annotation_kd, border_color=NA, fontsize = 10, scale="row",
     fontsize_row = 10, height=20)
```