%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%
% \VignetteIndexEntry{EFDR}
%\VignetteDepends{Biobase. golubEsets, ctest, multtest, bioclabs, modreg, wavethresh}
%\VignetteKeywords{Microarray}
\documentclass[12pt]{article}

\usepackage{amsmath,pstricks}
\usepackage[authoryear,round]{natbib}
\usepackage{hyperref}


\textwidth=6.2in
\textheight=8.5in
%\parskip=.3cm
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-.3in

\newcommand{\scscst}{\scriptscriptstyle}
\newcommand{\scst}{\scriptstyle}

\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}


\bibliographystyle{plainnat}


\title{EFDR: Some statistical methods for the identification of differentially
expressed genes}

\begin{document}

\maketitle

\section*{Introduction}

In this lab we introduce you to some extra statistical methods
for the identification of differentially expressed genes between
two populations. This lab is a complement of Lab 4, but relies
upon some ideas that come from wavelet analysis.

<<loadpacks>>=
library(Biobase)
library(golubEsets)
library(ctest)
library(multtest)
library(bioclabs)
library(modreg)
library(wavethresh)
library(Milan)
@

The lab comes in two parts, in the first we consider how wavelet
and thresholding procedures may be used for finding genes that
are differentially expressed,  while in the second part we consider FDR
based procedures that are available in R.

We use expression data from an
experiment where sixteen genes were spiked in at different known 
concentrations in different hybridizations and are thus differentially
expressed. Expression measures are stored in an \Robject{exprSet}
object available through the \Rpackage{bioclabs} package.
<<loaddata>>=
data("eset12")
genenames <- colnames(pData(eset12))[-1]
@

The concentrations of the sixteen genes in each of 24 hybridizations
are stored in the \Robject{phenoData} slot of the \Robject{exprSet}
object \Robject{eset12}.
<<phenodata>>=
pData(eset12)
@

Notice there are two populations and 12 replicates in each. We are
interested in identifying genes that are differentially expressed in
these two populations. 

Let us create a matrix containing for each of the 12626 genes on the
HGU95a chip (note this really is A and not Av2), its {\em A-value} or
average log intensity, its {\em M-value} or difference of log
intensities (log ratio), its two-sample $t$-statistic, and its nominal
$p$-value from the $t$-distribution. The rows of the \Robject{scores}
matrix correspond to genes and the columns to the four different types
of statistics.

The data have already been transformed to the log scale. Base 2
logarithms were used.
%'
<<scores,eval=TRUE>>=
Index1 <- which(eset12$population==0)
Index2 <- which(eset12$population==1)
scores <- esApply(eset12,1,function(x){
  tmp <- t.test(x[Index2],x[Index1],var.equal=TRUE)
  c(mean(tmp$estimate),-diff(tmp$estimate),tmp$statistic,tmp$p.value)
})
scores <- t(scores) ##put genes back in rows
colnames(scores) <- c("A","M","t.stat","p.value")
@


Before further discussing procedures designed to find 
differentially expressed genes, let us illustrate first
the decorrelation properties of the wavelet transform, which
may enhance multiple hypothesis based procedures if they
are applied in the wavelet domain. We select first 
a subseries of length 8192, a power of 2, from the 12626 
M-values. We next display an autocorrelation plot of the 
original $M$-values series, followed by autocorrelation plots
of their wavelet coefficients at several scales. Notice
the decorrelation properties of the discrete wavelet transform.


<<dwt,fig=TRUE>>=
yy<-scores[1:8192,2]
plot(acf(yy))
@

Let us look now the acf of the wavelet coefficients 
of the DWT of the previous series at several scales.

<<dwt2,fig=TRUE>>=
wds1 <- wd(yy)
par(mfcol = c(4, 2),mgp = c(5, 0.4, 0))
plot(acf(accessD(wds1,level=11)))
plot(acf(accessD(wds1,level=10)))
plot(acf(accessD(wds1,level=9)))
plot(acf(accessD(wds1,level=8)))
plot(acf(accessD(wds1,level=7)))
plot(acf(accessD(wds1,level=6)))
plot(acf(accessD(wds1,level=5)))
plot(acf(accessD(wds1,level=4)))
@

The above suggests that wavelet decompostion may enhance FDR multiple hypothesis
testing procedures by coming closer to independance and by reducing efficiently
the number of hypotheses being tested.

\bigskip

We are now going to use the empirical Bayes thresholding
procedure to find differentially expressed genes. There
are  12 replicates and we can use the $t$-test statistics to compare
the gene expressions taking into account their variability.
With 12 replicates but so many genes we expect a large number of false positives and we already
saw in a previous Lab that $t$-tests were not successful in identifying
the sixteen known genes.  The idea is to use the $t$-scores
as noisy observations and a natural approach to this problem, taking into
account the sparsity of significant scores is thresholding: 
if the absolute value of a particular score $Y_i$ exceeds some threshold $t$ then it 
is taken to correspond to a nonzero signal which is then estimated, 
most simply by $Y_i$ itself. If $|Y_i| < t$ then the corresponding gene
expression is estimated to be zero.

<<bayes,eval=TRUE>>=
y<-scores[,3]
fit<- ebayesthresh(y, prior = "cauchy", a = NA, bayesfac = T,sdev = NA, verbose = F, threshrule = "median")
aux<- abs(fit)
auxsv <- aux > 0
cat("Number of signicant genes: \n")
sum(as.integer(auxsv))
rownames(scores[aux>0,])
@


Let's see how we fared with the sixteen known to be differentially
expressed genes.
%'
<<results>>=
genenames
@
%'
The examples in this second part will rely on the data reported in
\citet{Golub99}. The basic comparisons here are between B-cell derived
ALL, T-cell derived ALL and AML. 

