rna-differential-expression-testing.Rmd

---
title: "RNA-seq differential expression testing"
output: github_document
---

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

So far in this workshop, we have been working with a subset of RNA-seq data that originates from chromosome 22 genes in order for the data processing to be efficient. However, for this section, we can analyse the data for all genes that has been pre-processed and saved in this GitHub repository. 

```{r, message=FALSE}
library(magrittr)
library(stringr)
library(edgeR)
library(ggplot2)
library(pheatmap)
library(gprofiler2)
library(wiggleplotr)
library(GenomicFeatures)
library(ensembldb)
library(EnsDb.Hsapiens.v86)
library(biomaRt)
```

To get started, load the sample data which has sample IDs and the BAM file identifiers in our counts matrix
```{r, cache=TRUE}
mdat <- read.csv("metadata/sample_data.csv")
mdat
```

And then load the gene and transcript counts data
```{r, cache=TRUE}
gene_counts <- readRDS("processed_data/Effector_CD4pos_ALL_gene_counts.Rds") 
tx_counts <- readRDS("processed_data/Effector_CD4pos_ALL_transcript_counts.Rds")
```

Next, let's ensure metadata and data labels are in the same order
```{r, cache=TRUE}
if (all(mdat$Bam != gene_counts$targets)){
    ind <- match(gene_counts$targets, mdat$Bam)
    mdat <- mdat[ind, ]
}

stopifnot(all(mdat$Bam == gene_counts$targets))
```

Create the DGElist object for using edgeR. 
```{r, cache=TRUE}
y <- DGEList(counts = gene_counts$counts,
             samples = mdat$Sample,
             group = mdat$Treatment)
```

### Design matrix for differential testing
```{r, cache=TRUE}
design <- model.matrix(~ 0 + y$samples$group)

# Clean up the column names (will make things easier later)
colnames(design) <- c("Stimulated", "Unstimulated")
design
```

Boxplot of counts-per-million (CPM) normalised expression for each sample. This allows us to inspect the distribution of library size normalised gene expression values. 
```{r, cache=TRUE, fig.height=5}
boxplot(cpm(y, log=TRUE), names=y$samples$samples)
```

Next we filter genes/transcripts to remove features with low counts that would not be useful in statistical testing. edgeR has a handy function for this called `filterByExpr` that filters features that do not have sufficient counts for differential expression testing. 
```{r, cache=TRUE}
keep <- filterByExpr(y)
table(keep)

y <- y[keep, , keep.lib.sizes=FALSE]
```

Boxplots of log2 counts for each sample now the data have been filtered
```{r, fig.height=5, cache=TRUE}
boxplot(cpm(y, log=TRUE), names=y$samples$samples)
```

Calculate normalisation factors and dispersions used in differential expression testing
```{r, cache=TRUE}
y <- calcNormFactors(y)
y$samples
```
```{r, cache=TRUE}
y <- estimateDisp(y, robust = TRUE)
```

**Plot samples using multi-dimensional scaling** How do the samples group based on the genes with the largest differences? Do they group by sample donor or by treatment group? In this type of plot, you want to look for outliers. 

```{r, cache=TRUE}
plotMDS(y, labels = y$samples$samples)
plotMDS(y, labels = y$samples$group)
```

Differential testing using the exact test in edgeR
```{r, cache=TRUE}
et <- exactTest(y)
tt <- topTags(et, n = nrow(y))
```

How many DE genes are there with FDR < 0.05 and logFC > 1?
```{r, cache=TRUE}
sum((tt$table$FDR < 0.05) & (abs(tt$table$logFC) > 1))
```

### Differential testing using Limma-Voom
```{r, cache=TRUE}
v <- voom(counts = y, design = design, plot = TRUE)
```

See how the mean-variance trend looks for unfiltered count data
```{r}
v2 <- voom(counts = gene_counts$counts, design = design, plot = TRUE)
```



```{r, fig.height=5, cache=TRUE}
boxplot(v$E, names=v$targets$samples)
```

If you want to normalise further using quantile normalisation
```{r, cache=TRUE}
vq <- voom(counts = y, design = design, plot = TRUE,
           normalize.method = "quantile")
```

```{r, fig.height=5, cache=TRUE}
boxplot(vq$E)
```
Notice how the quantile normalisation made the distribtions of expression extimates for each library look more similar?

> **_Challenge:_** From here, try to compare the results of using Limma Voom with and without quantile normalisation and see how this impacts your downstream results.

### Differential expression testing

Fitting linear models in Limma
```{r, cache=TRUE}
fit <- lmFit(object = v, design = design)
```

Make the contrasts matrix
```{r, cache=TRUE}
contr <- makeContrasts(Stimulated - Unstimulated, levels = colnames(coef(fit)))
contr
```

Estimate contrasts for each gene
```{r, cache=TRUE}
fit2 <- contrasts.fit(fit, contr)
```

