Question: Deseq2 result for a variable
gravatar for yoursbassanio
4 months ago by
yoursbassanio0 wrote:


The following is my deseq2 design and my commands

dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,directory = directory,design= ~Sex+V11+V12+V13+V14+V15+V16+V17+Visit,ignoreRank=FALSE)
dds <- DESeq(dds)
rld <- rlog(dds)
vsd <- varianceStabilizingTransformation(dds)

With the following commands I was able to get Deseq2 analysis results for contrast. 

ddssex <- results(dds,contrast=c("Sex","M","F"))write.table(ddSSsex, "ddSSsex.xls", sep="\t")

 [1] "Intercept"             "SexF"                  "SexM"                 
 [4] "V11"                   "V12"                   "V13"           
 [7] "V14"                   "V15"                   "V16"                  
[10] "V17"                   "VisitV1"  

I would like to have p-values for variable V11(coef 4). It seems lfcShrink() is not available with my version. Is there a way around for the same.?


R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
[3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
[7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
[9] LC_ADDRESS=C               LC_TELEPHONE=C            

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets
[8] methods   base     

other attached packages:
[1] pheatmap_1.0.7             DESeq2_1.12.4             
[3] SummarizedExperiment_1.4.0 Biobase_2.32.0            
[5] GenomicRanges_1.26.1       GenomeInfoDb_1.8.7        
[7] IRanges_2.8.0              S4Vectors_0.12.0          
[9] BiocGenerics_0.20.0       

loaded via a namespace (and not attached):
[1] Rcpp_0.12.8          RColorBrewer_1.1-2   plyr_1.8.4          
[4] XVector_0.12.1       tools_3.3.1          zlibbioc_1.18.0     
[7] rpart_4.1-10         RSQLite_1.0.0        annotate_1.50.0     
[10] tibble_1.2           gtable_0.2.0         lattice_0.20-34     
[13] Matrix_1.2-6         DBI_0.4-1            gridExtra_2.2.1     
[16] genefilter_1.56.0    cluster_2.0.5        locfit_1.5-9.1      
[19] grid_3.3.1           nnet_7.3-12          data.table_1.10.0   
[22] AnnotationDbi_1.36.0 XML_3.98-1.4         survival_2.39-4     
[25] BiocParallel_1.6.6   foreign_0.8-67       latticeExtra_0.6-28
[28] Formula_1.2-1        geneplotter_1.50.0   ggplot2_2.2.0       
[31] Hmisc_3.17-4         scales_0.4.1         splines_3.3.1       
[34] assertthat_0.1       colorspace_1.3-1     xtable_1.8-2        
[37] acepack_1.3-3.3      lazyeval_0.2.0       munsell_0.4.3
ADD COMMENTlink modified 4 months ago by Michael Love16k • written 4 months ago by yoursbassanio0

Hi ,

Is output from  

Res <- results(dds, name="V11")


ResultSexF <- lfcShrink(dds, coef=4)

from above resultsNames(dds) same?

ADD REPLYlink modified 4 months ago by Michael Love16k • written 4 months ago by yoursbassanio0

This is true, as you can tell from reading the help page for ?results. Or you can also confirm by inspecting the results tables.

ADD REPLYlink written 4 months ago by Michael Love16k
gravatar for Michael Love
4 months ago by
Michael Love16k
United States
Michael Love16k wrote:

Because you have version < 1.16, there is no function lfcShrink() in your version of DESeq2. You should just use results() with ‘contrast’. This will give shrunken LFC estimates. Or you can upgrade to latest version of Bioc and use the new function.

ADD COMMENTlink written 4 months ago by Michael Love16k

Thanks Mike

If I upgrade to latest Deseq2. should I need to rerun from the beginning(It took me 4 days to run Deseq2)?.  Doest contrast needs 3 argument. In my case V11 is continous data example





One more doubt: Does the below design accounts(corrects) for all other variables except Visit(variable of interest)

dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,directory = directory,design= ~Sex+V11+V12+V13+V14+V15+V16+V17+Visit,ignoreRank=FALSE)

ADD REPLYlink written 4 months ago by yoursbassanio0

Depending on what results table you build, having variables ~x + y + z, we would typically say we controlled for the other variables in the GLM when testing one, eg ‘z’ here.

ADD REPLYlink written 4 months ago by Michael Love16k

See ?results, for a single coefficient you should use ‘name’. 

How many samples do you have? That sounds like something isn’t right. 

With many samples, you can prefilter low count genes to speed things up, eg

dds <- estimateSizeFactors(dds)

keep <- rowSums(counts(dds, normalized=TRUE) >= 10) >= x

dds <- dds[keep,]

You can also use parallel=TRUE after you register multiple cores.

ADD REPLYlink modified 4 months ago • written 4 months ago by Michael Love16k

Hi Mike,

I think I found why it was taking this much time.

Firstly because of rlog transformation (I have 300+ samples) and I have changed the rlog to vst. Secondly I was not using parallel(multicore). I fixed  both of them. I have filtered the gene counts outside the deseq2.

One more doubt

AB <- results(dds,contrast=c("Condition","A","B"))

The LFC says direction (up/down) regulation of GeneN in ConditionB with respect to ConditonA? or is it other way?



ADD REPLYlink written 4 months ago by yoursbassanio0

See here, the explanation of the log2 fold change:

Also you can look up the help page for ?results, and the 'contrast' argument.

Additionally, it prints the comparison at the top of the results table when you print it out:

> res <- results(...)
> res
ADD REPLYlink written 4 months ago by Michael Love16k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 264 users visited in the last hour