What follows serves as a short tutorial on the q-value, as well as instructions for how to
use the q-value functions developed by Storey (2003) and available for use in R.
The q-value is similar to the well known p-value. It gives each hypothesis test a measure
of significance in terms of a certain error rate. The p-value of a test measures the
minimum false positive rate that is incurred when calling that test significant. Likewise,
the q-value of a test measures the minimum false discovery rate that is incurred when
calling that test significant.

We first load the golub data set.

<<loaddatag,fig=TRUE>>=
data(golub) # the golub data set 
classlabel<-golub.cl
dat=golub
o <- rep(T,nrow(dat))
for(i in 1:nrow(dat)) {
if(max(dat[i,])>20) {o[i]<-F}
}
dat <- dat[o,]
n <- ncol(dat)
m <- nrow(dat)
y=golub.cl
y=y+1
#makes two-sample t-statistics
ttest <- function(x,y,f=0) {
 m1 <- apply(x[,y==1],1,mean)
 m2 <- apply(x[,y==2],1,mean)
 v1 <- apply(x[,y==1],1,var)
 v2 <- apply(x[,y==2],1,var)
 s1 <- sqrt(v1/sum(y==1))
 s2 <- sqrt(v2/sum(y==2))
 tt <- (m2-m1)/(sqrt(s1^2 + s2^2)+f)
 return(tt=tt, ss=(sqrt(s1^2 + s2^2)+f))
 }
tt <- ttest(dat,y)$tt
hist(tt,breaks=100,col="red")
@
%'
Let us calculate null statistics and p-values.
<<pvalues,fig=TRUE>>=
tt0 <- 0
set.seed(123)
B <- 100
for(i in 1:B) {
v <- sample(y)
tt0 <- c(tt0,ttest(dat,v)$tt)
}
tt0 <- tt0[-1]
 
#form p-values
att <- abs(tt)
att0 <- abs(tt0)
 v <- c(rep(T,m),rep(F,m*B))
 v <- v[rev(order(c(att,att0)))]
 u <- 1:length(v)
 w <- 1:m
 p <- (u[v==TRUE]-w)/(B*m)
 p <- p[rank(-att)]
 hist(p,breaks=100,col="blue")
@

One can adjust the $p$-values to account for multiple hypothesis
testing using the appropriate functions.  

<<adjustedqvalues>>=
qobj <- qvalue(p)
#qobj$qval contains q-values
#qobj$pval contains p-values
#qobj$pi0 contains pi0 estimate

#calculate fold-change (on log base 2 scale)
m1 <- apply(dat[,y==1],1,mean)
m2 <- apply(dat[,y==2],1,mean)
fc <- (m2-m1)

qval <- rep(NA,length(o))
pval <- rep(NA,length(o))
fchng <- rep(NA,length(o))
qval[o] <- qobj$q
pval[o] <- qobj$pval
fchng[o] <- fc
for(i in 1:length(o)) {
cat(c(round(qval[i],6),pval[i],round(fchng[i],3),'\n'),file="qvalues.txt",append=T)
}
@ 

We now display the figures as they are given in Storey and Tibshirani's paper (2003).
Using a two-sample t-statistic, we calculated a p-value for each of 3051 genes under the null.
The figure 1 below  displays a density histogram of the 3051 p-values from the Golub et al. data. 
The dashed line is the density histogram we would expect if all genes were null 
(not differentially expressed), so it can be seen that many genes are
differentially expressed. The dotted line is at the height of the estimate of the proportion 
of null p-values.


<<Figure1,fig=TRUE>>=
hist(p,nclass=20,main=" ",xlab="Density of observed p-values",ylab=" ",freq=F,ylim=c(0,4))
abline(h=1)
abline(h=0.67,lty=2)
@

The next four figures display  successively (a) The q-values of the genes versus their
respective t-statistics. (b) The q-values versus their respective p-values. (c) The number of genes
occurring on the list up through each q-value versus the respective q-value. (d) The expected number
of false positive genes versus the total number of significant genes given by the q-values.

<<Figure2a,fig=TRUE>>=
plot(tt[order(tt)],qobj$q[order(tt)],type="l",xlim=c(-8,6),xlab="t-statistics",ylab="q-values",main="a")
@

<<Figure2b,fig=TRUE>>=
p2 <- qobj$pval[order(qobj$pval)]
q2 <- qobj$qval[order(qobj$pval)]
plot(p2,q2,type="l",xlab="p-values",ylab="q-values",main="b")
@

<<Figure2c,fig=TRUE>>=
plot(q2,1:m,type="l",xlim=c(0,0.5),ylim=c(0,3051),xlab="q-value",ylab="number of significant genes",main="c")
@

<<Figure2d,fig=TRUE>>=
plot(1:500,q2[1:500]*(1:500),type="l",xlab="number of significant genes",ylab="number of expected false positives",main="d")
@

The next figure displays the estimate $\hat{\pi}_0(\lambda)$ versus $\lambda$
for the Golub data. The solid line is a natural cubic spline
fit to these points to estimate $\hat{\pi}_0(\lambda)$, the proportion 
of non-differentially expressed genes.

<<Figure3,fig=TRUE>>=
lam <- seq(0,0.95,0.01)
pi0 <- rep(0,length(lam))
for(i in 1:length(lam)) {
pi0[i] <- mean(p>lam[i])/(1-lam[i])
}
spi0 <- smooth.spline(lam,pi0,df=3,w=(1-lam))
plot(lam,pi0,xlab=expression(lambda),ylab=expression(hat(pi)[0](lambda)))
lines(spi0)
@

We can easily repeat the whole procedure with the dataset comprising 12
replicates, by simply using \verb+data(esset12)+.

\end{document}
