same padj for all the genes after DEseq analysis

Question

[aiswaryabioinfo]

Hi everyone,

I have done a pair comparison with DEseq2 to find differentially expressed genes between two samples with 6 replicates, for the DEseq2 result, i got exactly same padj value for all the genes and it is not significant, is this normal ? I don't think the padj should be the same for all the genes. I have attached the code I used. Please go through it and help me sort out this issue.

The paper which handled the same dataset has published that they could identify ~1000 DEGs for the same samples using edgeR and DESeq2 Package.

Then why i am getting this kind of a result

Thank you in advance.

library(DESeq2)

countdata <- read.table("trt_16_32_64_matrix.txt", header=TRUE, row.names=1)

countdata

countdata <- as.matrix(countdata)

countdata
(condition <- factor(c('exp_Day16','exp_Day16','exp_Day16','exp_Day16','exp_Day16','exp_Day16','exp_Day32','exp_Day32','exp_Day32','exp_Day32','exp_Day32','exp_Day32',
'exp_Day64','exp_Day64','exp_Day64','exp_Day64','exp_Day64','exp_Day64')))
(coldata <- data.frame(row.names=colnames(countdata), condition))
dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)
dds
dds <- dds[ rowSums(counts(dds)) > 10, ]
dds
dds <- DESeq(dds, minReplicatesForReplace=Inf)
rld <- rlogTransformation(dds)
vst <- vst(dds)
resultsNames(dds)
log2cutoff <- 2
qvaluecutoff <- 0.05
day32_16 <- results(dds, contrast=c("condition", "exp_Day32", "exp_Day16"), cooksCutoff=FALSE, independentFiltering=FALSE)
write.table(day32_16,file="Deseq_results32_vs_16.txt", quote = FALSE, sep="\t")
day32_16
day64_32 <- results(dds, contrast=c("condition", "exp_Day64", "exp_Day32"), cooksCutoff=FALSE, independentFiltering=FALSE)
write.table(day64_32,file="Deseq_results64_vs_32.txt", quote = FALSE, sep="\t")
day64_32
day64_16 <- results(dds, contrast=c("condition", "exp_Day64", "exp_Day16"), cooksCutoff=FALSE, independentFiltering=FALSE)
write.table(day64_16,file="Deseq_results64_vs_16.txt", quote = FALSE, sep="\t")
day64_16
normal<-assay(rld)
write.table(normal,file="rlog_normalcount.txt", quote = FALSE, sep="\t")
rlogMat <- normal - rowMeans(normal)
write.table(rlogMat,file="rlogMatcount.txt", quote = FALSE, sep="\t")
normal<-assay(vst)
write.table(normal,file="vst_normalcount.txt", quote = FALSE, sep="\t")
vstlogMat <- normal - rowMeans(normal)
write.table(rlogMat,file="vstlogMatcount.txt", quote = FALSE, sep="\t")
filter <- day32_16[which(day32_16$padj<0.05),]
write.table(filter,"Deseq_results_pvalue_005_day32_16.txt", quote=F,sep="\t")
filter <- day32_16[which(day32_16$padj<0.05 & day32_16$log2FoldChange>=2),]
write.table(filter,"Deseq_results_pvalue_005_FC+2_day32_16.txt", quote=F,sep="\t")
filter <- day32_16[which(day32_16$padj<0.05 & day32_16$log2FoldChange<=-2),]
write.table(filter,"Deseq_results_pvalue_005_FC-2_day32_16.txt", quote=F,sep="\t")
filter <- day32_16[which(day32_16$padj<0.05 & day32_16$log2FoldChange<=-2 | day32_16$padj<0.05 & day32_16$log2FoldChange>=2),]
write.table(filter,"Deseq_results_pvalue_005_FC2_day32_16.txt", quote=F,sep="\t")
filter <- day32_16[which(day32_16$padj<0.05 & day32_16$log2FoldChange<=-1.5 | day32_16$padj<0.05 & day32_16$log2FoldChange>=1.5),]
write.table(filter,"Deseq_results_pvalue_005_FC1.5_day32_16.txt", quote=F,sep="\t")
filter <- day64_32[which(day64_32$padj<0.05),]
write.table(filter,"Deseq_results_pvalue_005_day64_32.txt", quote=F,sep="\t")
filter <- day64_32[which(day64_32$padj<0.05 & day64_32$log2FoldChange>=2),]
write.table(filter,"Deseq_results_pvalue_005_FC+2_day64_32.txt", quote=F,sep="\t")
filter <- day64_32[which(day64_32$padj<0.05 & day64_32$log2FoldChange<=-2),]
write.table(filter,"Deseq_results_pvalue_005_FC-2_day64_32.txt", quote=F,sep="\t")
filter <- day64_32[which(day64_32$padj<0.05 & day64_32$log2FoldChange<=-2 | day64_32$padj<0.05 & day64_32$log2FoldChange>=2),]
write.table(filter,"Deseq_results_pvalue_005_FC2_day64_32.txt", quote=F,sep="\t")
filter <- day64_32[which(day64_32$padj<0.05 & day64_32$log2FoldChange<=-1.5 | day64_32$padj<0.05 & day64_32$log2FoldChange>=1.5),]
write.table(filter,"Deseq_results_pvalue_005_FC1.5_day64_32.txt", quote=F,sep="\t")
filter <- day64_16[which(day64_16$padj<0.05),]
write.table(filter,"Deseq_results_pvalue_005_day64_16.txt", quote=F,sep="\t")
filter <- day64_16[which(day64_16$padj<0.05 & day64_16$log2FoldChange>=2),]
write.table(filter,"Deseq_results_pvalue_005_FC+2_day64_16.txt", quote=F,sep="\t")
filter <- day64_16[which(day64_16$padj<0.05 & day64_16$log2FoldChange<=-2),]
write.table(filter,"Deseq_results_pvalue_005_FC-2_day64_16.txt", quote=F,sep="\t")
filter <- day64_16[which(day64_16$padj<0.05 & day64_16$log2FoldChange<=-2 | day64_16$padj<0.05 & day64_16$log2FoldChange>=2),]
write.table(filter,"Deseq_results_pvalue_005_FC2_day64_16.txt", quote=F,sep="\t")
filter <- day64_16[which(day64_16$padj<0.05 & day64_16$log2FoldChange<=-1.5 | day64_16$padj<0.05 & day64_16$log2FoldChange>=1.5),]
write.table(filter,"Deseq_results_pvalue_005_FC1.5_day64_16.txt", quote=F,sep="\t")

