---
title: "BioQC: Detect tissue heterogeneity in gene expression data"
author: "Jitao David Zhang, Klas Hatje, Gregor Sturm, Clemens Broger, Martin Ebeling, and Laura Badi"
date: "`r Sys.Date()`"
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{BioQC: Detect tissue heterogeneity in gene expression data}
  %\usepackage[utf8]{inputenc}
output:
  rmarkdown::html_vignette:
    self_contained: no
  md_document:
    variant: markdown_phpextra
    preserve_yaml: TRUE 
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(out.extra='style="display:block; margin: auto"', fig.align="center")
```

Introduction
------------

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.



A dummy example
---------------
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.

```{r seed, include=FALSE}
set.seed(1887)
```

```{r lib, warning=FALSE, message=FALSE, results="hide"}
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.

```{r gmtReadin}
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.

```{r runBioQC}
dummyRes <- wmwTest(myEset, gmt, valType="p.greater", simplify=TRUE)
summary(p.adjust(dummyRes, "BH"))
```



### Using basic data structures
The 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.

```{r basicDs}
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)
```


Case study with real dataset
----------------------------
We 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. 

The data and the script used to perform the analysis can be found in the github repository [BioQC-example](https://github.com/Accio/BioQC-example). 



Benchmarking against R implementation
-------------------------------------
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.

```{r benchmark}
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))
system.time(Rres <- apply(exprs, 2, function(x)
                          sapply(ind, function(y)
                                 wmwTestInR(x, y, valType="p.less"))))

```

With 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.

Technical notes
----------------

### Background matters

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.

```{r backgroundMatters}
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))
(pHighExpGenes <- wmwTest(bgVec[1:10000], ind))
```

It is therefore recommended to report any filtering prior to the BioQC analysis to make the results reproducible.

Acknowledgement
----------------
We would like to thank Iakov Davydov for suggestions to improve the package.

R Session Info
----------------
```{r session_info}
sessionInfo()
```
