DESeq2 ANODEV with 3 sample groups
2
1
Entering edit mode
jenniewoo ▴ 70
@jenniewoo-8959
Last seen 11 months ago
United States

I'm trying to use DESeq2 to perform something like a one-way ANOVA on 3
groups of samples. I've seen the same question and answer here:

DESeq2 One-way ANOVA with 3 sample groups?

However, I can't seem to  get a single p value. I have my data input successfully, and I've been able to
perform pairwise comparisons using DESeq2 but I haven't been able to
extract one p-value across 3 different groups using Mike's kind reply code in the above post. Thanks for your help.

The code I used was 

group=factor(c("A","A","B","B","C","C"))
design(dds) = ~ group
dds = DESeq(dds, test = "LRT", reduced = ~ 1)
res=results(dds)

The still give the pairwise comparison, default C vs A p values. I can use contrast to get other pairwise p values but Is there anyway I can get a single p value per gene?

>sessionInfo()

R version 3.2.1 (2015-06-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252    LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                           LC_TIME=English_United States.1252    

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

other attached packages:
 [1] ggplot2_1.0.1             gplots_2.17.0             RColorBrewer_1.1-2        limma_3.24.15             BiocInstaller_1.18.4     
 [6] DESeq2_1.8.1              RcppArmadillo_0.5.400.2.0 Rcpp_0.12.0               GenomicRanges_1.20.6      GenomeInfoDb_1.4.2       
[11] IRanges_2.2.7             S4Vectors_0.6.5           BiocGenerics_0.14.0      

deseq2 • 8.7k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 46 minutes ago
United States

"This still gives the pairwise comparison, default C vs A p values"

No, the p-value here is for the comparison of the full model and the reduced model, so all 3 levels.

For the LRT, there is only one p-value per gene and it is the one you are getting back from results(). You will notice when you print the object that the LFC is associated with a comparison, but the p-value says "~condition" vs "~1". So is the visual indication to the user that the p-value is a LRT p-value, not a pairwise comparison of 2 levels. You can see either LFC by using the 'name' argument. But this doesn't change the fact that the p-value you see is associated with all 3 levels.

Please find the paragraph in the Details section in ?results about likelihood ratio test p-values.

ADD COMMENT
0
Entering edit mode

Thanks very much for the kind reply! It is very helpful. You are right. I did see that the object returned by results() says "LTR p-value '~ condition' vs '~ 1'. That is clear. The log2 fold change (MLE) is C vs A. My question is, if I do res2=results(dds, contrast=c("condition", "B", "A")), I got another set of p values and LFC for B vs. A. same for C vs B.

For ANODEV single p value comparing all three groups at the same time, which p value should I report?

 

> head(res)
log2 fold change (MLE): condition C vs A 
LRT p-value: '~ condition' vs '~ 1' 
DataFrame with 63677 rows and 6 columns
                    baseMean log2FoldChange     lfcSE      stat     pvalue      padj
                   <numeric>      <numeric> <numeric> <numeric>  <numeric> <numeric>
ENSG00000273423 9.268989e-02   -0.075613392 5.2677216 0.8703565  0.8325749        NA
ENSG00000110514 1.503328e+03    0.284893906 0.1693066 6.3021485  0.0978005 0.3797831
ENSG00000268358 1.337136e+01   -0.313182594 0.7293037 0.4572943  0.9281683        NA
ENSG00000086015 1.425711e+03    0.045650942 0.1450839 1.8636136  0.6011907 0.8423642
ENSG00000272373 1.925744e+01   -0.008622687 0.5492149 0.3001698  0.9599965 0.9888905
...                      ...            ...       ...       ...        ...       ...
ADD REPLY
0
Entering edit mode
@mikelove
Last seen 46 minutes ago
United States

"My question is, if I do res2=results(dds, contrast=c("condition", "B", "A")), I got another set of p values and LFC for B vs. A. same for C vs B."

I believe that you get a new LFC but the p-values should be exactly the same and they should say "LRT p-value...". Can you check this?

 

ADD COMMENT
1
Entering edit mode

Thanks Mike for the prompt reply. Yes I checked and the baseMean, test stats, pvalue and padj are all the same! That makes sense. Thank you so much for your help! Our community benefit tremendously from your responses. 

ADD REPLY
0
Entering edit mode

Hi Michael,  it is easy to understand that the LRT p-value is for the comparison of the 3 levels. However, how can we retrieve the p-value for the comparison of condition B vs A? 

Thank you!

ADD REPLY
0
Entering edit mode

You can tell results() to generate a Wald test p-value for a contrast with results(dds, ..., test="Wald"). See ?results for more.

ADD REPLY
0
Entering edit mode

Thanks Michael. This is exactly the answer I wanted!

ADD REPLY

Login before adding your answer.

Traffic: 706 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