Deseq2 padj value very low
1
0
Entering edit mode
@jarod_v6liberoit-6654
Last seen 5.1 years ago
Italy

I have perform differential analysis using deseq2.  In the result I have many genes have very low padj. It is possible?

 

> res[which.min(res$padj),]
log2 fold change (MAP): condition Myoepithelioma vs EMC
Wald test p-value: condition Myoepithelioma vs EMC
DataFrame with 1 row and 11 columns
                 baseMean log2FoldChange     lfcSE      stat       pvalue         padj
                <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
ENSG00000150361  1289.267      -12.98314 0.7835716 -16.56919 1.163906e-61 1.875285e-57
                        ensembl    entrez hgnc_symbol  chromosome        band
                    <character> <integer> <character> <character> <character>
ENSG00000150361 ENSG00000150361     57626       KLHL1          13      q21.33

subset(res,res$padj< 1*10^-20)
log2 fold change (MAP): condition Myoepithelioma vs EMC
Wald test p-value: condition Myoepithelioma vs EMC
DataFrame with 4 rows and 11 columns
                 baseMean log2FoldChange     lfcSE      stat       pvalue         padj
                <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
ENSG00000150361 1289.2669     -12.983144 0.7835716 -16.56919 1.163906e-61 1.875285e-57
ENSG00000184905 1047.7598       9.413017 0.9076872  10.37033 3.383477e-25 1.362864e-21
ENSG00000196850  972.6995      -4.283253 0.3794026 -11.28947 1.479300e-29 7.944829e-26
ENSG00000197696  482.4788      -7.368074 0.5813144 -12.67485 8.151262e-37 6.566656e-33
                        ensembl    entrez hgnc_symbol  chromosome        band
                    <character> <integer> <character> <character> <character>
ENSG00000150361 ENSG00000150361     57626       KLHL1          13      q21.33
ENSG00000184905 ENSG00000184905    140597      TCEAL2           X       q22.1
ENSG00000196850 ENSG00000196850    160760       PPTC7          12      q24.11
ENSG00000197696 ENSG00000197696      4828         NMB          15       q25.3

 

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.3 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=it_IT.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=it_IT.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=it_IT.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=it_IT.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] gplots_3.0.1               genefilter_1.56.0          limma_3.30.13             
 [4] biomaRt_2.30.0             reshape2_1.4.3             RColorBrewer_1.1-2        
 [7] ggplot2_2.2.1              pheatmap_1.0.8             DESeq2_1.14.1             
[10] SummarizedExperiment_1.4.0 Biobase_2.34.0             GenomicRanges_1.26.4      
[13] GenomeInfoDb_1.10.3        IRanges_2.8.2              S4Vectors_0.12.2          
[16] BiocGenerics_0.20.0        RTCGAToolbox_2.4.0        

loaded via a namespace (and not attached):
 [1] tidyr_0.7.2          bit64_0.9-7          splines_3.3.2        gtools_3.5.0        
 [5] Formula_1.2-2        assertthat_0.2.0     latticeExtra_0.6-28  blob_1.1.0          
 [9] pillar_1.1.0         RSQLite_2.0          backports_1.1.2      lattice_0.20-35     
[13] glue_1.2.0           digest_0.6.14        XVector_0.14.1       checkmate_1.8.5     
[17] QoRTs_1.1.8          colorspace_1.3-2     htmltools_0.3.6      Matrix_1.2-10       
[21] plyr_1.8.4           XML_3.98-1.9         pkgconfig_2.0.1      zlibbioc_1.20.0     
[25] purrr_0.2.4          xtable_1.8-2         RCircos_1.2.0        scales_0.5.0        
[29] gdata_2.18.0         BiocParallel_1.8.2   htmlTable_1.11.0     tibble_1.4.1        
[33] annotate_1.52.1      nnet_7.3-12          lazyeval_0.2.1       survival_2.41-3     
[37] RJSONIO_1.3-0        magrittr_1.5         memoise_1.1.0        foreign_0.8-69      
[41] tools_3.3.2          data.table_1.10.4-3  stringr_1.2.0        munsell_0.4.3       
[45] locfit_1.5-9.1       cluster_2.0.6        AnnotationDbi_1.36.2 bindrcpp_0.2        
[49] caTools_1.17.1       rlang_0.1.6          grid_3.3.2           RCurl_1.95-4.10     
[53] rstudioapi_0.7       htmlwidgets_0.9      labeling_0.3         bitops_1.0-6        
[57] base64enc_0.1-3      gtable_0.2.0         DBI_0.7              R6_2.2.2            
[61] gridExtra_2.3        knitr_1.18           dplyr_0.7.4          bit_1.1-12          
[65] bindr_0.1            Hmisc_4.1-1          KernSmooth_2.23-15   stringi_1.1.6       
[69] Rcpp_0.12.14         geneplotter_1.52.0   rpart_4.1-12         acepack_1.4.1
deseq2 • 1.8k views
ADD COMMENT
0
Entering edit mode
@mikelove
Last seen 3 hours ago
United States

