Question: Deseq2 padj value very low
0
gravatar for jarod_v6@libero.it
18 months ago by
Italy
jarod_v6@libero.it40 wrote:

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 • 429 views
ADD COMMENTlink modified 18 months ago by Michael Love24k • written 18 months ago by jarod_v6@libero.it40
Answer: Deseq2 padj value very low
0
gravatar for Michael Love
18 months ago by
Michael Love24k
United States
Michael Love24k wrote:

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 COMMENTlink written 18 months ago by Michael Love24k

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 REPLYlink written 18 months ago by jarod_v6@libero.it40

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

ADD REPLYlink written 18 months ago by Michael Love24k

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

Thanks s much!

ADD REPLYlink written 18 months ago by jarod_v6@libero.it40
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 283 users visited in the last hour