%\VignetteIndexEntry{BioQC: The kidney expression example}
%\VignetteDepends{Biobase}
%\VignettePackage{BioQC}

\documentclass[11pt]{article}

\usepackage{times}
\usepackage{hyperref}
\usepackage{geometry}
\usepackage{longtable}
\usepackage[pdftex]{graphicx}
\SweaveOpts{keep.source=TRUE,eps=FALSE,pdf=TRUE,prefix=TRUE,eval=TRUE} 

% R part
\newcommand{\R}[1]{{\textsf{#1}}}
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}
\newcommand{\Metas}[1]{{\texttt{#1}}}

\begin{document}
\setkeys{Gin}{width=0.8\textwidth}
\title{Case studies of BioQC}
\author{Jitao David Zhang}

\maketitle

\begin{abstract}
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 \emph{C} routine to allow fast computation in scenarios where a numeric vector (\emph{e.g.} expression values of all genes) is compared against a collection of subsets of the vector (\emph{e.g.} expression values of genes belonging to a collection of gene sets). We demonstrate the use of the package with the \emph{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.
\end{abstract}

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 \textit{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 by follow-up experiments.
  
\section{A dummy example}

We demonstrate the basic use of the package by a dummy example. First we load BioQC library and the tissue signatures into the R session. 

<<setSeed, echo=FALSE>>=
set.seed(1887)
@ 

<<libLoad>>=
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.
<<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.
<<runBioQC>>=
dummyRes <- wmwTest(myEset, gmt, alternative="greater")
summary(p.adjust(dummyRes, "BH"))
@ 

\subsection{Using basic data structures}
The dummy example above shows how to run BioQC algorithm with an \textit{ExpressionSet} and a list read from GMT files (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 to interested user to the documentation of \textit{wmwTest} function.

<<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)
@ 

\section{Case study with a real data set}
As another example, we have applied BioQC to a real gene expression profiling dataset. BioQC helped to generate hypotheses about potential contamination of rat kidney examples by pancreas tissues, which was confirmed by further qRT-PCR experiments. 

The data and script used to perform the analysis can be found at \url{https://github.com/Accio/BioQC-example}.They are not included in the package due to size limitations.

\section{Benchmarking against R implementation}
In the core of \textit{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 \textit{stats} (\textit{wilcox.text}) and \textit{limma} (\textit{rankSumTestWithCorrelation}) packages, the \textit{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. 

<<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, alternative="less"))
wmwTestR <- function(matrix, index, alternative, stat) {
  sub <- rep(FALSE, length(matrix))
  sub[index]=TRUE
  return(wilcox.test(matrix[sub], matrix[!sub], alternative=alternative)$p.value)
}
system.time(Rres <- apply(exprs, 2, function(x)
                          sapply(ind, function(y) 
                                 wmwTestR(x, y, alternative="less", stat=FALSE))))

@ 

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.

\section{Session Info}
The script runs within the following session:
<<sessionInfo, echo=FALSE, results=verbatim>>=
sessionInfo()
@ 

\end{document}