Hello everyone, I'm trying to figure out whether the DESeq package and its GLM function will work for a complex mRNA-seq time course. I have a time course of five time points, with treated and untreated samples at each time point, and multiple replicates. The preliminary results from this program make me sure I haven't set up the contrast correctly, because the results look wrong to me (more than half of the genome is differentially expressed at an FDR <1%). Am I missing something? I note that in the DESeq package vignette, the example case involves binary conditions (treated/untreated vs paired/single end). Will DESeq even work for a time course of five time points, or am I wasting my time? My script, with output, is below. I apologize in advance - it is quite lengthy! any help on this would be so fantastic! Elena > options(digits=3) > setwd("C:\\Rdata\\DESeq\\GLMs") > library(DESeq) > > # input time course data > countsTable<-read.delim("input2_mut_ALL.txt",header=TRUE, stringsAsFactors=TRUE) > head(countsTable) gene untreat.1 untreat.2 D.1h.1 D.1h.2 U.1h.1 U1h.2 D.2h.1 D.2h.2 U.2h.1 U.2h.2 D.6h.1 D.6h.2 U.6h.1 U.6h.2 D.18h.1 D.18h.2 U.18h.1 U.18h.2 1 2L52.1 98 55 100 103 98 69 70 64 73 36 59 49 33 32 92 83 47 47 2 2RSSE.1 423 261 374 459 327 380 216 215 191 223 146 155 117 110 196 133 110 110 3 2RSSE.2 8 9 8 16 4 8 4 3 2 3 2 1 2 0 2 1 1 1 4 2RSSE.3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 3R5.1 74 38 96 46 69 52 56 41 38 32 31 14 28 31 38 33 81 81 6 3R5.2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > rownames(countsTable) <- countsTable$gene > countsTable <- countsTable [ , -1] > > # make data.frame > design <- read.table("C:\\Rdata\\DESeq\\puf8lip1\\GLMs\\design2.txt", header = T, row.names=1) > design<- as.data.frame(design) > is.data.frame(design) [1] TRUE > > design treat.type time.hr untreat.1 none 0 untreat.2 none 0 D.1h.1 vehicle 1 D.1h.2 vehicle 1 U.1h.1 drug 1 U.1h.2 drug 1 D.2h.2 vehicle 2 D.2h.1 vehicle 2 U.2h.1 drug 2 U.2h.2 drug 2 D.6h.1 vehicle 6 D.6h.2 vehicle 6 U.6h.1 drug 6 U.6h.2 drug 6 D.18h.1 vehicle 18 D.18h.2 vehicle 18 U.18h.1 drug 18 U.18h.2 drug 18 > > # create a count data object > cds <- newCountDataSet( countsTable, design) > head(counts(cds)) untreat.1 untreat.2 D.1h.1 D.1h.2 U.1h.1 U1h.2 D.2h.1 D.2h.2 U.2h.1 U.2h.2 D.6h.1 D.6h.2 U.6h.1 U.6h.2 D.18h.1 D.18h.2 U.18h.1 U.18h.2 2L52.1 98 55 100 103 98 69 70 64 73 36 59 49 33 32 92 83 47 47 2RSSE.1 423 261 374 459 327 380 216 215 191 223 146 155 117 110 196 133 110 110 2RSSE.2 8 9 8 16 4 8 4 3 2 3 2 1 2 0 2 1 1 1 2RSSE.3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3R5.1 74 38 96 46 69 52 56 41 38 32 31 14 28 31 38 33 81 81 3R5.2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > > # estimate library size of count data set > cds <- estimateSizeFactors ( cds ) > sizeFactors( cds) untreat.1 untreat.2 D.1h.1 D.1h.2 U.1h.1 U1h.2 D.2h.1 D.2h.2 U.2h.1 U.2h.2 D.6h.1 D.6h.2 U.6h.1 U.6h.2 D.18h.1 D.18h.2 2.307 1.149 1.435 1.301 1.370 1.230 1.152 1.040 1.079 0.926 0.602 0.593 0.620 0.596 0.926 0.694 U.18h.1 U.18h.2 0.887 0.887 > > # estimate dispersion (biological variation) > # this is done for each condition/factor > cds <- estimateDispersions( cds , method = "pooled") > > ### Check the fit ###### > > plotDispEsts <- function( cds ) + { + plot( + rowMeans( counts( cds, normalized=TRUE ) ), + fitInfo(cds)$perGeneDispEsts, + pch = 8, log="xy",xlab = "Per Gene Average Expression Level (in counts)", ylab= "Per Gene Variance Estimate") + xg <- 10^seq( -.5, 5, length.out=300 ) + lines( xg, fitInfo(cds)$dispFun( xg ), col="red", pch=16 ) + } > > plotDispEsts(cds) Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log) : 6941 x values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 6968 y values <= 0 omitted from logarithmic plot > head(fData(cds)) #dispersion values are stored in the feature data slot of cds disp_pooled 2L52.1 0.0245 2RSSE.1 0.0192 2RSSE.2 0.3018 2RSSE.3 Inf 3R5.1 0.0308 3R5.2 Inf > str(fitInfo(cds)) #verify that the displ_pooled List of 5 $ perGeneDispEsts: num [1:27924] 0.00641 0.0192 0.10196 NaN 0.02612 ... $ dispFunc :function (q) ..- attr(*, "coefficients")= Named num [1:2] 0.00894 1.04347 .. ..- attr(*, "names")= chr [1:2] "asymptDisp" "extraPois" ..- attr(*, "fitType")= chr "parametric" $ fittedDispEsts : num [1:27924] 0.0245 0.0137 0.3018 Inf 0.0308 ... $ df : int 9 $ sharingMode : chr "maximum" > > ###### Fit data to Gen. Lin. Model ##### > > fit1 <- fitNbinomGLMs(cds, count ~ treat.type + time.hr) .................. There were 50 or more warnings (use warnings() to see the first 50) > fit0 <- fitNbinomGLMs(cds, count ~ treat.type) .................. Warning messages: 1: glm.fit: algorithm did not converge 2: glm.fit: algorithm did not converge 3: glm.fit: algorithm did not converge 4: glm.fit: algorithm did not converge 5: glm.fit: algorithm did not converge > > ###DIFF EXPRESSION ########## > > str(fit1) 'data.frame': 27924 obs. of 6 variables: $ (Intercept) : num 6.24 8.16 2.67 NA 5.28 ... $ treat.typenone : num -0.757 -0.484 -0.225 NA -0.255 ... $ treat.typeU0126: num -0.532 -0.295 -0.713 NA 0.226 ... $ time.hr : num 0.0171 -0.0392 -0.1143 NA 0.035 ... $ deviance : num 13.95 22.92 8.29 NA 23.29 ... $ converged : logi TRUE TRUE TRUE NA TRUE NA ... - attr(*, "df.residual")= Named num 14 ..- attr(*, "names")= chr "2L52.1" > pvalsGLM <- nbinomGLMTest (fit1, fit0) > padjGLM <- p.adjust (pvalsGLM, method = "BH") # correct for mult.test > head(fit1) # hope to see fully-converged, ie TRUE (Intercept) treat.typenone treat.typeU0126 time.hr deviance converged 2L52.1 6.24 -0.757 -0.532 0.0171 13.95 TRUE 2RSSE.1 8.16 -0.484 -0.295 -0.0392 22.92 TRUE 2RSSE.2 2.67 -0.225 -0.713 -0.1143 8.29 TRUE 2RSSE.3 NA NA NA NA NA NA 3R5.1 5.28 -0.255 0.226 0.0350 23.29 TRUE 3R5.2 NA NA NA NA NA NA > head(padjGLM) [1] 1.73e-01 1.39e-05 3.66e-02 NA 6.29e-03 NA > > hist(pvalsGLM, breaks=100, col="skyblue", border="slateblue", main="") > > > ######### Filter for significant genes at p < 0.05 ######## > res <- cbind(fit1,pvalsGLM,padjGLM) > resSig <- res[ res$padjGLM < 0.05, ] > > ######### LAST PART: WRITE RESULTS TO .CSV FILE ########### > write.table(resSig, file="DESeqOutput.txt", quote=FALSE) > > ############ > sessionInfo() R version 2.14.0 (2011-10-31) Platform: i386-pc-mingw32/i386 (32-bit) locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C [5] LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] pasilla_0.2.10 DESeq_1.6.1 locfit_1.5-6 lattice_0.20-0 akima_0.5-4 Biobase_2.14.0 loaded via a namespace (and not attached): [1] annotate_1.32.0 AnnotationDbi_1.16.10 DBI_0.2-5 genefilter_1.36.0 geneplotter_1.32.1 grid_2.14.0 IRanges_1.12.5 [8] RColorBrewer_1.0-5 RSQLite_0.11.1 splines_2.14.0 survival_2.36-10 tools_2.14.0 xtable_1.6-0 >

