BioQC Algorithm: Speeding up the Wilcoxon-Mann-Whitney Test
Gregor Sturm and Jitao David Zhang
2018-10-30
Supplementary Information for “Detect issue heterogenity in gene expression data with BioQC” (Jitao David Zhang, Klas Hatje, Gregor Sturm, Clemens Broger, Martin Ebeling, Martine Burtin, Fabiola Terzi, Silvia Ines Pomposiello and Laura Badi)
In this vignette, we explain the underlaying algorithmic details of our implementation of the Wilcoxon-Mann-Whitney test. The source code used to produce this document can be found in the github repository BioQC.
BioQC is a R/Bioconductor package to detect tissue heterogeneity from high-throughput gene expression profiling data. It implements an efficient Wilcoxon-Mann-Whitney test, and offers tissue-specific gene signatures that are ready to use ‘out of the box’.
Algorithmic improvements
The Wilcoxon-Mann-Whitney (WMW) test is a non-parametric statistical test to test if median values of two population are equal or not. Unlike the t-test, it does not require the assumption of normal distributions, which makes it more robust against noise.
We improved the computational efficiency of the Wilcoxon-Mann-Whitney test in comparison to the native R implementation based on three modifications:
- The approximative WMW-statistic (Zar, J. H. (1999). Biostatistical analysis. Pearson Education India. pp. 174-177) is used. The differences to the exact version are negligible for high-throughput gene expression data.
- The core algorithm is implemented in C instead of R as programming language.
- BioQC avoids futile expensive sorting operations.
While (1) and (2) are straightforward, we elaborate (3) in the following.
Let \(W_{a,b}\) be the approximative WMW test of two gene vectors \(a,b\), where \(a\) is the gene set of interest, typically containing less than a few hundreds of genes, and \(b\) is the set of all genes outside the gene set (background genes) typically containing \(>10000\) genes. In the context of BioQC, the gene sets are referred to as tissue signatures.
Given an \(m \times n\) input matrix of gene expression data with \(m\) genes and \(n\) samples \(s_1, \dots, s_n\), and \(k\) gene sets \(d_1, \dots, d_k\), the WMW-test needs to be applied for each sample \(s_i, i \in 1..n\) and each gene set \(d_j, j \in 1..k\). The runtime of the WMW-test is essentially determined by the sorting operation on the two input vectors. Using native R wilcox.test, the vectors \(a\) and \(b\) are sorted individually for each gene set. However, in the context of gene set analysis, this is futile, as the (large) background set changes insignificantly in relation to the (small) gene set, when testing different gene sets on the same sample.
Therefore, we approximate the WMW-test by extending \(b\) to all genes in the sample, keeping the background unchanged when testing multiple gene sets. Like this, \(b\) has to be sorted only once per sample. The individual gene sets still need to be sorted, which is not a major issue, as they are small in comparison to the set of background genes.
Figure 1: BioQC speeds up the Wilcoxon-Mann-Whitney test by avoiding futile sorting operations on the same sample.
Time benchmark
To demonstrate BioQC’s superior performance, we apply both BioQC and the native R wilcox.test to random expression matrices and measure the runtime.
We setup random expression matrices of 20186 human protein-coding genes of 1, 5, 10, 50, or 100 samples. Genes are \(i.i.d\) distributed following \(\mathcal{N}(0,1)\). The native R and the BioQC implementations of the Wilcoxon-Mann-Whitney test are applied to the matrices respectively.
The numeric results of both implementations, bioqcNumRes (from BioQC) and rNumRes (from R), are equivalent, as shown by the next command.
The BioQC implementation is more than 500 times faster: while it takes about one second for BioQC to calculate enrichment scores of all 155 signatures in 100 samples, the native R implementation takes about 20 minutes:
Conclusion
We have shown that BioQC achieves identical results as the native implementation in two orders of magnitude less time. This renders BioQC a highly efficient tool for quality control of large-scale high-throughput gene expression data.
R Session Info
## R version 3.5.1 Patched (2018-07-12 r74967)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.8-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.8-bioc/R/lib/libRlapack.so
##
## 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
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] rbenchmark_1.0.0 gplots_3.0.1 gridExtra_2.3
## [4] latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-35
## [7] hgu133plus2.db_3.2.3 org.Hs.eg.db_3.7.0 AnnotationDbi_1.44.0
## [10] IRanges_2.16.0 S4Vectors_0.20.0 BioQC_1.10.0
## [13] Biobase_2.42.0 BiocGenerics_0.28.0 Rcpp_0.12.19
## [16] testthat_2.0.1 knitr_1.20
##
## loaded via a namespace (and not attached):
## [1] highr_0.7 compiler_3.5.1 bitops_1.0-6
## [4] tools_3.5.1 digest_0.6.18 bit_1.1-14
## [7] RSQLite_2.1.1 evaluate_0.12 memoise_1.1.0
## [10] gtable_0.2.0 pkgconfig_2.0.2 rlang_0.3.0.1
## [13] DBI_1.0.0 yaml_2.2.0 stringr_1.3.1
## [16] caTools_1.17.1.1 gtools_3.8.1 rprojroot_1.3-2
## [19] bit64_0.9-7 grid_3.5.1 R6_2.3.0
## [22] rmarkdown_1.10 gdata_2.18.0 blob_1.1.1
## [25] magrittr_1.5 backports_1.1.2 htmltools_0.3.6
## [28] KernSmooth_2.23-15 stringi_1.2.4
