same padj for all the genes after DEseq analysis
1
0
Entering edit mode
10 months ago

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")
NGS DEGs DESeq2 RNA-seq padj • 810 views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

P value distribution

Padj distribution

ADD REPLY
1
Entering edit mode
10 months ago

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.

ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

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