A low adjusted p-value is not a problem. It means the null can be rejected among this and genes with lower adjusted p-values. If you don't like very small p-values, you can test against a banded null, instead of a point null at LFC=0, which is often trivial to reject for well-powered experiments. See the DESeq2 paper, and the vignette section on the use of positive lfcThreshold.

ADD COMMENT
0
Entering edit mode

thanks for the help..

the problem that I  have all the data are in that situation...

 

 

 

 

res<-results(dds)
> head(res)
log2 fold change (MLE): Intercept
Wald test p-value: Intercept
DataFrame with 6 rows and 6 columns
                 baseMean log2FoldChange     lfcSE      stat    pvalue
                <numeric>      <numeric> <numeric> <numeric> <numeric>
ENSG00000000003     132.3           7.05    0.1643      42.9   0.0e+00
ENSG00000000419      80.1           6.32    0.0991      63.8   0.0e+00
ENSG00000000457     109.2           6.77    0.0622     108.8   0.0e+00
ENSG00000000460      50.8           5.67    0.0931      60.8   0.0e+00
ENSG00000000938      22.6           4.50    0.1746      25.8  2.2e-146
ENSG00000000971    1729.9          10.76    0.2008      53.6   0.0e+00
                     padj
                <numeric>
ENSG00000000003  0.00e+00
ENSG00000000419  0.00e+00
ENSG00000000457  0.00e+00
ENSG00000000460  0.00e+00
ENSG00000000938 3.28e-146
ENSG00000000971  0.00e+00
> plotCounts(dds,gene=which.min(res$padj))
> res["ENSG00000000003",]
log2 fold change (MLE): Intercept
Wald test p-value: Intercept
DataFrame with 1 row and 6 columns
                 baseMean log2FoldChange     lfcSE      stat    pvalue
                <numeric>      <numeric> <numeric> <numeric> <numeric>
ENSG00000000003       132           7.05     0.164      42.9         0
                     padj
                <numeric>
ENSG00000000003         0
> summary(res)

out of 19632 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)     : 19624, 100%
LFC < 0 (down)   : 0, 0%
outliers [1]     : 0, 0%
low counts [2]   : 0, 0%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.3 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=it_IT.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=it_IT.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=it_IT.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=it_IT.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] genefilter_1.56.0          rafalib_1.0.0             
 [3] ggplot2_2.2.1              limma_3.30.13             
 [5] RColorBrewer_1.1-2         gplots_3.0.1              
 [7] org.Hs.eg.db_3.4.0         annotate_1.52.1           
 [9] XML_3.98-1.9               AnnotationDbi_1.36.2      
[11] biomaRt_2.30.0             pheatmap_1.0.8            
[13] tximportData_1.2.0         tximport_1.2.0            
[15] DESeq2_1.14.1              SummarizedExperiment_1.4.0
[17] Biobase_2.34.0             GenomicRanges_1.26.4      
[19] GenomeInfoDb_1.10.3        IRanges_2.8.2             
[21] S4Vectors_0.12.2           BiocGenerics_0.20.0       
[23] RTCGAToolbox_2.4.0        

loaded via a namespace (and not attached):
 [1] bit64_0.9-7         splines_3.3.2       gtools_3.5.0       
 [4] Formula_1.2-2       latticeExtra_0.6-28 blob_1.1.0         
 [7] pillar_1.1.0        RSQLite_2.0         backports_1.1.2    
[10] lattice_0.20-35     digest_0.6.15       XVector_0.14.1     
[13] checkmate_1.8.5     QoRTs_1.1.8         colorspace_1.3-2   
[16] htmltools_0.3.6     Matrix_1.2-12       plyr_1.8.4         
[19] pkgconfig_2.0.1     zlibbioc_1.20.0     xtable_1.8-2       
[22] RCircos_1.2.0       scales_0.5.0        gdata_2.18.0       
[25] BiocParallel_1.8.2  htmlTable_1.11.2    tibble_1.4.2       
[28] nnet_7.3-12         lazyeval_0.2.1      survival_2.41-3    
[31] RJSONIO_1.3-0       magrittr_1.5        memoise_1.1.0      
[34] foreign_0.8-69      tools_3.3.2         data.table_1.10.4-3
[37] stringr_1.2.0       munsell_0.4.3       locfit_1.5-9.1     
[40] cluster_2.0.6       caTools_1.17.1      rlang_0.1.6        
[43] grid_3.3.2          RCurl_1.95-4.10     rstudioapi_0.7     
[46] htmlwidgets_1.0     labeling_0.3        bitops_1.0-6       
[49] base64enc_0.1-3     gtable_0.2.0        DBI_0.7            
[52] gridExtra_2.3       knitr_1.19          bit_1.1-12         
[55] Hmisc_4.1-1         KernSmooth_2.23-15  stringi_1.1.6      
[58] Rcpp_0.12.15        geneplotter_1.52.0  rpart_4.1-12       
[61] acepack_1.4.1 
ADD REPLY
0
Entering edit mode

You are testing the Intercept coefficient it says above. Use a design eg ~condition

ADD REPLY
0
Entering edit mode

thanks It was an error on the preparation on DESeqDataSetFromHTSeqCount so they not change that parameter.

Thanks s much!

ADD REPLY

Login before adding your answer.

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