Dear Elena thank you for the detailed report. This is perhaps not the answer you are expecting, but perhaps would it be better to first figure out _what_ you want to achieve before discussing the _how_? When in a dataset like yours you test against the null hypothesis of 'no effect of time', it is pretty plausible that the tests reject very many genes. For such time courses, I have often found it more instructive to cluster the response patterns (e.g. using k-means on the profiles themselves or on the differences between treatment and control) and then, if needed, to set up statistical tests for a specific hypothesis induced from that data exploration. These tools are useful in that context: - variance stabilising transformation of the count data (see DESeq vignette) to bring them on a scale where the clustering metric (e.g. Euclidean) is more appropriate, - variance filtering, to first remove genes that do not chance appreciably at all and thus to simplify the clustering, - heatmaps. Hope this helps Wolfgang Jan/5/12 12:33 AM, Elena Sorokin scripsit:: > Hello everyone, > > I'm trying to figure out whether the DESeq package and its GLM function > will work for a complex mRNA-seq time course. I have a time course of > five time points, with treated and untreated samples at each time point, > and multiple replicates. The preliminary results from this program make > me sure I haven't set up the contrast correctly, because the results > look wrong to me (more than half of the genome is differentially > expressed at an FDR <1%). Am I missing something? I note that in the > DESeq package vignette, the example case involves binary conditions > (treated/untreated vs paired/single end). Will DESeq even work for a > time course of five time points, or am I wasting my time? My script, > with output, is below. I apologize in advance - it is quite lengthy! > > any help on this would be so fantastic! > Elena > > > options(digits=3) > > setwd("C:\\Rdata\\DESeq\\GLMs") > > library(DESeq) > > > > # input time course data > > countsTable<-read.delim("input2_mut_ALL.txt",header=TRUE, > stringsAsFactors=TRUE) > > head(countsTable) > gene untreat.1 untreat.2 D.1h.1 D.1h.2 U.1h.1 U1h.2 D.2h.1 D.2h.2 U.2h.1 > U.2h.2 D.6h.1 D.6h.2 U.6h.1 U.6h.2 D.18h.1 D.18h.2 U.18h.1 U.18h.2 > 1 2L52.1 98 55 100 103 98 69 70 64 73 36 59 49 33 32 92 83 47 47 > 2 2RSSE.1 423 261 374 459 327 380 216 215 191 223 146 155 117 110 196 > 133 110 110 > 3 2RSSE.2 8 9 8 16 4 8 4 3 2 3 2 1 2 0 2 1 1 1 > 4 2RSSE.3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > 5 3R5.1 74 38 96 46 69 52 56 41 38 32 31 14 28 31 38 33 81 81 > 6 3R5.2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > > > > rownames(countsTable) <- countsTable$gene > > countsTable <- countsTable [ , -1] > > > > # make data.frame > > design <- read.table("C:\\Rdata\\DESeq\\puf8lip1\\GLMs\\design2.txt", > header = T, row.names=1) > > design<- as.data.frame(design) > > is.data.frame(design) > [1] TRUE > > > > design > treat.type time.hr > untreat.1 none 0 > untreat.2 none 0 > D.1h.1 vehicle 1 > D.1h.2 vehicle 1 > U.1h.1 drug 1 > U.1h.2 drug 1 > D.2h.2 vehicle 2 > D.2h.1 vehicle 2 > U.2h.1 drug 2 > U.2h.2 drug 2 > D.6h.1 vehicle 6 > D.6h.2 vehicle 6 > U.6h.1 drug 6 > U.6h.2 drug 6 > D.18h.1 vehicle 18 > D.18h.2 vehicle 18 > U.18h.1 drug 18 > U.18h.2 drug 18 > > > > # create a count data object > > cds <- newCountDataSet( countsTable, design) > > head(counts(cds)) > untreat.1 untreat.2 D.1h.1 D.1h.2 U.1h.1 U1h.2 D.2h.1 D.2h.2 U.2h.1 > U.2h.2 D.6h.1 D.6h.2 U.6h.1 U.6h.2 D.18h.1 D.18h.2 U.18h.1 U.18h.2 > 2L52.1 98 55 100 103 98 69 70 64 73 36 59 49 33 32 92 83 47 47 > 2RSSE.1 423 261 374 459 327 380 216 215 191 223 146 155 117 110 196 133 > 110 110 > 2RSSE.2 8 9 8 16 4 8 4 3 2 3 2 1 2 0 2 1 1 1 > 2RSSE.3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > 3R5.1 74 38 96 46 69 52 56 41 38 32 31 14 28 31 38 33 81 81 > 3R5.2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > > > > # estimate library size of count data set > > cds <- estimateSizeFactors ( cds ) > > sizeFactors( cds) > untreat.1 untreat.2 D.1h.1 D.1h.2 U.1h.1 U1h.2 D.2h.1 D.2h.2 U.2h.1 > U.2h.2 D.6h.1 D.6h.2 U.6h.1 U.6h.2 D.18h.1 D.18h.2 > 2.307 1.149 1.435 1.301 1.370 1.230 1.152 1.040 1.079 0.926 0.602 0.593 > 0.620 0.596 0.926 0.694 > U.18h.1 U.18h.2 > 0.887 0.887 > > > > # estimate dispersion (biological variation) > > # this is done for each condition/factor > > cds <- estimateDispersions( cds , method = "pooled") > > > > ### Check the fit ###### > > > > plotDispEsts <- function( cds ) > + { > + plot( > + rowMeans( counts( cds, normalized=TRUE ) ), > + fitInfo(cds)$perGeneDispEsts, > + pch = 8, log="xy",xlab = "Per Gene Average Expression Level (in > counts)", ylab= "Per Gene Variance Estimate") > + xg <- 10^seq( -.5, 5, length.out=300 ) > + lines( xg, fitInfo(cds)$dispFun( xg ), col="red", pch=16 ) > + } > > > > plotDispEsts(cds) > Warning messages: > 1: In xy.coords(x, y, xlabel, ylabel, log) : > 6941 x values <= 0 omitted from logarithmic plot > 2: In xy.coords(x, y, xlabel, ylabel, log) : > 6968 y values <= 0 omitted from logarithmic plot > > head(fData(cds)) #dispersion values are stored in the feature data > slot of cds > disp_pooled > 2L52.1 0.0245 > 2RSSE.1 0.0192 > 2RSSE.2 0.3018 > 2RSSE.3 Inf > 3R5.1 0.0308 > 3R5.2 Inf > > str(fitInfo(cds)) #verify that the displ_pooled > List of 5 > $ perGeneDispEsts: num [1:27924] 0.00641 0.0192 0.10196 NaN 0.02612 ... > $ dispFunc :function (q) > ..- attr(*, "coefficients")= Named num [1:2] 0.00894 1.04347 > .. ..- attr(*, "names")= chr [1:2] "asymptDisp" "extraPois" > ..- attr(*, "fitType")= chr "parametric" > $ fittedDispEsts : num [1:27924] 0.0245 0.0137 0.3018 Inf 0.0308 ... > $ df : int 9 > $ sharingMode : chr "maximum" > > > > ###### Fit data to Gen. Lin. Model ##### > > > > fit1 <- fitNbinomGLMs(cds, count ~ treat.type + time.hr) > .................. > There were 50 or more warnings (use warnings() to see the first 50) > > fit0 <- fitNbinomGLMs(cds, count ~ treat.type) > .................. > Warning messages: > 1: glm.fit: algorithm did not converge > 2: glm.fit: algorithm did not converge > 3: glm.fit: algorithm did not converge > 4: glm.fit: algorithm did not converge > 5: glm.fit: algorithm did not converge > > > > ###DIFF EXPRESSION ########## > > > > str(fit1) > 'data.frame': 27924 obs. of 6 variables: > $ (Intercept) : num 6.24 8.16 2.67 NA 5.28 ... > $ treat.typenone : num -0.757 -0.484 -0.225 NA -0.255 ... > $ treat.typeU0126: num -0.532 -0.295 -0.713 NA 0.226 ... > $ time.hr : num 0.0171 -0.0392 -0.1143 NA 0.035 ... > $ deviance : num 13.95 22.92 8.29 NA 23.29 ... > $ converged : logi TRUE TRUE TRUE NA TRUE NA ... > - attr(*, "df.residual")= Named num 14 > ..- attr(*, "names")= chr "2L52.1" > > pvalsGLM <- nbinomGLMTest (fit1, fit0) > > padjGLM <- p.adjust (pvalsGLM, method = "BH") # correct for mult.test > > head(fit1) # hope to see fully-converged, ie TRUE > (Intercept) treat.typenone treat.typeU0126 time.hr deviance converged > 2L52.1 6.24 -0.757 -0.532 0.0171 13.95 TRUE > 2RSSE.1 8.16 -0.484 -0.295 -0.0392 22.92 TRUE > 2RSSE.2 2.67 -0.225 -0.713 -0.1143 8.29 TRUE > 2RSSE.3 NA NA NA NA NA NA > 3R5.1 5.28 -0.255 0.226 0.0350 23.29 TRUE > 3R5.2 NA NA NA NA NA NA > > head(padjGLM) > [1] 1.73e-01 1.39e-05 3.66e-02 NA 6.29e-03 NA > > > > hist(pvalsGLM, breaks=100, col="skyblue", border="slateblue", main="") > > > > > > ######### Filter for significant genes at p < 0.05 ######## > > res <- cbind(fit1,pvalsGLM,padjGLM) > > resSig <- res[ res$padjGLM < 0.05, ] > > > > ######### LAST PART: WRITE RESULTS TO .CSV FILE ########### > > write.table(resSig, file="DESeqOutput.txt", quote=FALSE) > > > > ############ > > sessionInfo() > R version 2.14.0 (2011-10-31) > Platform: i386-pc-mingw32/i386 (32-bit) > > locale: > [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United > States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C > [5] LC_TIME=English_United States.1252 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] pasilla_0.2.10 DESeq_1.6.1 locfit_1.5-6 lattice_0.20-0 akima_0.5-4 > Biobase_2.14.0 > > loaded via a namespace (and not attached): > [1] annotate_1.32.0 AnnotationDbi_1.16.10 DBI_0.2-5 genefilter_1.36.0 > geneplotter_1.32.1 grid_2.14.0 IRanges_1.12.5 > [8] RColorBrewer_1.0-5 RSQLite_0.11.1 splines_2.14.0 survival_2.36-10 > tools_2.14.0 xtable_1.6-0 > > > > _______________________________________________ > Bioconductor mailing list > Bioconductor at r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor -- Best wishes Wolfgang Wolfgang Huber EMBL http://www.embl.de/research/units/genome_biology/huber

