airway_salmon_edgeR.Rmd

---
title: "Airway: DE analysis based on transcript counts and DESeq2"
author: "Lieven Clement"
date: "statOmics, Ghent University (https://statomics.github.io)"
output:
    html_document:
      code_download: true
      theme: cosmo
      toc: true
      toc_float: true
      highlight: tango
      number_sections: true
---

# Background
The data used in this workflow comes from an RNA-seq experiment where airway smooth muscle cells were treated with dexamethasone, a synthetic glucocorticoid steroid with anti-inflammatory effects (Himes et al. 2014). Glucocorticoids are used, for example, by people with asthma to reduce inflammation of the airways. In the experiment, four human airway smooth muscle cell lines were treated with 1 micromolar dexamethasone for 18 hours. For each of the four cell lines, we have a treated and an untreated sample.
For more description of the experiment see the article, PubMed entry 24926665, and for raw data see the GEO entry GSE52778.

Many parts of this tutorial are based on parts of a published RNA-seq workflow
available via Love et al. 2015 [F1000Research](http://f1000research.com/articles/4-1070) and as a [Bioconductor
package](https://www.bioconductor.org/packages/release/workflows/html/rnaseqGene.html) and on Charlotte Soneson's material from the [bss2019](https://uclouvain-cbio.github.io/BSS2019/rnaseq_gene_summerschool_belgium_2019.html) workshop.

# Data

We will use the salmon quantification files here. This is a lightweight transcript level mapper. 

## Meta data

As before, we will retrieve the meta data by linking the SRA files to the meta data on GEO. 

```{r}
library(edgeR)
library(tidyverse)
library(GEOquery)

download.file("https://github.com/statOmics/SGA/archive/airwaySeqData.zip","SGA-airwaySeqData.zip")
unzip("SGA-airwaySeqData.zip", exdir = "./")
```

### GEO data

```{r}
gse<-getGEO("GSE52778")
pdata<-pData(gse[[1]])
pdata$SampleName<-rownames(pdata) %>% as.factor
head(pdata)[1:6,]
```
### SRA info 
Download SRA info. To link sample info to info sequencing: Go to corresponding SRA page and save the information via the “Send to: File button” This file can also be used to make a script to download sequencing files from the web. Note that sra files can be converted to fastq files via the fastq-dump function of the sra-tools.

File is also available on course website. 

```{r}
sraInfo<-read.csv("https://raw.githubusercontent.com/statOmics/SGA/airwaySeqData/SraRunInfo.csv")
sraInfo$SampleName <- as.factor(sraInfo$SampleName)
levels(pdata$SampleName)==levels(sraInfo$SampleName)
```

SampleNames are can be linked.

```{r}
pdata<-merge(pdata,sraInfo,by="SampleName")
pdata$Run
```



We do not have the Albuterol samples

```{r}
rownames(pdata) <- pdata$Run
pdata <- pdata[-grep("Albuterol",pdata[,"treatment:ch1"]),]
pdata[,grep(":ch1",colnames(pdata))]
```

```{r}
pdata$cellLine <- pdata[, grep("cell line:ch1", colnames(pdata))] %>%
  as.factor
pdata$treatment<-pdata[, grep("treatment:ch1", colnames(pdata))] %>%
  as.factor
pdata$treatment<-relevel(pdata$treatment, "Untreated")
```



## Alignment-free transcript quantification

Software such as *kallisto*
[@Bray2016Near], *Salmon* [@Patro2017Salmon] and *Sailfish*
[@Patro2014Sailfish], as well as other transcript quantification methods like
*Cufflinks* [@Trapnell2010Cufflinks; @Trapnell2013Cufflinks2] and *RSEM*
[@Li2011RSEM], differ from the counting methods introduced in the previous tutorials in that they
provide quantifications (usually both as counts and as TPMs) for each
*transcript*. These can then be summarized on the gene level by adding all
values for transcripts from the same gene. A simple way to import results from
these packages into R is provided by the `tximport` and
`tximeta` packages. Here, *tximport* reads the quantifications into a
list of matrices, and *tximeta* aggregates the information into a
*SummarizedExperiment* object, and also automatically adds additional
annotations for the features. Both packages can return quantifications on the
transcript level or aggregate them on the gene level. They also calculate
average transcript lengths for each gene and each sample, which can be used as
offsets to improve the differential expression analysis by accounting for
differential isoform usage across samples [@Soneson2015Differential].
aturecounts object


### Libraries

```{r}
suppressPackageStartupMessages({
  if(!"tximeta" %in% installed.packages()[,1]) BiocManager::install("tximeta")
  if(!"org.Hs.eg.db" %in% installed.packages()[,1]) BiocManager::install("org.Hs.eg.db")
  library(tximeta)
  library(org.Hs.eg.db)
  library(SummarizedExperiment)
  library(ggplot2)
  library(tidyverse)
})
```

### Importing Salmon abundance quantification

The code below imports the *Salmon* quantifications into R using the *tximeta*
package. Note how the transcriptome that was used for the quantification is
automatically recognized and used to annotate the resulting data object. In
order for this to work, *tximeta* requires that the output folder structure from
Salmon is retained, since it reads information from the associated log files in
addition to the quantified abundances themselves. With the `addIds()` function,
we can add additional annotation columns.

```{r}

## List all quant.sf output files from Salmon
salmonfiles <- paste0("SGA-airwaySeqData/salmon/", pdata$Run, "/quant.sf")
names(salmonfiles) <- pdata$Run
stopifnot(all(file.exists(salmonfiles)))

## Add a column "files" to the metadata table. This table must contain at least
## two columns: "names" and "files"
coldata <- cbind(pdata, files = salmonfiles, stringsAsFactors = FALSE)
coldata$names <- coldata$Run

## Import quantifications on the transcript level
st <- tximeta::tximeta(coldata,importer=read.delim)

## Summarize quantifications on the gene level
sg <- tximeta::summarizeToGene(st)

## Add gene symbols
sg <- tximeta::addIds(sg, "SYMBOL", gene = TRUE)
sg
```

Note that *Salmon* returns *estimated* counts, which are not necessarily
integers. They may need to be rounded before they are passed to count-based
statistical methods (e.g. *DESeq2*).
This is not necesary for *edgeR*.


# DE analysis at the gene level

## Setup count object in edgeR

```{r}
suppressPackageStartupMessages({
  library(edgeR)
})
dge <- makeDGEList(sg)
```


## Filtering and normalisation

```{r}
design <- model.matrix(~treatment+cellLine,data=colData(sg))
keep <- filterByExpr(dge, design)
dge <- dge[keep, ,keep.lib.sizes=FALSE]
dge <- calcNormFactors(dge)
dge$samples
```

Note, that the recalculation of libsizes and normalisation factors are overruled because we already have calculated user defined offsets.

## Data exploration

One way to reduce dimensionality is the use of multidimensional scaling (MDS). For MDS, we first have to calculate all pairwise distances between our objects (samples in this case), and then create a (typically) two-dimensional representation where these pre-calculated distances are represented as accurately as possible. This means that depending on how the pairwise sample distances are defined, the two-dimensional plot can be very different, and it is important to choose a distance that is suitable for the type of data at hand.

edgeR contains a function plotMDS, which operates on a DGEList object and generates a two-dimensional MDS representation of the samples. The default distance between two samples can be interpreted as the "typical" log fold change between the two samples, for the genes that are most different between them (by default, the top 500 genes, but this can be modified). We generate an MDS plot from the DGEList object dge, coloring by the treatment and using different plot symbols for different cell lines.

```{r}
colnames(dge) <- paste0(substr(pdata$cellLine,1,3),"_",substr(pdata$treatment,1,3))
plotMDS(dge, top = 500,col=as.double(colData(sg)$cellLine))
```



## Model

### Estimate the overdispersion.

```{r}
dge <- estimateDisp(dge, design)
plotBCV(dge)
```

## Fit the genelevel quasi-NB model.

```{r}
fit <- glmQLFit(dge, design)
ftest <- glmQLFTest(fit, coef = "treatmentDexamethasone")
ttAll <- topTags(ftest, n = nrow(dge)) # all genes
hist(ttAll$table$PValue)
tt <- topTags(ftest, n = nrow(dge), p.value = 0.05) # genes with adj.p<0.05
head(tt)
nrow(tt)
```


## Plots

We first make a volcanoplot and an MA plot.

```{r}
library(ggplot2)
volcano <-  ggplot(ttAll$table,aes(x=logFC,y=-log10(PValue),color=FDR<0.05)) + geom_point() + scale_color_manual(values=c("black","red"))
volcano
plotSmear(ftest, de.tags = row.names(tt$table))
heatmap(cpm(dge,log=TRUE)[rownames(tt$table)[1:30],])
```

# Transcript level analysis 

## Setup count table

```{r}
dte <- DGEList(assays(st)$counts)
design <- model.matrix(~treatment+cellLine,data=colData(st))
keep <- filterByExpr(dte, design)
dte <- dte[keep, ,keep.lib.sizes=FALSE]
dte <- calcNormFactors(dte)
dte$samples
```

## Data Exploration

```{r}
colnames(dte) <- paste0(substr(pdata$cellLine,1,3),"_",substr(pdata$treatment,1,3))
plotMDS(dte, top = 500,col=as.double(colData(st)$cellLine))
```

## Model

### Estimate transcript level dispersion

```{r}
dte <- estimateDisp(dte, design)
plotBCV(dte)
```

Note that we observe problems with respect to the dispersion estimation! 

This is probably due to over excess zero's with respect to the NB distribution.

We have seen similar patterns in the counts of plate based single cell RNA-seq data that are not using UMI's (unique molecular identifiers). 

### Fit the transcript level quasi-NB model.

```{r}
fitt <- glmQLFit(dte, design)
ftestt <- glmQLFTest(fitt, coef = "treatmentDexamethasone")
tttAll <- topTags(ftestt, n = nrow(dge)) # all genes
hist(tttAll$table$PValue)
ttt <- topTags(ftestt, n = nrow(dge), p.value = 0.05) # genes with adj.p<0.05
head(ttt)
nrow(ttt)
```

## Plots

We first make a volcanoplot and an MA plot.

```{r}
library(ggplot2)
volcano <-  ggplot(tttAll$table,aes(x=logFC,y=-log10(PValue),color=FDR<0.05)) + geom_point() + scale_color_manual(values=c("black","red"))
volcano
plotSmear(ftestt, de.tags = row.names(tt$table))
heatmap(cpm(dte,log=TRUE)[rownames(ttt$table)[1:30],])
```