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. Omdat de orginele count table erg groot is gebruikt het script een miniatuur versie gemaakt van 1000 willekeurige features.

Preliminaries

load all required packages

library("DESeq2")
library("gplots")

load data

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

#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.

head(colnames(CT))
## [1] "probe"                "X101_NORM.htseq.txt"  "X102_NORM.htseq.txt" 
## [4] "X104_NORM.htseq.txt"  "X105_NORM_.htseq.txt" "X106_NORM.htseq.txt"
head(db$rnaseq.id)
## [1] X9_0_T06_90304   X613_0_T07_90064 X377_0_T06_90091 X535_0_T07_90052
## [5] X433_0_T06_90135 X25_0_T06_90258 
## 232 Levels: X101_NORM_ X102_NORM_ X104_NORM_ X105_NORM_ ... X99_NORM_

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:

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).

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.

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.

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.

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)
## log2 fold change (MLE): asthma.ics PersA no ICS vs ClinR 
## LRT p-value: '~ gender + asthma.ics' vs '~ gender' 
## DataFrame with 4 rows and 6 columns
##                         baseMean     log2FoldChange              lfcSE
##                        <numeric>          <numeric>          <numeric>
## ENSG00000130779 1789.38912546738 0.0478949207887351 0.0654612878978894
## ENSG00000119535 215.440807373427  0.206590053580329  0.465871144277395
## ENSG00000182957 562.782379494973  0.097333381899491  0.100818494331977
## ENSG00000100604 16.9398710766147   -0.5828701888104  0.239816945028832
##                             stat               pvalue                 padj
##                        <numeric>            <numeric>            <numeric>
## ENSG00000130779 58.7554612212198 5.29625676133409e-12 3.49023320571917e-09
## ENSG00000119535 52.7191486758725 9.75644624431816e-11 3.21474903750283e-08
## ENSG00000182957 48.8988210026882  6.1296215354578e-10 1.34647353062223e-07
## ENSG00000100604 43.5464068704941 7.96970085754099e-09 1.31300821627988e-06
summary(res)
## 
## out of 983 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 75, 7.6%
## LFC < 0 (down)     : 63, 6.4%
## outliers [1]       : 0, 0%
## low counts [2]     : 326, 33%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
sum(res$padj < 0.05, na.rm=T)
## [1] 94
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 ENSG00000130779, maar waarschijnlijk ben je geinteresseerd in een specifiek contrast, bijvoorbeeld asthma versus healthy:

dds.wald <- DESeq(dds, test="Wald", betaPrior = TRUE)
## using pre-existing size factors
## estimating dispersions
## found already estimated dispersions, replacing these
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 33 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
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 ENSG00000119535. Hierna zul je zelf je weg moeten vinden. Meer informatie vind je hier en hier.

session info

Packages used

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
## 
## locale:
##  [1] LC_CTYPE=nl_NL.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=nl_NL.UTF-8        LC_COLLATE=nl_NL.UTF-8    
##  [5] LC_MONETARY=nl_NL.UTF-8    LC_MESSAGES=nl_NL.UTF-8   
##  [7] LC_PAPER=nl_NL.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=nl_NL.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] gplots_3.0.1.1              DESeq2_1.24.0              
##  [3] SummarizedExperiment_1.14.0 DelayedArray_0.10.0        
##  [5] BiocParallel_1.18.0         matrixStats_0.54.0         
##  [7] Biobase_2.44.0              GenomicRanges_1.36.0       
##  [9] GenomeInfoDb_1.20.0         IRanges_2.18.0             
## [11] S4Vectors_0.22.0            BiocGenerics_0.30.0        
## 
## loaded via a namespace (and not attached):
##  [1] bit64_0.9-7            splines_3.6.1          gtools_3.8.1          
##  [4] Formula_1.2-3          assertthat_0.2.1       latticeExtra_0.6-28   
##  [7] blob_1.1.1             GenomeInfoDbData_1.2.1 yaml_2.2.0            
## [10] pillar_1.4.0           RSQLite_2.1.1          backports_1.1.4       
## [13] lattice_0.20-38        glue_1.3.1             digest_0.6.18         
## [16] RColorBrewer_1.1-2     XVector_0.24.0         checkmate_1.9.3       
## [19] colorspace_1.4-1       htmltools_0.3.6        Matrix_1.2-18         
## [22] plyr_1.8.4             XML_3.98-1.19          pkgconfig_2.0.2       
## [25] genefilter_1.66.0      zlibbioc_1.30.0        purrr_0.3.2           
## [28] xtable_1.8-4           scales_1.0.0           gdata_2.18.0          
## [31] htmlTable_1.13.1       tibble_2.1.1           annotate_1.62.0       
## [34] ggplot2_3.1.1          nnet_7.3-12            lazyeval_0.2.2        
## [37] survival_3.1-7         magrittr_1.5           crayon_1.3.4          
## [40] memoise_1.1.0          evaluate_0.13          foreign_0.8-72        
## [43] tools_3.6.1            data.table_1.12.2      stringr_1.4.0         
## [46] locfit_1.5-9.1         munsell_0.5.0          cluster_2.1.0         
## [49] AnnotationDbi_1.46.0   compiler_3.6.1         caTools_1.17.1.2      
## [52] rlang_0.3.4            grid_3.6.1             RCurl_1.95-4.12       
## [55] rstudioapi_0.10        htmlwidgets_1.3        labeling_0.3          
## [58] bitops_1.0-6           base64enc_0.1-3        rmarkdown_1.12        
## [61] gtable_0.3.0           DBI_1.0.0              R6_2.4.0              
## [64] gridExtra_2.3          knitr_1.22             dplyr_0.8.1           
## [67] bit_1.1-14             Hmisc_4.2-0            KernSmooth_2.23-16    
## [70] stringi_1.4.3          Rcpp_1.0.1             geneplotter_1.62.0    
## [73] rpart_4.1-15           acepack_1.4.1          tidyselect_0.2.5      
## [76] xfun_0.7