BioQC Algorithm: Speeding up the Wilcoxon-Mann-Whitney Test

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:

1. 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.
2. The core algorithm is implemented in C instead of R as programming language.
3. 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.

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 63891 human 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.

expect_equal(bioqcNumRes, rNumRes)

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:

Figure 2: Time benchmark results of BioQC and R implementation of Wilcoxon-Mann-Whitney test. Left panel: elapsed time in seconds (logarithmic Y-axis). Right panel: ratio of elapsed time by two implementations. All results achieved by a single thread on in a RedHat Linux server.

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

sessionInfo()

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows Server x64 (build 17763)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                          
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] rbenchmark_1.0.0     gplots_3.1.1         gridExtra_2.3       
##  [4] latticeExtra_0.6-29  lattice_0.20-45      org.Hs.eg.db_3.14.0 
##  [7] AnnotationDbi_1.56.0 IRanges_2.28.0       S4Vectors_0.32.0    
## [10] testthat_3.1.0       limma_3.50.0         RColorBrewer_1.1-2  
## [13] BioQC_1.22.0         Biobase_2.54.0       BiocGenerics_0.40.0 
## [16] knitr_1.36          
## 
## loaded via a namespace (and not attached):
##  [1] KEGGREST_1.34.0        gtools_3.9.2           locfit_1.5-9.4        
##  [4] xfun_0.27              bslib_0.3.1            vctrs_0.3.8           
##  [7] htmltools_0.5.2        yaml_2.2.1             blob_1.2.2            
## [10] rlang_0.4.12           jquerylib_0.1.4        withr_2.4.2           
## [13] DBI_1.1.1              bit64_4.0.5            jpeg_0.1-9            
## [16] GenomeInfoDbData_1.2.7 stringr_1.4.0          zlibbioc_1.40.0       
## [19] Biostrings_2.62.0      gtable_0.3.0           caTools_1.18.2        
## [22] evaluate_0.14          memoise_2.0.0          fastmap_1.1.0         
## [25] GenomeInfoDb_1.30.0    highr_0.9              Rcpp_1.0.7            
## [28] KernSmooth_2.23-20     edgeR_3.36.0           cachem_1.0.6          
## [31] desc_1.4.0             pkgload_1.2.3          jsonlite_1.7.2        
## [34] XVector_0.34.0         bit_4.0.4              png_0.1-7             
## [37] digest_0.6.28          stringi_1.7.5          rprojroot_2.0.2       
## [40] grid_4.1.1             tools_4.1.1            bitops_1.0-7          
## [43] magrittr_2.0.1         sass_0.4.0             RCurl_1.98-1.5        
## [46] RSQLite_2.2.8          crayon_1.4.1           pkgconfig_2.0.3       
## [49] rmarkdown_2.11         httr_1.4.2             R6_2.5.1              
## [52] compiler_4.1.1