Question

Deseq2 result for a variable

0

Entering edit mode

yoursbassanio ▴ 10

@yoursbassanio-12717

Last seen 4.3 years ago

Hi,

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

resultsNames(dds)
 [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.?

sessionInfo()

R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
[3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
[7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
[9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=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

deseq2 differential expression results • 1.6k views

ADD COMMENT • link updated 7.2 years ago by Michael Love 43k • written 7.2 years ago by yoursbassanio ▴ 10

0

Entering edit mode

Hi ,

Is output from

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

and

ResultSexF <- lfcShrink(dds, coef=4)

from above resultsNames(dds) same?

ADD REPLY • link updated 7.2 years ago by Michael Love 43k • written 7.2 years ago by yoursbassanio ▴ 10

0

Entering edit mode

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 REPLY • link 7.2 years ago Michael Love 43k

score 1 · Answer 1 · 2017-10-15

1

Entering edit mode

Michael Love 43k

@mikelove

Last seen 5 days ago

United States

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 COMMENT • link 7.2 years ago Michael Love 43k

0

Entering edit mode

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

-0.08279180

0.32866904

0.19715990

0.02425360

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 REPLY • link 7.2 years ago yoursbassanio ▴ 10

1

Entering edit mode

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 REPLY • link 7.2 years ago Michael Love 43k

0

Entering edit mode

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 REPLY • link 7.2 years ago Michael Love 43k

0

Entering edit mode

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 REPLY • link 7.2 years ago yoursbassanio ▴ 10

0

Entering edit mode

See here, the explanation of the log2 fold change:

http://www.bioconductor.org/help/workflows/rnaseqGene/#building-the-results-table

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 REPLY • link 7.2 years ago Michael Love 43k