workflow.Rmd

---
title: "Example workflow for bulk RNA-seq analysis"
author: "Mark Ziemann"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
    fig_width: 7
    fig_height: 7
theme: cosmo
---

Source: https://github.com/markziemann/background

## Intro

Here we are performing an analysis of some gene expression data to demonstrate the difference between ORA and FCS methods and to highlight the differences caused by improper background gene set use.

The dataset being used is SRP128998 and we are comparing the cells grown in normal glucose condition (control) to the high glucose condition (case).

Data are obtained from http://dee2.io/

```{r,begin}

suppressPackageStartupMessages({
  library("kableExtra")
  library("eulerr")
  library("gplots")
  library("getDEE2")
  library("DESeq2")
})

```

## Get expression data and make an MDS plot

```{r,getdata}

name = "SRP128998"
mdat <- getDEE2Metadata("hsapiens")
samplesheet <- mdat[grep("SRP128998",mdat$SRP_accession),]
samplesheet <- samplesheet[order(samplesheet$SRR_accession),]
samplesheet$trt <- as.factor(c(1,1,1,1,1,1,0,0,0,0,0,0))
samplesheet$VPA <- as.factor(c(0,0,0,1,1,1,0,0,0,1,1,1))
s1 <- subset(samplesheet,VPA==0)

s1 %>%
  kbl(caption = "sample sheet") %>%
  kable_paper("hover", full_width = F)

w <- getDEE2("hsapiens", samplesheet$SRR_accession,
  metadata=mdat,legacy = TRUE)

x <- Tx2Gene(w)
x <- x$Tx2Gene

# table of gene symbols
gt <- w$GeneInfo[,1,drop=FALSE]
gt$accession <- rownames(gt)

# fix gene symbols
rownames(x) <- sapply(strsplit(rownames(x),"\\."),"[[",1)
x  <- merge(gt,x,by=0)
rownames(x) <- paste(x$Row.names, x$GeneSymbol)
x <- x[,-c(1:3)]

# counts
x1 <- x[,which(colnames(x) %in% s1$SRR_accession)]

colnames(x1) <- c("HG1","HG2","HG3","NG1","NG2","NG3")

head(x1) %>%
  kbl(caption = "counts") %>%
  kable_paper("hover", full_width = F)

```

Here show the number of genes in the annotation set, and those detected above the detection threshold.

```{r,filter}

# filter out lowly expressed genes
x1 <- x1[which(rowSums(x1)/ncol(x1)>=(10)),]
nrow(x)
nrow(x1)

```

Now multidimensional scaling (MDS) plot to show the correlation between the datasets.
If the control and case datasets are clustered separately, then it is likely that there will be many differentially expressed genes with FDR<0.05.

```{r,mds}

plot(cmdscale(dist(t(x1))), xlab="Coordinate 1", ylab="Coordinate 2", pch=19, col=s1$trt, main="MDS")

```

## Differential expression

Now run DESeq2 for control vs case.

```{r,deseq2}

y <- DESeqDataSetFromMatrix(countData = round(x1), colData = s1, design = ~ trt)
y <- DESeq(y)
de <- results(y)
de <- as.data.frame(de[order(de$pvalue),])
#rownames(de) <- sapply(strsplit(rownames(de),"\\."),"[[",1)
head(de) %>% kbl() %>% kable_paper("hover", full_width = F)

```

Now let's have a look at some of the charts showing differential expression.
In particular, an MA plot and volcano plot.

```{r,deplots,fig.width=7,fig.height=7}

maplot <- function(de,contrast_name) {
  sig <-subset(de, padj < 0.05 )
  up <-rownames(subset(de, padj < 0.05 & log2FoldChange > 0))
  dn <-rownames(subset(de, padj < 0.05 & log2FoldChange < 0))
  GENESUP <- length(up)
  GENESDN <- length(dn)
  DET=nrow(de)
  SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down", DET, "detected")
  ns <-subset(de, padj > 0.05 )
  plot(log2(de$baseMean),de$log2FoldChange, 
       xlab="log2 basemean", ylab="log2 foldchange",
       pch=19, cex=0.5, col="dark gray",
       main=contrast_name, cex.main=0.7)
  points(log2(sig$baseMean),sig$log2FoldChange,
         pch=19, cex=0.5, col="red")
  mtext(SUBHEADER,cex = 0.7)
}

make_volcano <- function(de,name) {
    sig <- subset(de,padj<0.05)
    N_SIG=nrow(sig)
    N_UP=nrow(subset(sig,log2FoldChange>0))
    N_DN=nrow(subset(sig,log2FoldChange<0))
    DET=nrow(de)
    HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
    plot(de$log2FoldChange,-log10(de$padj),cex=0.5,pch=19,col="darkgray",
        main=name, xlab="log2 FC", ylab="-log10 pval", xlim=c(-6,6))
    mtext(HEADER)
    grid()
    points(sig$log2FoldChange,-log10(sig$padj),cex=0.5,pch=19,col="red")
}

maplot(de,name)
make_volcano(de,name)

```

Heatmap of top genes.

```{r,heat,fig.width=7,fig.height=7}

colfunc <- colorRampPalette(c("blue", "white", "red"))
topgenes <- rownames(head(de,30))
rpm <- apply(x1,2,function(x) {x / sum(x) * 1e6 } )
top <- rpm[which(rownames(rpm) %in% topgenes),]
heatmap.2(top,trace="none",col=colfunc(25),scale="row", margins = c(10,20))

```

## Save table

```{r,save}

saveRDS(de,"workflow.Rds")

```

## Session information

```{r,session}

sessionInfo()

```