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