Search
Question: All same padj of contrast results in the DESeq2
0
gravatar for jjinhyoungkim
3.3 years ago by
Canada
jjinhyoungkim20 wrote:

Hi,

I am trying to run DESeq2 with the duplicated raw counts from 4 groups and 2 conditions. In the comparison using "result ( contrast)", some comparisons had all same padj (0.9999874). Is it normal or anything wrong? FYI, the code are below. Thank you in advance for your time.

>library("DESeq2")
>countMatrix = read.table("~/Desktop/aaa.txt",header = T,row.names = 1)
>coldata = data.frame(row.names =colnames(countMatrix),group =rep(c("gt1","gt2","gt3","gt4"),2,each = 2),treatment = rep(c("control","treated"),each= 8))
>coldata$treatment = factor(x = coldata$treatment,levels = c('control','treated'))>dds = DESeqDataSetFromMatrix(countData =countMatrix, colData = coldata, design = ~ group + treatment +group:treatment)
>dds = DESeq(dds)
> resultsNames(dds)
 [1] "Intercept"       "groupgt1"       "groupgt2"       "groupgt3"                 [5] "groupgt4"   "treatmentcontrol"    "treatmenttreated" "groupgt1.treatmentcontrol"
 [9] "groupgt2.treatmentcontrol" "groupgt3.treatmentcontrol" "groupgt4.treatmentcontrol" "groupgt1.treatmenttreated"
[13] "groupgt2.treatmenttreated" "groupgt3.treatmenttreated" "groupgt4.treatmenttreated"
> res <-results(dds, contrast=c("treatment", "treated", "control"))
> res

log2 fold change (MAP): treatment treated vs control
Wald test p-value: treatment treated vs control
DataFrame with 35584 rows and 6 columns
                 baseMean log2FoldChange     lfcSE      stat       pvalue         padj
                <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric>
unigene22422654  82.11115       2.753315 0.2880473  9.558550 1.194176e-21 1.863751e-17
unigene22392883 110.42190       2.742819 0.3119236  8.793239 1.453092e-18 1.133920e-14
unigene22410436  94.91773       2.054046 0.2427170  8.462723 2.612059e-17 1.306140e-13
unigene22398806 133.11805       2.269261 0.2690689  8.433751 3.347576e-17 1.306140e-13
unigene22391050 159.42752       1.589496 0.1900555  8.363325 6.097477e-17 1.903266e-13
...                   ...            ...       ...       ...          ...          ...
unigene22403397         0             NA        NA        NA           NA           NA
unigene22417358         0             NA        NA        NA           NA           NA
unigene22397386         0             NA        NA        NA           NA           NA
unigene22415322         0             NA        NA        NA           NA           NA
unigene22402022         0             NA        NA        NA           NA           NA

> res1 <- results(dds, contrast=list("groupgt2.treatmenttreated", "groupgt2.treatmentcontrol"))

> res1
log2 fold change (MAP): groupgt2.treatmenttreated vs groupgt2.treatmentcontrol
Wald test p-value: groupgt2.treatmenttreated vs groupgt2.treatmentcontrol
DataFrame with 35584 rows and 6 columns
                 baseMean log2FoldChange     lfcSE       stat    pvalue      padj
                <numeric>      <numeric> <numeric>  <numeric> <numeric> <numeric>
unigene22408676 23.641057     -0.2285822 0.5213050 -0.4384808 0.6610378 0.9999874
unigene22408675  1.032055     -0.6212150 1.2989812 -0.4782325 0.6324848 0.9999874
unigene22408674 11.503303     -0.1884673 0.5643655 -0.3339456 0.7384206 0.9999874
unigene22408673  1.849157      1.5855804 1.1773601  1.3467251 0.1780688 0.9999874
unigene22408672  2.717044     -1.2518670 1.0399068 -1.2038261 0.2286568 0.9999874
...                   ...            ...       ...        ...       ...       ...
unigene22403397         0             NA        NA         NA        NA        NA
unigene22417358         0             NA        NA         NA        NA        NA
unigene22397386         0             NA        NA         NA        NA        NA
unigene22415322         0             NA        NA         NA        NA        NA
unigene22402022         0             NA        NA         NA        NA        NA

ADD COMMENTlink modified 3.2 years ago by Michael Love20k • written 3.3 years ago by jjinhyoungkim20
4
gravatar for Michael Love
3.2 years ago by
Michael Love20k
United States
Michael Love20k wrote:

Repeated values for the adjusted p-value is a direct consequence of the method of Benjamini and Hochberg. The paper is quite accessible and worth a look:

http://www.jstor.org/stable/2346101

See formula (1) under "False Discovery Rate Controlling Procedure".

All the p-values p_i up to p_k also get rejected at the FDR q.

For a visualization of this formula, take a few sorted p-values with adjusted p-value written above:

p <- sort(c(.01,.2,.21,.22,.5,.51,.52,.8,.9))
padj <- p.adjust(p, method="BH")
plot(p,ylim=c(0,1),xlim=c(0,length(p)))
text(seq_along(p),p,round(padj,3),pos=3)
abline(0,padj[4]/length(p))

To translate from the formula to the code above, k = 4, m=length(p), and q = padj[4].

ADD COMMENTlink modified 3.2 years ago • written 3.2 years ago by Michael Love20k
0
gravatar for Ryan C. Thompson
3.3 years ago by
The Scripps Research Institute, La Jolla, CA
Ryan C. Thompson7.0k wrote:

If all the adjusted p-values are 1, then that means the algorithm didn't find any significantly differentially expressed genes for that contrast.

ADD COMMENTlink written 3.3 years ago by Ryan C. Thompson7.0k
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 2.2.0
Traffic: 121 users visited in the last hour