same padj for all the genes after DEseq analysis
1
0
Entering edit mode
@aiswaryabioinfo-15154
Last seen 3.3 years ago

I have done a pair comparison with DEseq2 (wald test) 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.

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

Then why i am getting this kind of a result. The author has used edgeR package etLRT test whereas I am using DESeq2 wald test. Is the difference in DEGs is due to the usage of two packages and two type of test

Please go through it and help me sort out this issue.

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)
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

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_India.1252  LC_CTYPE=English_India.1252    LC_MONETARY=English_India.1252 LC_NUMERIC=C                  
[5] LC_TIME=English_India.1252    

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DESeq2_1.26.0               SummarizedExperiment_1.16.0 DelayedArray_0.12.0         BiocParallel_1.20.0        
 [5] matrixStats_0.55.0          Biobase_2.46.0              GenomicRanges_1.38.0        GenomeInfoDb_1.22.0        
 [9] IRanges_2.20.1              S4Vectors_0.24.0            BiocGenerics_0.32.0        

loaded via a namespace (and not attached):
 [1] bit64_0.9-7            splines_3.6.1          Formula_1.2-3          assertthat_0.2.1       latticeExtra_0.6-29   
 [6] blob_1.2.1             GenomeInfoDbData_1.2.2 pillar_1.4.6           RSQLite_2.1.2          backports_1.1.5       
[11] lattice_0.20-38        glue_1.3.1             digest_0.6.22          RColorBrewer_1.1-2     XVector_0.26.0        
[16] checkmate_1.9.4        colorspace_1.4-1       htmltools_0.4.0        Matrix_1.2-17          XML_3.98-1.20         
[21] pkgconfig_2.0.3        genefilter_1.68.0      zlibbioc_1.32.0        purrr_0.3.3            xtable_1.8-4          
[26] scales_1.0.0           jpeg_0.1-8.1           htmlTable_2.0.1        tibble_2.1.3           annotate_1.64.0       
[31] ggplot2_3.3.2          ellipsis_0.3.1         nnet_7.3-12            survival_3.2-3         magrittr_1.5          
[36] crayon_1.3.4           memoise_1.1.0          foreign_0.8-71         tools_3.6.1            data.table_1.13.0     
[41] lifecycle_0.2.0        stringr_1.4.0          locfit_1.5-9.1         munsell_0.5.0          cluster_2.1.0         
[46] AnnotationDbi_1.48.0   compiler_3.6.1         rlang_0.4.7            grid_3.6.1             RCurl_1.95-4.12       
[51] rstudioapi_0.11        htmlwidgets_1.5.1      bitops_1.0-6           base64enc_0.1-3        gtable_0.3.0          
[56] DBI_1.1.0              R6_2.4.1               gridExtra_2.3          knitr_1.29             dplyr_0.8.3           
[61] bit_1.1-14             Hmisc_4.3-0            stringi_1.4.3          Rcpp_1.0.2             geneplotter_1.64.0    
[66] vctrs_0.3.2            rpart_4.1-15           acepack_1.4.1          png_0.1-7              tidyselect_0.2.5      
[71] xfun_0.16             
>

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

gene072340  1.10232459596773    -1.29929603334227   1.69869436533412    -0.764879227162627  0.444343464578242   1

gene072350  0.0946211792117389  -0.925869002733872  3.67358965348958    -0.252033866072761  0.801014887195908   1

gene072380  0.382201543382355   -2.46595641962353   2.34574179393414    -1.05124802141491   0.293144693024846   1

gene072420  0.727765396747566   -1.07434573556786   1.30955623458107    -0.820389157180061  0.411994294679308   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
etLRT WALD edgeR DESEq • 1.9k views
ADD COMMENT
0
Entering edit mode

Cross posted to Biostars https://www.biostars.org/p/9484960/

ADD REPLY
0
Entering edit mode
@mikelove
Last seen 1 day ago
United States

Answer from the forum:

adj.P metod = "BH" many exact the same value

ADD COMMENT

Login before adding your answer.

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