Limma adj.p.value the same for many different p.values
1
0
Entering edit mode
walkerrh95 • 0
@a92e11c5
Last seen 7 weeks ago
United Kingdom

I have been using limma to analyse some statistical comparisons, however I have then checked my data and found that instead of proper adj.p.values many different p.values have the same adj.p.value.

Here is the head of the data that I am using and the problem that I am getting (I couldn't work out how to upload an attachment with it but would if I could.


Output:

Identifier logFC AveExpr t P.Value adj.P.Val B log.P.Value log.adj.P.Value Threshold Threshold_adj.P.Val pvaluecheck qvalue Sample 1A Sample1B Sample1C Sample2B Sample2C ID 2734 0.691697932 17.32119982 2.293548327 0.049347051 0.2588444 -4.052126814 1.306738799 0.586961226 Sig Increased Not Significant 0.049347051 0.174471411 17.16420046 16.22682986 16.92039146 17.2788217 17.64552202 ID 4852 0.640589147 19.90720357 2.294737142 0.049253556 0.2588444 -4.050361884 1.307562411 0.586961226 Sig Increased Not Significant 0.049253556 0.174471411 19.16438921 19.80824867 19.43875164 19.82073962 20.40136502 ID 5496 1.580758704 16.5507023 2.987776287 0.048986077 0.2588444 -3.462916788 1.309927335 0.586961226 Sig Increased Not Significant 0.048986077 0.174471411 15.36993643 NA 16.08303656 NA 17.3072452 ID 5926 0.669486048 16.52943882 2.352074189 0.048847084 0.2588444 -3.934596674 1.311161359 0.586961226 Sig Increased Not Significant 0.048847084 0.174471411 16.13102241 NA 16.47239228 17.17680579 16.765581 ID 876 0.981091651 16.20035906 2.301930364 0.048691575 0.2588444 -4.039679151 1.312546179 0.586961226 Sig Increased Not Significant 0.048691575 0.174471411 15.83553995 15.32407361 15.80576499 17.3689849 15.90345077 ID 1736 1.120753382 15.60917806 2.702128465 0.048328399 0.258734098 -3.586944058 1.315797593 0.587146332 Sig Increased Not Significant 0.048328399 0.174397063 14.64995302 NA NA 15.93223475 15.60917806 ID 6504 1.725697353 16.0821604 2.538702326 0.04816353 0.258356442 -3.677392124 1.317281694 0.587780704 Sig Increased Not Significant 0.04816353 0.174142508 16.21553825 NA 16.02509432 17.84601364 NA ID 3532 0.61867589 15.85210422 2.316599176 0.04756518 0.257576777 -4.017875669 1.322710851 0.589093295 Sig Increased Not Significant 0.04756518 0.173616982 15.81271282 15.20612771 15.62447851 16.10525568 16.22764213 ID 1082 0.824659566 17.43663583 2.744051677 0.046140154 0.255455613 -3.542963732 1.335920962 0.592684551 Sig Increased Not Significant 0.046140154 0.172187233 16.65630453 16.8713518 NA 17.58848773 NA


design_paired_group1 = model.matrix(~0+Group, data=metadata)
#corfit=duplicateCorrelation(log2intensity.filtered, design=design_paired_group1)

## fit linear model to the data
fit1_paired = lmFit(batch.patient, design = design_paired_group1)

## pairwise comparison between control vs treatment
cont_paired <- makeContrasts(GroupTreatment-GroupControl, levels = design_paired_group1)
fit2_paired = contrasts.fit(fit1_paired, contrasts = cont_paired)

## This is the final results array- the standard errors are moderated using a simple empirical Bayes model using eBayes. The eBayes command will compute the consensus pooled variance, and then use it to compute the empirical Bayes (moderated) pooled variance for each gene. This also adjusts the degrees of freedom for the contrast t-tests.The command also computes the t-tests and associated p-values.


fit3_paired <- eBayes(fit2_paired)
finalresult_paired <- topTable(fit3_paired, adjust.method = "BH", number = Inf)

setDT(finalresult_paired, keep.rownames ="Identifier")
setDTControlvsTreatment, keep.rownames ="Identifier")
names(finalresult_paired)[1]<-"Identifier"
finalresult_paired$log.P.Value = -log10(finalresult_paired$P.Value)
finalresult_paired$log.adj.P.Value = -log10(finalresult_paired$adj.P.Val)
merged_finalresult <- merge(finalresult_paired,Protein_Names,by="Identifier")# merge with original data set by Identifier to have all additional info e.g. protein names, etc.
merged_finalresult_ControlvsTreatment <- merge(merged_finalresult,ControlvsTreatment, by="Identifier") # merge with normalised log2 intensities that were used for the differential testing

