Question

Deseq2 padj value very low

0

Entering edit mode

jarod_v6@libero.it ▴ 40

@jarod_v6liberoit-6654

Last seen 5.7 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 • 2.0k views

ADD COMMENT • link updated 6.7 years ago by Michael Love 42k • written 6.7 years ago by jarod_v6@libero.it ▴ 40

score 0 · Answer 1 · 2018-02-12

0

Entering edit mode

Michael Love 42k

@mikelove

Last seen 18 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 • link 6.7 years ago Michael Love 42k

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 • link 6.7 years ago jarod_v6@libero.it ▴ 40

0

Entering edit mode

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

ADD REPLY • link 6.7 years ago Michael Love 42k

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 • link 6.7 years ago jarod_v6@libero.it ▴ 40