Empirical Bayes smoothing of standard errors (shrinks standard errors that are much larger or smaller than those from other genes towards the average standard error) (see https://www.degruyter.com/doi/10.2202/1544-6115.1027)
```{r, cache=TRUE}
eb <- eBayes(fit2)
```

Check out the top differential expressed genes
```{r, cache=TRUE}
topTable(eb)
```

Put all DE testing results into an object for plotting and saving
```{r, cache=TRUE}
tt2 <- topTable(eb, number = nrow(v))
tt2$Significant <- ((abs(tt2$logFC) > 1) & (tt2$adj.P.Val < 0.05))
```

How many DE genes are there?
```{r, cache=TRUE}
table(tt2$Significant)
```

### Differential expression summary plots

Volcano plot
```{r, cache=TRUE}
gg_volcano <- ggplot(tt2, aes(x = logFC, y = -log10(P.Value),
                             fill=Significant, colour=Significant)) +
    geom_point(alpha=0.5) +
    geom_vline(xintercept = c(-1, 1), linetype="dashed")
gg_volcano
```

MA plot
```{r, cache=TRUE}
gg_ma <- ggplot(tt2, aes(x = AveExpr, y = logFC, colour=Significant)) +
    geom_point(alpha=0.5) +
    geom_vline(xintercept = c(1), linetype="dashed") +
    geom_hline(yintercept = c(-1, 1), linetype="dashed")
gg_ma
```

### Significant gene plots

DE gene heatmap of DE genes
```{r, fig.height=1.5, fig.width=1, cache=TRUE}
top_de_cpm <- v$E[rownames(v$E) %in% rownames(tt2)[tt2$Significant == TRUE], ]
colnames(top_de_cpm) <- v$targets$samples

pheatmap(top_de_cpm, scale = "row", show_rownames = FALSE)
```

Plot normalised gene expression for top genes
```{r, cache=TRUE}
plot_gene <- function(gene_ids){
    
    # Get the normalised expression data for the selected genes
    dat <- v$E[rownames(v$E) %in% gene_ids, ] %>% 
        data.frame()
    colnames(dat) <- v$targets$samples
    
    # Re-shape the data into the ggplot preferred long format
    dat$gene_id <- rownames(dat)
    dat <- reshape2::melt(dat)
    
    # Add the group data
    dat$group <- v$targets$group[match(dat$variable, v$targets$samples)]
    
    # plot the gene expression vaues
    ggplot(data = dat, aes(x = group, y = value)) + 
        geom_point(size=3, alpha=0.5) +
        facet_wrap(~gene_id, scales = "free_y") +
        ylab("Normalised gene expression") +
        theme_bw()
        
}

my_top_genes <- rownames(tt2)[tt2$Significant == TRUE][1:6]

plot_gene(gene_ids = my_top_genes)

```

### Genome browser track plots

First, we download annotation data for our genes/transcripts from Ensembl
```{r, cache=TRUE}
ensembl_mart <- useMart("ENSEMBL_MART_ENSEMBL",
                        host = "jan2020.archive.ensembl.org")

ensembl_dataset <- useDataset("hsapiens_gene_ensembl", mart=ensembl_mart)
ensembl_dataset

attributes <- listAttributes(ensembl_dataset)
head(attributes)

selected_attributes <- c("ensembl_transcript_id", "ensembl_gene_id", 
                        "external_gene_name", "strand", 
                        "gene_biotype", "transcript_biotype",
                        "chromosome_name", "start_position",
                        "end_position")

bm_data <- getBM(attributes = selected_attributes,
                 mart = ensembl_dataset)
head(bm_data)

# rename the columns to suit the plotting functions
bm_data <- dplyr::rename(bm_data, 
                     transcript_id = ensembl_transcript_id, 
                     gene_id = ensembl_gene_id, 
                     gene_name = external_gene_name)
head(bm_data)
```

And annotate the DE gene table
```{r, cache=TRUE}
ind <- match(rownames(tt2), bm_data$gene_id)
tt_annotated <- cbind(tt2, bm_data[ind, ])

head(tt_annotated)
```


And get the postions of exons for plotting
```{r, cache=TRUE}
exons <- exonsBy(EnsDb.Hsapiens.v86, by = "gene")
```

Then create a table of our bigwig files that we will use for plotting. Remember, our bigwig files we created earlier were only for chromosome 22 mapped reads, so we will have to limit our plotting to chromosome 22 differential genes here. 
```{r, cache=TRUE}
bw_files <- list.files(path = "aligned_data/",
                       pattern = "bigwig$",
                       full.names = TRUE)

ids <- basename(bw_files) %>%
    str_remove("_filt_fastp_sorted.bigwig") %>%
    str_remove("Effector_CD4pos_")

condition <- ifelse(test = grepl(pattern = "_T-S_", x = bw_files),
                    yes = "Stimulated", "Unstimulated")

bw_df <- data.frame(sample_id=ids, track_id=ids,
                    scaling_factor=1, bigWig=bw_files, colour_group=condition)
bw_df
```

Get the top DE genes from chromosome 22
```{r, cache=TRUE}
chr22_de_genes <- tt_annotated$gene_id[tt_annotated$Significant == TRUE & tt_annotated$chromosome_name == "22"]
```

And then use this information to plot the coverage of differential expressed genes for each sample
```{r, cache=TRUE}
plotCoverage(exons = exons[chr22_de_genes[1]], track_data = bw_df,
             transcript_annotations = bm_data, 
             rescale_introns = TRUE)
```

```{r, cache=TRUE}
plotCoverage(exons = exons[chr22_de_genes[3]], track_data = bw_df,
             transcript_annotations = bm_data,
             rescale_introns = FALSE)
```

### Ontology testing

Finally, we can use the `gprofiler` package to test for enrichment of gene ontology terms (among other things) in our differential expressed gene list as follows
```{r, cache=TRUE}
ontology_result <- gprofiler2::gost(query = rownames(tt2)[tt2$Significant])
head(ontology_result$result)
```

Plot ontology result
```{r, cache=TRUE}
gprofiler2::gostplot(ontology_result)
```

===

```{r}
sessionInfo()
```