write.table(merged_finalresult_ControlvsTreatment,"mergedfinalresult_paired_ControlvsTreatment.txt",sep = "\t",
            row.names = F,quote=F)

sig_group1 <- finalresult_paired[merged_finalresult_ControlvsTreatment$P.Value <0.05,]
write.table(sig_group1, file = "sigProts_ControlvsTreatment.txt", col.names = TRUE, sep = "\t")

Metadata:

Sample Group Bio Sample1A Control 1 Sample1B Control 2 Sample1C Control 3 Sample2B Treatment 2

Sample2C Treatment 3

sessioninfo()

R version 4.2.1 (2022-06-23 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22000)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.utf8  LC_CTYPE=English_United Kingdom.utf8   
[3] LC_MONETARY=English_United Kingdom.utf8 LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.utf8    

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

other attached packages:
 [1] rgl_0.110.2            plotly_4.10.0          ggfortify_0.4.14      
 [4] factoextra_1.0.7       pRolocdata_1.34.0      MSnbase_2.22.0        
 [7] ProtGenerics_1.28.0    S4Vectors_0.34.0       mzR_2.30.0            
[10] Rcpp_1.0.9             Biobase_2.56.0         BiocGenerics_0.42.0   
[13] ggforce_0.3.4          EnhancedVolcano_1.14.0 qvalue_2.28.0         
[16] pheatmap_1.0.12        ReactomePA_1.40.0      pwr_1.3-0             
[19] plyr_1.8.7             ggrepel_0.9.1          yaml_2.3.5            
[22] statmod_1.4.37         kableExtra_1.3.4       ggridges_0.5.4        
[25] clusterProfiler_4.4.4  edgeR_3.38.4           ggpubr_0.4.0          
[28] ggsci_2.9              RColorBrewer_1.1-3     BBmisc_1.12           
[31] gplots_3.1.3           forcats_0.5.2          stringr_1.4.1         
[34] dplyr_1.0.10           purrr_0.3.4            readr_2.1.2           
[37] tidyr_1.2.1            tibble_3.1.8           tidyverse_1.3.2       
[40] data.table_1.14.2      corrplot_0.92          ggplot2_3.3.6         
[43] BiocManager_1.30.18    limma_3.52.3          