Hi Elena With this here fit1 <- fitNbinomGLMs(cds, count ~ treat.type + time.hr) fit0 <- fitNbinomGLMs(cds, count ~ treat.type) you are testing against the null hypothesis that a gene's expression is independent of the time of the sample, i.e., you will get a hit whenever you can reject the claim that the gene has the same expression in all five time points. This is a pretty narrow null hypothesis. It seems quite plausible to me that around 50% of the genes change noticeable over time. Note also that treatment is a factor in both your full model and in your null model, so you are not testing whether the treatment has an effect. I suspect that this might not be what you want. In a simple comparison of treatment vs control the question one wants to test is usually fairly obvious: does the treatment influence the expression or not. For time course data one usually does not have such an obvious yes-or-no question, and hence, you may need to tell us more about what exactly you want to test before we can advise you. Simon

Dear Wolfgang and Simon, Thanks very much for your replies and help! OK, let me try to clarify: what I'm attempting is typical for a time course where there's a control and an experimental group. I'm attempting to find patterns of expression that are specifically associated with the experimental group, not the control. I have paired treatments of control treatment and drug treatment. What I'd like to end up with is a single FDR or p-value for each row/gene of my heatmap - a statistic indicating which patterns are significant, and which are due to chance. My question is, can DESeq perform an ANOVA-like test with five means? If so, how would you recommend I set up the test? My previous attempt was clearly flawed... Many thanks for your time and help with this! Best wishes, Elena Message: 18 Date: Thu, 05 Jan 2012 09:38:21 +0100 From: Wolfgang Huber<whuber@embl.de> To:bioconductor@r-project.org Subject: Re: [BioC] DESeq for an mRNA-seq time course Message-ID:<4F05617D.40507@embl.de> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Dear Elena thank you for the detailed report. This is perhaps not the answer you are expecting, but perhaps would it be better to first figure out _what_ you want to achieve before discussing the _how_? When in a dataset like yours you test against the null hypothesis of 'no effect of time', it is pretty plausible that the tests reject very many genes. For such time courses, I have often found it more instructive to cluster the response patterns (e.g. using k-means on the profiles themselves or on the differences between treatment and control) and then, if needed, to set up statistical tests for a specific hypothesis induced from that data exploration. These tools are useful in that context: - variance stabilising transformation of the count data (see DESeq vignette) to bring them on a scale where the clustering metric (e.g. Euclidean) is more appropriate, - variance filtering, to first remove genes that do not chance appreciably at all and thus to simplify the clustering, - heatmaps. Hope this helps Wolfgang ------------------------------ Message: 19 Date: Thu, 05 Jan 2012 10:25:57 +0100 From: Simon Anders<anders@embl.de> To:bioconductor@r-project.org Subject: Re: [BioC] DESeq for an mRNA-seq time course Message-ID:<4F056CA5.7050004@embl.de> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Hi Elena With this here fit1<- fitNbinomGLMs(cds, count ~ treat.type + time.hr) fit0<- fitNbinomGLMs(cds, count ~ treat.type) you are testing against the null hypothesis that a gene's expression is independent of the time of the sample, i.e., you will get a hit whenever you can reject the claim that the gene has the same expression in all five time points. This is a pretty narrow null hypothesis. It seems quite plausible to me that around 50% of the genes change noticeable over time. Note also that treatment is a factor in both your full model and in your null model, so you are not testing whether the treatment has an effect. I suspect that this might not be what you want. In a simple comparison of treatment vs control the question one wants to test is usually fairly obvious: does the treatment influence the expression or not. For time course data one usually does not have such an obvious yes-or-no question, and hence, you may need to tell us more about what exactly you want to test before we can advise you. Simon [[alternative HTML version deleted]]