Source: https://github.com/markziemann/10mistakes
.libPaths("~/.cache/R/")
#detach("package:mitch", unload = TRUE)
#remove.packages("mitch")
install.packages("mitch2_chatgpt", repos = NULL, type="source")## Installing package into '/mnt/data/mdz/.cache/R'
## (as 'lib' is unspecified)library("mitch")
library("kableExtra")
library("getDEE2")
library("DESeq2")## Loading required package: S4Vectors## Loading required package: stats4## Loading required package: BiocGenerics## Loading required package: generics##
## Attaching package: 'generics'## The following objects are masked from 'package:base':
##
## as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
## setequal, union##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique,
## unsplit, which.max, which.min##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## Loading required package: IRanges## Loading required package: GenomicRanges## Loading required package: GenomeInfoDb## Loading required package: SummarizedExperiment## Loading required package: MatrixGenerics## Loading required package: matrixStats##
## Attaching package: 'MatrixGenerics'## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.##
## Attaching package: 'Biobase'## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedianslibrary("tictoc")##
## Attaching package: 'tictoc'## The following object is masked from 'package:SummarizedExperiment':
##
## shift## The following object is masked from 'package:GenomicRanges':
##
## shift## The following object is masked from 'package:IRanges':
##
## shiftlibrary("parallel")
library("plyr")##
## Attaching package: 'plyr'## The following object is masked from 'package:matrixStats':
##
## count## The following object is masked from 'package:IRanges':
##
## desc## The following object is masked from 'package:S4Vectors':
##
## renameRhpcBLASctl::blas_set_num_threads(1)Get gene sets and expression data.
gs <- gmt_import("ReactomePathways.gmt")
if ( ! file.exists("de.Rds") ) {
x <- getDEE2_bundle(species="hsapiens", col="SRP_accession",query="SRP128998", counts="GeneCounts",legacy=TRUE)
names(x)
# set up gene table
gt <- x$GeneInfo
gt <- gt[,1,drop=FALSE]
gt$id <- rownames(gt)
saveRDS(gt,"gt.Rds")
# get sample sheet
ss <- x$MetadataSummary
ss <- ss[grep("VPA",ss$Experiment_title,invert=TRUE),]
ss$case <- grepl("high",ss$Experiment_title)
# get gene counts
x <- x[["GeneCounts"]]
x <- x[,which(colnames(x) %in% rownames(ss))]
dim(x)
head(x)
}Now use deseq2
if ( ! file.exists("de.Rds") ) {
dds <- DESeqDataSetFromMatrix(countData=x, colData = ss, design = ~ case)
dds <- DESeq(dds)
de <- DESeq2::results(dds)
de <- de[order(de$pvalue),]
head(de)
saveRDS(object=de,file="de.Rds")
}First using the recommended code from ChatGPT.
de <- readRDS("de.Rds")
gt <- readRDS("gt.Rds")
library("Matrix")##
## Attaching package: 'Matrix'## The following object is masked from 'package:S4Vectors':
##
## expandm <- mitch_import(x=data.frame(de),DEtype="deseq2",geneTable=gt)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 58302## Note: no. genes in output = 33611## Note: estimated proportion of input genes in output = 0.576m2 <- m
m2$x2<-sample(x=m2$x,size=nrow(m2),replace=FALSE)
tic()
mres <- mitch_calc(x=m2,genesets=gs,minsetsize=5,priority="effect",cores=8)
mresx <- mres$enrichment_result
mresx <- mresx[order(mresx$pMANOVA),]
toc() # old:47s new:48.6(1)## 7.056 sec elapsedhead(mresx) |>
kbl(caption="Top result for ChatGPT script") |>
kable_paper("hover", full_width = F)| set | setSize | pMANOVA | x | x2 | s.dist | SD | p.adjustMANOVA | |
|---|---|---|---|---|---|---|---|---|
| Diseases of Base Excision Repair | Diseases of Base Excision Repair | 5 | 1 | 0.8382670 | -0.3203773 | 0.8974036 | 0.5793221 | 1 |
| NFE2L2 regulating ER-stress associated genes | NFE2L2 regulating ER-stress associated genes | 5 | 1 | -0.8527525 | -0.2812117 | 0.8979236 | 0.2857704 | 1 |
| Unwinding of DNA | Unwinding of DNA | 12 | 1 | 0.8548717 | 0.1206211 | 0.8633395 | 0.3671253 | 1 |
| Lanosterol biosynthesis | Lanosterol biosynthesis | 14 | 1 | 0.8343857 | 0.1682229 | 0.8511748 | 0.3330814 | 1 |
| Cholesterol biosynthesis via desmosterol (Bloch pathway) | Cholesterol biosynthesis via desmosterol (Bloch pathway) | 10 | 1 | 0.8300348 | -0.1819232 | 0.8497375 | 0.5059790 | 1 |
| Activation of C3 and C5 | Activation of C3 and C5 | 7 | 1 | 0.7420630 | -0.3882191 | 0.8374793 | 0.5651411 | 1 |
Now using the existing mitch code.
detach("package:mitch", unload = TRUE)
remove.packages("mitch")## Removing package from '/mnt/data/mdz/.cache/R'
## (as 'lib' is unspecified)BiocManager::install("mitch")## 'getOption("repos")' replaces Bioconductor standard repositories, see
## 'help("repositories", package = "BiocManager")' for details.
## Replacement repositories:
## CRAN: https://cloud.r-project.org## Bioconductor version 3.21 (BiocManager 1.30.27), R 4.5.2 (2025-10-31)## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'mitch'## Installation paths not writeable, unable to update packages
## path: /mnt/data/rlibs/site-library
## packages:
## actuar, ada, altdoc, antiword, arrow, BayesXsrc, BB, broom.mixed, clue,
## connectcreds, countrycode, Cubist, DBI, dm, doFuture, effects, emmeans, ez,
## factoextra, fairml, flashClust, forecast, forestplot, fs, future.apply,
## geosphere, ggiraph, ggpubr, ggrepel, ggstats, graphlayouts, grf, highr,
## ivreg, LaplacesDemon, later, lessR, lhs, libcoin, lme4, lmerTest, lobstr,
## lucid, magick, marginaleffects, mirai, mlbench, mlr3misc, mlr3viz,
## multcompView, nanonext, nleqslv, openssl, openxlsx2, osqp, paws.common,
## pls, PoissonBinomial, R2OpenBUGS, R2WinBUGS, ragg, RcppML, readr, REBayes,
## renv, rgl, rucrdtw, rvg, scatterplot3d, semTools, shiny, shinyWidgets, sn,
## survey, systemfonts, testit, textshaping, tm, tseries, tweedie, vcfR,
## vegan, xtable, xtslibrary("mitch")Original package.
m <- mitch_import(x=data.frame(de),DEtype="deseq2",geneTable=gt)## The input is a single dataframe; one contrast only. Converting
## it to a list for you.## Note: Mean no. genes in input = 58302## Note: no. genes in output = 33611## Note: estimated proportion of input genes in output = 0.576m2 <-m
m2$x2<-sample(x=m2$x,size=nrow(m2),replace=FALSE)
tic()
mres2 <- mitch_calc(x=m2,genesets=gs,minsetsize=5,priority="effect",cores=8)## Note: Enrichments with large effect sizes may not be
## statistically significant.toc() # old:47s new:48.6(1)## 11.607 sec elapsedhead(mres2$enrichment_result) |>
kbl(caption="Top result for original script") |>
kable_paper("hover", full_width = F)| set | setSize | pMANOVA | s.x | s.x2 | p.x | p.x2 | s.dist | SD | p.adjustMANOVA | |
|---|---|---|---|---|---|---|---|---|---|---|
| 2004 | Unwinding of DNA | 12 | 0.0000004 | 0.8548717 | -0.2869108 | 0.0000003 | 0.0850664 | 0.9017336 | 0.8073621 | 0.0000052 |
| 1072 | NFE2L2 regulating ER-stress associated genes | 5 | 0.0030901 | -0.8527525 | -0.2100577 | 0.0009522 | 0.4157126 | 0.8782432 | 0.4544538 | 0.0138686 |
| 134 | Assembly and Release of Dengue Virus Virions | 8 | 0.0001975 | -0.8071898 | 0.2386989 | 0.0000761 | 0.2420767 | 0.8417438 | 0.7395550 | 0.0012594 |
| 462 | Diseases of Base Excision Repair | 5 | 0.0051008 | 0.8382670 | 0.0239541 | 0.0011616 | 0.9260537 | 0.8386092 | 0.5758062 | 0.0208725 |
| 910 | Lanosterol biosynthesis | 14 | 0.0000004 | 0.8343857 | 0.0108853 | 0.0000001 | 0.9437448 | 0.8344567 | 0.5823027 | 0.0000054 |
| 276 | Cholesterol biosynthesis via desmosterol (Bloch pathway) | 10 | 0.0000312 | 0.8300348 | -0.0419214 | 0.0000054 | 0.8183340 | 0.8310928 | 0.6165661 | 0.0002509 |
#detach("package:mitch", unload = TRUE)
#remove.packages("mitch")For reproducibility.
sessionInfo()## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_AU.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_AU.UTF-8 LC_COLLATE=en_AU.UTF-8
## [5] LC_MONETARY=en_AU.UTF-8 LC_MESSAGES=en_AU.UTF-8
## [7] LC_PAPER=en_AU.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_AU.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] mitch_1.20.0 Matrix_1.7-4
## [3] plyr_1.8.9 tictoc_1.2.1
## [5] DESeq2_1.48.2 SummarizedExperiment_1.38.1
## [7] Biobase_2.68.0 MatrixGenerics_1.20.0
## [9] matrixStats_1.5.0 GenomicRanges_1.60.0
## [11] GenomeInfoDb_1.44.3 IRanges_2.42.0
## [13] S4Vectors_0.46.0 BiocGenerics_0.54.1
## [15] generics_0.1.4 getDEE2_1.18.0
## [17] kableExtra_1.4.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 viridisLite_0.4.3 dplyr_1.2.0
## [4] farver_2.1.2 S7_0.2.1 bitops_1.0-9
## [7] fastmap_1.2.0 GGally_2.4.0 promises_1.5.0
## [10] digest_0.6.39 mime_0.13 lifecycle_1.0.5
## [13] magrittr_2.0.4 compiler_4.5.2 rlang_1.1.7
## [16] sass_0.4.10 tools_4.5.2 yaml_2.3.12
## [19] knitr_1.51 S4Arrays_1.8.1 htmlwidgets_1.6.4
## [22] DelayedArray_0.34.1 xml2_1.5.2 RColorBrewer_1.1-3
## [25] BiocParallel_1.42.2 abind_1.4-8 KernSmooth_2.23-26
## [28] purrr_1.2.1 grid_4.5.2 caTools_1.18.3
## [31] xtable_1.8-4 ggplot2_4.0.2 scales_1.4.0
## [34] gtools_3.9.5 MASS_7.3-65 dichromat_2.0-0.1
## [37] cli_3.6.5 rmarkdown_2.30 crayon_1.5.3
## [40] otel_0.2.0 rstudioapi_0.18.0 httr_1.4.8
## [43] reshape2_1.4.5 cachem_1.1.0 stringr_1.6.0
## [46] network_1.20.0 BiocManager_1.30.27 XVector_0.48.0
## [49] vctrs_0.7.1 jsonlite_2.0.0 echarts4r_0.5.0
## [52] beeswarm_0.4.0 systemfonts_1.3.1 locfit_1.5-9.12
## [55] tidyr_1.3.2 jquerylib_0.1.4 glue_1.8.0
## [58] statnet.common_4.13.0 codetools_0.2-20 ggstats_0.12.0
## [61] stringi_1.8.7 gtable_0.3.6 later_1.4.6
## [64] UCSC.utils_1.4.0 htm2txt_2.2.2 tibble_3.3.1
## [67] pillar_1.11.1 htmltools_0.5.9 gplots_3.3.0
## [70] GenomeInfoDbData_1.2.14 R6_2.6.1 textshaping_1.0.4
## [73] evaluate_1.0.5 shiny_1.12.1 lattice_0.22-9
## [76] RhpcBLASctl_0.23-42 httpuv_1.6.16 bslib_0.10.0
## [79] Rcpp_1.1.1 SparseArray_1.8.1 svglite_2.2.2
## [82] coda_0.19-4.1 gridExtra_2.3 xfun_0.56
## [85] pkgconfig_2.0.3