%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%
% \VignetteIndexEntry{Lab 4}
%\VignetteDepends{Biobase,ctest,multtest,bioclabs}
%\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{Lab 4: Differential Gene Expression}

\begin{document}

\maketitle

In this lab, we demonstrate how to use R to find genes that are
differentially expressed in two populations. We demonstrate two useful
plots, the MA-plot and the volcano plot. For a more formal assessment,
we use the \verb+multtest+ package for obtaining adjusted $p$-values.
<<loadpacks>>=
library(Biobase)
library(ctest)
library(multtest)
library(bioclabs)
@

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("eset3")
genenames <- colnames(pData(eset3))[-1]
@
%%FIXME: couldn't we have better varLabels?
The concentrations of the sixteen genes in each of six hybridizations
are stored in the \Robject{phenoData} slot of the \Robject{exprSet}
object \Robject{eset3}.
<<phenodata>>=
pData(eset3)
@
Notice there are two populations and 3 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(eset3$population==0)
Index2 <- which(eset3$population==1)
scores <- esApply(eset3,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")
@

The following commands produce an {\em MA-plot} of the differences of
log intensities in the two populations, M, vs. the average log
intensities, A.
<<MvA,fig=TRUE>>=
plot(scores[,1],scores[,2],xlab="A",ylab="M",pch=".")
#points(scores[genenames,1],scores[genenames,2],pch=16,col=1:16)
text(scores[genenames,1],scores[genenames,2],genenames,col=1:16)
abline(h=c(-1,1)) ##fc = 2
@

In the MA-plot, points with large vertical deviations (absolute
M-values) suggest differentially expressed genes.  The horizontal
lines show the typical two-fold-change cutoff (because the expression
data are on a log 2 scale).  The colored symbols correspond to the
sixteen genes that are truly differentially expressed, i.e., were
spiked in at different concentrations. Other points with large
M-values correspond to false positives. One of the sixteen known genes
has a low M-value, corresponding to a false negative ({\tt 1597\_at}).
Based on our knowledge of the concentrations for each of the sixteen
genes, one expects only one of these genes ({\tt 684\_at}) to have a
negative M-value, i.e., higher expression measure in population 0 than
1 (infinitely more abundant in population 0). This is indeed the point
with the very low M-value.


Should we take the variability of the estimates into account? There
are only 3 replicates but we can try a $t$-test. The following is a
so-called {\em volcano plot} of the $t$-statistic vs. the numerator of
the $t$-statistic, or M-value. Genes in the top left and right corners
of the plot correspond to genes with both large absolute differences
and large relative (to standard error) differences in expression
between the two populations. Although the sixteen known genes tend to
have large absolute $t$-statistics, they standout more in terms of their
M-values.

In other situations you may have seen the volcano plot defined as a
plot of minus the base 10 logarithm of the data versus fold
change. The two versions (ours and this one) are equivalent since all
$t$-statistics are based on the same number of degrees of freedom.
<<volcano,fig=TRUE>>=
plot(scores[,2],abs(scores[,3]),xlab="M",ylab="t.stat",pch=".")
points(scores[genenames,2],abs(scores[genenames,3]),pch=16,col=1:16)
abline(v=c(-1,1))
a <- qt(.975,4)
abline(h=a)
@

Note that the gene with small M value really distorts our visual
impression of the data. In order to remove that effect we replot the
data with new limits.  For a close-up, simply change the \verb+xlim+
and \verb+ylim+.
<<volcanocloseup,fig=TRUE>>=
plot(scores[,2],abs(scores[,3]),xlab="M",
       ylab="t.stat",pch=".",xlim=c(-1.5,1.5),ylim=c(0,15))
text(scores[genenames,2],abs(scores[genenames,3]),genenames,col=1:16)
abline(v=c(-1,1)); a <- qt(.975,4); abline(h=a)
@

How many genes have $p$-values less than 0.05? How about 0.01?
<<pvalues>>=
sum(scores[,"p.value"]<=0.05)
sum(scores[,"p.value"]<=0.01)
@

One can adjust the $p$-values to account for multiple hypothesis
testing using the 
\Rpackage{multtest} package.  The function
\Robject{mt.rawp2adjp} gives adjusted $p$-values according to various
methods using only the raw $p$-values.
<<adjustedpvalues>>=
tmp <- mt.rawp2adjp(scores[,4])
adj.p.values <- tmp$adjp[order(tmp$index),-1]
scores <- cbind(scores,adj.p.values)
@ 
The function \Rfunction{maxT} computes permutation adjusted
$p$-values for the Westfall \& Young step-down maxT procedure. One
needs to use the original data again. Similarly, for the step-down
minP procedure, one can use the function \verb+mt.minP+.
<<maxT>>=
tmp <- mt.maxT(exprs(eset3),eset3$population)
scores <- cbind(scores,tmp$adjp[order(tmp$index)])
colnames(scores)[12] <- "maxT"
@
%$ -- to keep fontlocking happier

Let's see how we fared with the sixteen known to be differentially
expressed genes.
%'
<<results>>=
round(scores[genenames,],2)
@

One can make some pretty substantial arguments against this
procedure. Make an adjustment for all 12626 probes is probably not
appropriate in any situation and is certainly hard to support
here. These data arose from a designed experiment where we know which
genes are going to change. Not all 12626 were candidates for
change. In a real experiment you are likely to find that approximately
40\% of the genes are not expressed in the tissue that you are
studying. In that case you have corrected for something which should
have been excluded from the analysis. There is of course an issue in
determining which genes are not differentially expressed is very
difficult but the gains from an approximate solution are likely to be
quite large.

With 3 replicates we don't expect to have much power. Should we even
use $t$-statistics over the more simple fold change estimates?  Let's
see which one does better at ranking the sixteen truly differentially
expressed genes. Ideally, one would like the sixteen genes to have
ranks 1 through 16. It seems like the simple M-value or fold change
measure was more successful at identifying the sixteen known genes.
%'
<<ranks>>=
m.ranks <- rank(-abs(scores[,2])) ##by fold change
names(m.ranks) <- rownames(scores)
t.ranks <- rank(-abs(scores[,3])) ##by fold change
names(t.ranks) <- rownames(scores)
cbind(m.ranks,t.ranks)[genenames,]
@

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

\end{document}







