The Wilcoxon-Mann-Whitney rank sum test is one of the most commonly used algorithms for non-parametric tests. The current implementations in R however are not optimal for multiple hypothesis testing. The BioQC package implements the algorithm in a native C routine to allow fast computation in scenarios where a numeric vector (e.g. expression values of all genes) is compared against a collection of subsets of the vector (e.g. expression values of genes belonging to a collection of gene sets). We demonstrate the use of the package with the BioQC algorithm, which uses the Wilcoxon-Mann-Whitney rank sum test and a collection of tissue-preferentially gene signatures to detect tissue heterogeneity in high-throughput gene expression studies.
In this vignette we demonstrate the use of BioQC by a simulated expression dataset, and compare the performance of Wilcoxon-Mann-Whitney rank sum test algorithm implemented in BioQC to other implementations available in R. To use BioQC, the users only need to provide an expression dataset, in the form of a numeric matrix, or an ExpressionSet object. The BioQC package provides tissue-specific genes that can be used directly with the algorithm. The output is one score for each tissue type and each sample. The ranks of the score within each sample can be compared with prior knowledge about the sample to infer tissue heterogeneity. The hypotheses generated by BioQC should be further tested in follow-up experiments.
We demonstrate the basic use of the package with a dummy example. First, we load BioQC library and the tissue signatures into the R session.
library(Biobase)
library(BioQC)
gmtFile <- system.file("extdata/exp.tissuemark.affy.roche.symbols.gmt", package="BioQC")
gmt <- readGmt(gmtFile)Next, we synthesize an ExpressionSet object, using randomly generated data.
Nrow <- 2000L
Nsample <- 5L
gss <- unique(unlist(sapply(gmt, function(x) x$genes)))
myEset <- new("ExpressionSet",
exprs=matrix(rnorm(Nrow*Nsample), nrow=Nrow),
featureData=new("AnnotatedDataFrame",
data.frame(GeneSymbol=sample(gss, Nrow))))Finally we run the BioQC algorithm and print the summary of the results. As expected, no single tissue scored significantly after multiple correction.
dummyRes <- wmwTest(myEset, gmt, valType="p.greater", simplify=TRUE)
summary(p.adjust(dummyRes, "BH"))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.9215 0.9347 1.0000 0.9743 1.0000 1.0000The dummy example above shows how to run BioQC algorithm with an ExpressionSet and a list read from a GMT file (a file format capturing gene sets). Users can also use more basic data structures (e.g. matrices and list of integer indexes) to run the algorithm as shown by the following example. Other data structures will be coerced into these basic data structures; we refer the interested user to the documentation of the wmwTest function.
myMatrix <- matrix(rnorm(Nrow*Nsample),
ncol=Nsample,
dimnames=list(NULL, LETTERS[1:Nsample]))
myList <- list(signature1=sample(1:Nrow, 100),
signature2=sample(1:Nrow, 50),
signature3=sample(1:Nrow, 200))
wmwTest(myMatrix, myList, valType="p.greater", simplify=TRUE)## A B C D E
## signature1 0.1529726 0.39838925 0.44942433 0.0827896 0.8544985
## signature2 0.1521229 0.63237801 0.01732592 0.5765494 0.3946468
## signature3 0.8593754 0.07167512 0.16760316 0.0523685 0.3693992We have applied BioQC to a real gene expression profiling dataset. BioQC helped to generated hypotheses about potential contamination of rat kidney examples by pancreas tissues, which could be confirmed with qRT-PCR.
See the BioQC-kidney vignette for more details.
In the core of wmwTest, an efficient C-implementation makes it feasible to run a large number of Wilcoxon-Mann-Whitney tests on large-scale expression profiling datasets and with many signature lists. Compared to native R implementations in stats (wilcox.text) and limma (rankSumTestWithCorrelation) packages, the BioQC implementation requires less memory and avoids repetitive statistical ranking of data.
The following code, though in a small scale, demonstrates the difference between the performances of two implementations.
bm.Nrow <- 22000
bw.Nsample <- 5
bm.Ngs <- 5
bm.Ngssize <- sapply(1:bm.Ngs, function(x) sample(1:bm.Nrow/2, replace=TRUE))
ind <- lapply(1:bm.Ngs, function(i) sample(1:bm.Nrow, bm.Ngssize[i]))
exprs <- matrix(round(rnorm(bm.Nrow*bw.Nsample),4), nrow=bm.Nrow)
system.time(Cres <- wmwTest(exprs, ind, valType="p.less", simplify=TRUE))## user system elapsed
## 0.067 0.000 0.067system.time(Rres <- apply(exprs, 2, function(x)
sapply(ind, function(y)
wmwTestInR(x, y, valType="p.less"))))## user system elapsed
## 0.939 0.000 0.939With 22000 genes, five samples, and five gene sets, the BioQC implementation is about 20x faster than the R implementation (dependent on individual machines and settings). Our benchmark shows that with the same number of genes, 2000 samples and 200 gene sets (similar to the total number of tissues collected in the BioQC signature list), the BioQC implementation can be about 1000x faster than the R implementation.
Even using the same set of signatures, the p-value of enrichment may vary depending on the background. A toy example is shown below: assuming that in a gene expression profile with values of 20000 genes, there are half of them lower expressed than the other half, mimicking the commonly observed patterns that some genes are almost not or very low expressed; depending on whether including these lowly expressed genes or not, the p-value reported by the Wilcoxon’s test varies.
bgVec <- rnorm(20000)
bgVec[1:10000] <- bgVec[1:10000] + 2
bgVec[1:10] <- bgVec[1:10] + 1
ind <- c(1:10)
(pAllGenes <- wmwTest(bgVec, ind))## [1] 4.2681e-06## [1] 0.0002268254It is therefore recommended to report any filtering prior to the BioQC analysis to make the results reproducible.
We would like to thank Iakov Davydov for suggestions to improve the package.
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 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_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] rbenchmark_1.0.0 gplots_3.2.0 gridExtra_2.3
## [4] latticeExtra_0.6-30 lattice_0.22-6 org.Hs.eg.db_3.20.0
## [7] AnnotationDbi_1.69.0 IRanges_2.39.2 S4Vectors_0.43.2
## [10] testthat_3.2.1.1 limma_3.61.12 RColorBrewer_1.1-3
## [13] BioQC_1.35.0 Biobase_2.65.1 BiocGenerics_0.51.3
## [16] knitr_1.48 rmarkdown_2.28
##
## loaded via a namespace (and not attached):
## [1] KEGGREST_1.45.1 gtable_0.3.6 xfun_0.48
## [4] bslib_0.8.0 caTools_1.18.3 bitops_1.0-9
## [7] vctrs_0.6.5 tools_4.4.1 RSQLite_2.3.7
## [10] highr_0.11 blob_1.2.4 pkgconfig_2.0.3
## [13] KernSmooth_2.23-24 desc_1.4.3 lifecycle_1.0.4
## [16] GenomeInfoDbData_1.2.13 compiler_4.4.1 deldir_2.0-4
## [19] Biostrings_2.73.2 brio_1.1.5 statmod_1.5.0
## [22] GenomeInfoDb_1.41.2 htmltools_0.5.8.1 sys_3.4.3
## [25] buildtools_1.0.0 sass_0.4.9 yaml_2.3.10
## [28] crayon_1.5.3 jquerylib_0.1.4 cachem_1.1.0
## [31] gtools_3.9.5 locfit_1.5-9.10 digest_0.6.37
## [34] maketools_1.3.1 rprojroot_2.0.4 fastmap_1.2.0
## [37] grid_4.4.1 cli_3.6.3 magrittr_2.0.3
## [40] edgeR_4.3.21 UCSC.utils_1.1.0 bit64_4.5.2
## [43] XVector_0.45.0 httr_1.4.7 jpeg_0.1-10
## [46] bit_4.5.0 interp_1.1-6 png_0.1-8
## [49] memoise_2.0.1 evaluate_1.0.1 rlang_1.1.4
## [52] Rcpp_1.0.13 glue_1.8.0 DBI_1.2.3
## [55] pkgload_1.4.0 jsonlite_1.8.9 R6_2.5.1
## [58] zlibbioc_1.51.2