[Carlo Yague]

Can you show us the pvalue distribution (like a histogram) and the p.adj distribution for day32_16 ?

[aiswaryabioinfo]

This is the pvalue and padj I am getting for day 32_16

baseMean    log2FoldChange  lfcSE   stat    pvalue  padj
gene071790  0.00364470990823371 0   3.67603702902271    0   1   1
gene071910  2.52670110106324    -0.925869377753292  3.6735896523344 -0.252033968237292  0.801014808228701   1
gene071920  0.0017684908046795  0   3.67603703366083    0   1   1
gene072140  0.0132338542310079  -1.63447209000514   3.67524930232428    -0.444724141290627  0.656519120921536   1
gene072170  0.0233134553727202  -1.69569030301467   3.67543399998756    -0.461357843188154  0.644541891644605   1
gene072280  0.0090162636710289  -1.71864674299001   3.67550534890818    -0.467594678783569  0.640074470953938   1
gene072290  0.00282681108602329 0   3.67603703167104    0   1   1
gene072330  0.00396696145517866 0   3.67603702802475    0   1   1
gene072340  1.10232459596773    -1.29929603334227   1.69869436533412    -0.764879227162627  0.444343464578242   1
gene072350  0.0946211792117389  -0.925869002733872  3.67358965348958    -0.252033866072761  0.801014887195908   1
gene072360  0.0225718617774871  0   3.67603698364365    0   1   1
gene072380  0.382201543382355   -2.46595641962353   2.34574179393414    -1.05124802141491   0.293144693024846   1
gene072390  0.00915777006668266 0   3.67603701205712    0   1   1
gene072420  0.727765396747566   -1.07434573556786   1.30955623458107    -0.820389157180061  0.411994294679308   1
gene072430  0.185235533353433   0   3.67603695271477    0   1   1
gene072460  0.00556966720967944 -1.83503712195606   3.67588448263864    -0.499209681540053  0.617631674812785   1
gene072500  0.227903343157695   -2.42970097182479   2.63053686527404    -0.923652127403915  0.355667464392678   1
gene072540  0.00323182607045902 0   3.67603703034583    0   1   1

Answer 1

[Carlo Yague]

The pvalue distribution is a bit weird, somewhere in between the uniform and conservative types, meaning that there are very few significant hits. So after multiple-testing correction, it is normal that all the p.adj increase to 1 (the highest value), so no significantly differentially expressed genes.

Why do you have much lower power than in the published analysis (1000 DEG) ? I can not say, but you should check if you included all the data properly, and also did not mixed or mislabelled the conditions. For instance you could make a PCA to see if your samples cluster in an unexpected way.

[aiswaryabioinfo]

The author has used edgeR package etLRT test where as I am using DESeq2 wald test. Is the difference in DEGs is due to this

[Carlo Yague]

It is unlikely that such a big difference comes from this. That being said, you can use LRT in DESeq2 too.