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")
•
written
by
aiswaryabioinfo •
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.
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
It is unlikely that such a big difference comes from this. That being said, you can use LRT in DESeq2 too.
Log in to answer this question.
Can you show us the pvalue distribution (like a histogram) and the p.adj distribution for
day32_16?This is the pvalue and padj I am getting for day 32_16