132 lines
5.1 KiB
Plaintext
132 lines
5.1 KiB
Plaintext
---
|
|
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()
|
|
```
|