system_genetics/Tessa_sample_workflow/DESeq2_example.Rmd

---
title: "Sample workflow DESeq2 for Tessa"
author: "Corneel Vermeulen"
date: "`r Sys.Date()`"
output: html_document
---

### GSK asthma remission project

Hoi Tesssa, Dit is een aangepaste versie van het tutorial script dat ik voor Ilse maakte, vanwege een [verandering in DESeq2 VERSION 1.16.0](# https://support.bioconductor.org/p/95695/). Omdat de orginele count table erg groot is gebruikt het script een miniatuur versie gemaakt van 1000 willekeurige features.    
   
### Preliminaries

load all required packages

```{r load_packgs, warning=FALSE, message=FALSE}
library("DESeq2")
library("gplots")
```

### load data

Hier laad je de count table and de database respectievelijk. Je kunt hier evt. je werkmap instellen.

```{r setwd_dataload}
#setwd("/path/to/your/wd")
CT<-read.delim("sample.htseq.txt.table", header = TRUE, sep = "\t", 
               quote = "", fill = TRUE, comment.char = "")
db <-read.csv("db.csv",header=T,sep=";",dec=",")
```

### Prepare data 

Dit is het lastigste stuk. De sample_IDs zijn gemangeld tijdens het verwerken van de samples en dat maakt het matchen aan de klinische data een tijdrovende klus. 

```{r showids}
head(colnames(CT))
head(db$rnaseq.id)
```

Ik heb al wat werk verzet door in de db overal een "X" voor te plakken en alle streepjes te vervangen door underscores. Nu moeten we nog de ".htseq.txt" kwijtraken. Zie ook dat sommige samples "NORM" en anderen "NORM_" genoemd zijn. Verder zijn er meerdere controle samples. Veel lelijke code hier, en waarschijnlijk niet bruikbaar voor jou:

```{r prepdata}
colnames(CT)[1] <- ""
cn.CT <- colnames(CT)
cn.CT <- cn.CT[2:length(cn.CT)] #drop the empty first cell
cn.CT <- gsub(".htseq.txt", "", cn.CT) # remove the suffix
cn.CT <- gsub("NORM.{0,}", "NORM_", cn.CT) # add one trailing underscore
cn.all <- as.character(db$rnaseq.id) 
idx <- match(cn.CT,cn.all) #match the clinical data to the correct samples
idx[which(is.na(idx))] <- which(cn.all == "X263_110_T04_90226") # multiple ref samples
coldata <- db[idx,]
rownames(coldata) <- cn.CT
countdata <- CT[,c(2:ncol(CT))]
rownames(countdata) <- CT[,1]
coldata$batch.extract[which(coldata$rnaseq.id == "X263_110_T04_90226")] <- seq(1:6)
coldata$batch.extract<-as.factor(coldata$batch.extract) #batch nrs are factor levels
```

### QC-filter

We gebruiken alleen samples die de QC doorstaan en de beste reference sample (#2).

```{r qcfilter, message=FALSE}
qc.pass <- scan("rnaseqids_qc4_pass_plusref2.txt", what = "character")
idx.qc.pass <- which(coldata$rnaseq.id %in% qc.pass)
idx.ref2 <- which(coldata$rnaseq.id == "X263_110_T04_90226" & coldata$batch.extract == 2)
idx.qc.pass <- c(idx.qc.pass,idx.ref2)
coldata <- coldata[idx.qc.pass,]
countdata <- countdata[,idx.qc.pass]
colnames(countdata) <- NULL # fixes error
```
### Exploratory analysis

Nu we onze data klaar hebben, wordt alles veel eenvoudiger: We converteren onze count table naar een DESeqDataSet. We geven nog geen model op.

```{r convert}
dds <- DESeqDataSetFromMatrix(countData = countdata,
colData = coldata,
design = ~ 1)
```

Dan kijken we in een PCA of onze dataset zich gedraagt. Het is noodzakelijk de ruwe count data te transformeren daarvoor. Omdat we teveel samples hebben om een rlog transformatie te doen gebruiken we een [trucje](https://support.bioconductor.org/p/77122/).

```{r pca, message = FALSE}
dds <- estimateSizeFactors(dds)
baseMean <- rowMeans(counts(dds, normalized=TRUE))
idx.sub <- sample(which(baseMean > 5), 100)
dds.sub <- dds[idx.sub, ]
dds.sub <- estimateDispersions(dds.sub)
dispersionFunction(dds) <- dispersionFunction(dds.sub)
vsd <- varianceStabilizingTransformation(dds, blind=FALSE)
plotPCA(vsd, intgroup = "gender")
```

### Differential expression

Tenslotte het echte werk. We voeren ons model in en creeren een nieuwe DESeqDataSet (de oude wordt daarbij overschreven). Dan testen we of er een overall verschil is tussen de levels van asthma. Vervolgens kijken we naar een specifiek contrast.

```{r difexp, message = FALSE}
dds <- DESeqDataSetFromMatrix(countData = countdata,
colData = coldata,
design = ~ gender + asthma.ics)
# Overall test (takes long time to run)
dds <- DESeq(dds, test="LRT", full=~ gender + asthma.ics, reduced=~ gender)
res <- results(dds)
res$symbol <- mcols(dds)$symbol
head(res[order(res$pvalue),],4)
summary(res)
sum(res$padj < 0.05, na.rm=T)
topGene <- rownames(res)[which.min(res$padj)]
plotCounts(dds, gene=topGene, intgroup=c("asthma.ics"))
```

Het feature met de meeste variatie tussen de groepen is `r topGene`, maar waarschijnlijk ben je geinteresseerd in een specifiek contrast, bijvoorbeeld asthma versus healthy:

```{r contrast}
dds.wald <- DESeq(dds, test="Wald", betaPrior = TRUE)
res.c <- results(dds.wald, contrast = c("asthma.ics","PersA_no_ICS","H"))
topGene.c <- rownames(res.c)[which.min(res.c$padj)]
```

Het feature dat het meest verschillend tot expressie komt in dit contrast is `r topGene.c`. Hierna zul je zelf je weg moeten vinden. Meer informatie vind je [hier](http://www.bioconductor.org/help/workflows/rnaseqGene/) en [hier](https://www.bioconductor.org/help/course-materials/2015/LearnBioconductorFeb2015/B02.1.1_RNASeqLab.html).

## session info

Packages used

```{r sessionInfo, include=TRUE, echo=TRUE, results='markup'}
sessionInfo()
```