loaded via a namespace (and not attached):
  [1] utf8_1.2.2             tidyselect_1.1.2       RSQLite_2.2.17        
  [4] AnnotationDbi_1.58.0   htmlwidgets_1.5.4      grid_4.2.1            
  [7] BiocParallel_1.30.3    scatterpie_0.1.8       munsell_0.5.0         
 [10] preprocessCore_1.58.0  codetools_0.2-18       withr_2.5.0           
 [13] colorspace_2.0-3       GOSemSim_2.22.0        knitr_1.40            
 [16] rstudioapi_0.14        ggsignif_0.6.3         DOSE_3.22.1           
 [19] mzID_1.34.0            labeling_0.4.2         GenomeInfoDbData_1.2.8
 [22] polyclip_1.10-0        bit64_4.0.5            farver_2.1.1          
 [25] downloader_0.4         vctrs_0.4.1            treeio_1.20.2         
 [28] generics_0.1.3         xfun_0.32              doParallel_1.0.17     
 [31] R6_2.5.1               GenomeInfoDb_1.32.4    clue_0.3-61           
 [34] graphlayouts_0.8.1     locfit_1.5-9.6         MsCoreUtils_1.8.0     
 [37] bitops_1.0-7           cachem_1.0.6           fgsea_1.22.0          
 [40] gridGraphics_0.5-1     assertthat_0.2.1       scales_1.2.1          
 [43] ggraph_2.0.6           enrichplot_1.16.2      googlesheets4_1.0.1   
 [46] gtable_0.3.1           affy_1.74.0            tidygraph_1.2.2       
 [49] rlang_1.0.6            systemfonts_1.0.4      splines_4.2.1         
 [52] rstatix_0.7.0          lazyeval_0.2.2         impute_1.70.0         
 [55] gargle_1.2.1           broom_1.0.1            checkmate_2.1.0       
 [58] reshape2_1.4.4         abind_1.4-5            modelr_0.1.9          
 [61] crosstalk_1.2.0        backports_1.4.1        tools_4.2.1           
 [64] ggplotify_0.1.0        affyio_1.66.0          ellipsis_0.3.2        
 [67] base64enc_0.1-3        zlibbioc_1.42.0        RCurl_1.98-1.8        
 [70] viridis_0.6.2          cluster_2.1.4          haven_2.5.1           
 [73] fs_1.5.2               magrittr_2.0.3         DO.db_2.9             
 [76] reprex_2.0.2           reactome.db_1.81.0     pcaMethods_1.88.0     
 [79] googledrive_2.0.0      hms_1.1.2              patchwork_1.1.2       
 [82] evaluate_0.16          XML_3.99-0.10          readxl_1.4.1          
 [85] IRanges_2.30.1         gridExtra_2.3          compiler_4.2.1        
 [88] ncdf4_1.19             KernSmooth_2.23-20     crayon_1.5.1          
 [91] shadowtext_0.1.2       htmltools_0.5.3        ggfun_0.0.7           
 [94] tzdb_0.3.0             aplot_0.1.7            lubridate_1.8.0       
 [97] DBI_1.1.3              tweenr_2.0.2           dbplyr_2.2.1          
[100] MASS_7.3-58.1          rappdirs_0.3.3         Matrix_1.5-1          
[103] car_3.1-0              cli_3.3.0              vsn_3.64.0            
[106] parallel_4.2.1         igraph_1.3.5           pkgconfig_2.0.3       
[109] foreach_1.5.2          MALDIquant_1.21        xml2_1.3.3            
[112] ggtree_3.4.2           svglite_2.1.0          webshot_0.5.4         
[115] XVector_0.36.0         rvest_1.0.3            yulab.utils_0.0.5     
[118] digest_0.6.29          graph_1.74.0           Biostrings_2.64.1     
[121] rmarkdown_2.16         cellranger_1.1.0       fastmatch_1.1-3       
[124] tidytree_0.4.1         gtools_3.9.3           graphite_1.42.0       
[127] lifecycle_1.0.2        nlme_3.1-159           jsonlite_1.8.0        
[130] carData_3.0-5          viridisLite_0.4.1      fansi_1.0.3           
[133] pillar_1.8.1           lattice_0.20-45        KEGGREST_1.36.3       
[136] fastmap_1.1.0          httr_1.4.4             GO.db_3.15.0          
[139] glue_1.6.2             iterators_1.0.14       png_0.1-7             
[142] bit_4.0.4              stringi_1.7.8          blob_1.2.3            
[145] caTools_1.18.2         memoise_2.0.1          ape_5.6-2
limmaGUI limma adj.p.value StatisticalMethod • 150 views
ADD COMMENT
1
Entering edit mode
Guido Hooiveld ★ 3.4k
@guido-hooiveld-2020
Last seen 1 day ago
Wageningen University, Wageningen, the …

I did not go through your code, but genes (or proteins) having the same adj.P.Val despite having a different P.Value is a known phenomenon. See e.g. here for more details on this: Limma Voom produces many identical adjusted p-values

ADD COMMENT
0
Entering edit mode

Thanks so much - I was trying to find another post about it but obviously didn't search for the right thing - this is exactly what I wanted!

ADD REPLY

Login before adding your answer.

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