DESeq2 different stats
2
0
Entering edit mode
@catalina-aguilar-hurtado-6554
Last seen 4.1 years ago
United States

Hi Michael and Simon,
 

I wanted to ask if there was a resent update that could be changing my pvalues in a drastic way?

I ran my analysis on the 18 of August and was getting  2885 DEGs (padj <0.01), today the same .R script (using my same input files)  I get 1559 DEGs (padj <0.01). The output stat, pvalues and padj are very different. I also checked with other datasets that I have and got the same different results.

ddsLRT <- DESeq (dds, test ="LRT", reduced= ~Subject)

resLRT <- results(ddsLRT)

resOrdered <- resLRT[order(resLRT$log2FoldChange),]

before (August):

                           baseMean    log2FoldChange    lfcSE                 stat              pvalue         padj
1.2.23884.m1    6.44702248    -21.40753771    243.210067    22.99579811    1.62E-06    2.81E-05
1.2.17710.m1    6.409085561    -21.39654014    243.2647266    20.83683881    5.00E-06    7.52E-05
1.2.15201.m1    5.287495588    -21.13663198    244.5660123    17.45099162    2.95E-05    NA
1.2.16626.m1    3.685224198    -20.64942877    247.0622647    12.50468092    0.000405934    NA
1.2.23059.m1    1.842612099    -19.71611395    252.0642049    6.687733347    0.009707856    NA
1.2.20081.m1    1.760838826    -19.62249465    252.5823911    8.624799606    0.003316169    NA
1.2.17405.m1    1.522157821    -19.45937386    253.4939475    5.638415871    0.017571076    NA

Today:

                          baseMean          log2FoldChange    lfcSE            stat                pvalue           padj
1.2.23884.m1    6.44702248    -21.40611559    243.254081    1.541647987    0.214372672    NA
1.2.17710.m1    6.409085561    -21.39506324    243.3105312    1.09081661    0.296290345    NA
1.2.15201.m1    5.287495588    -21.13518038    244.6113169    1.058806029    0.30348766    NA
1.2.16626.m1    3.685224198    -20.64802728    247.1065218    0.99876432    0.317609691    NA
1.2.23059.m1    1.842612099    -19.71599394    252.0680863    3.763912397    0.052369906    NA
1.2.20081.m1    1.760838826    -19.62096702    252.593144    5.97981623    0.01447051    NA
1.2.17405.m1    1.522157821    -19.45922906    253.4986628    3.12559404    0.077071777    NA

My problem is that my old computer broke so I can not check the DESeq2 version that I was using two months ago, but it should have been up to date.

Do you have any idea of what is going on?

Thanks,

Catalina

deseq2 • 2.4k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 2 days ago
United States
Can you post your sessionInfo()?
ADD COMMENT
0
Entering edit mode
@catalina-aguilar-hurtado-6554
Last seen 4.1 years ago
United States

> sessionInfo()
R version 3.1.0 (2014-04-10)
Platform: x86_64-apple-darwin10.8.0 (64-bit)

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

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

other attached packages:
 [1] geneplotter_1.42.0        annotate_1.42.1           AnnotationDbi_1.26.1      lattice_0.20-29           Biobase_2.24.0           
 [6] sSeq_1.2.0                RColorBrewer_1.0-5        caTools_1.17.1            DESeq2_1.4.5              RcppArmadillo_0.4.450.1.0
[11] Rcpp_0.11.3               GenomicRanges_1.16.4      GenomeInfoDb_1.0.2        IRanges_1.22.10           BiocGenerics_0.10.0      

loaded via a namespace (and not attached):
 [1] bitops_1.0-6      DBI_0.3.1         genefilter_1.46.1 grid_3.1.0        locfit_1.5-9.1    RSQLite_0.11.4    splines_3.1.0    
 [8] stats4_3.1.0      survival_2.37-7   tools_3.1.0       XML_3.98-1.1      xtable_1.7-4      XVector_0.4.0    
 

ADD COMMENT
0
Entering edit mode

hi Catalina,

There was not a recent update, i.e. within the current release of DESeq2 v1.4, that would have changed results at all.

I would guess you were using DESeq2 v1.0 or 1.2 previously. It's hard to say what changed without knowing the version which you ran previously. If you installed R between April 2013 and Oct 2013, it would have likely been DESeq2 version 1.0, whereas if you installed R between Oct 2013 and April 2014, it would have likely been version 1.2. As 1.0 was the first release, there were many changes between 1.0 and 1.2. All the changes which affect end users are recorded in the NEWS file, but these might not be very informative for you.

Note that with test="LRT", there is not shrinkage of log fold changes, so sorting by log2 fold change will not necessarily be the genes with the smallest p-values (whereas with shrinkage of LFC these are typically a similar set to the genes with smallest p-value).

ADD REPLY
0
Entering edit mode

Hi Michael,

Thanks for your reply. But I guess my previous results will still be ok, depending on which genes of interest I will lose between the two different analysis?

As for the LRT test, what I understand from the manual, is to use LRT if I need to define a reduced formula? Which I should if I want to remove the subject effect.

ADD REPLY
0
Entering edit mode

It's not only the LRT that can account for the subject effect. You can also account for the subject effect using the standard steps which produce a Wald test. This is described in the section of the vignette: "Multi-factor designs".

design(dds) = ~ subject + condition
​dds = DESeq(dds)

 

ADD REPLY

Login before adding your answer.

Traffic: 764 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6