%\VignetteDepends{ALL, genefilter, hgu95av2.db, limma}
%\VignetteIndexEntry{Microarray Analysis: Self-Study Exercise}
%\VignettePackage{SeattleIntro2010}

\documentclass{article}

\usepackage{BioconductorVignette}
\SweaveOpts{keep.source=TRUE,width=9,height=8,eps=FALSE,include=FALSE}

\title{Microarray Analysis: Self-Study Exercises}
\author{Chao-Jen Wong \\ Fred Hutchinson Cancer Research Center}
\date{9-10 December, 2010}

\begin{document}
\maketitle

<<setup, echo=FALSE>>=
options(width = 40)
@ 


%% keep this code chuck here in case the package won't build!
<<get the dataset, echo=FALSE>>=
library(ALL) ## data package
library(hgu95av2.db) ## annotation package
data(ALL)
## subsetting
bcell <- grep("^B", as.character(ALL$BT))
types <- c("NEG", "BCR/ABL")
moltyp <- which(as.character(ALL$mol.biol) %in% types)
ALL_bcrneg <- ALL[, intersect(bcell, moltyp)]

ALL_bcrneg$mol.biol <- factor(ALL_bcrneg$mol.biol)

## nonspecific filtering
library(genefilter)
filt_bcrneg <- nsFilter(ALL_bcrneg, require.entrez=TRUE,
                        remove.dupEntrez=TRUE,
                        var.func=IQR,
                        var.cutoff=0.5, feature.exclude="^AFFX")
## get the real dataset
ALLfilt_bcrneg <- filt_bcrneg$eset
@ 


For this lab activity we use a subset of \Robject{ALL} called
\Robject{ALLfilt\_bcrneg} for differential analysis. It  was
constructed to obtained samples from B-cell 
tumors harboring the BCR/ABL mutation and from B-cell tumors with no observed
cytogenetic abnormalities (NEG). The goal is to find genes that are
differentially expressed between these two distinct molecular subtypes of B-cell
leukemia. 

\section{Loading data}
\Robject{ALLfilt\_bcrneg} is an instance of \Robject{ExpressionSet} and is
available in the package \Rpackage{SeattleIntro2010}. First of all, we load
the dataset to the workspace:


<<load dataset, eval=FALSE>>=
library(SeattleIntro2010)
data(ALLfilt_bcrneg)
@ 
\section{$t$-tests}
We can now perform gene-by-gene tests for differential expression using a simple
two-class $t$-tests. The function \Rfunction{rowttest} in the
\Rpackage{genefilter} package uses the $t$-test, row by row, to detect
significant differences in the location of the distribution of expression values
of two groups of samples defined by a factor variable.

\begin{Exercise}
  Use \Rfunction{rowttest} to perform $t$-tests for each genes.
\end{Exercise}

\begin{Solution}
<<t-test>>=
library(genefilter)
table(ALLfilt_bcrneg$mol.biol)
tt <- rowttests(ALLfilt_bcrneg, "mol.biol")
@ 
\end{Solution}

\begin{Exercise}
  Consult the manual for \Rfunction{rowttest} for the meaning of the four
  elements of the return value \Robject{tt}.
\end{Exercise}

\begin{Exercise}
  Create a histogram of the resulting $p$-values.
\end{Exercise}

\section{Multiple testing correction}
The \Rpackage{multtest} package (by K. Pollard, Y. Ge and S. Dudoit) provides a
wide variety of $p$-value correction methods for controlling a broad class of
Type I error rates. These methods including single-step procedures,
FWER-controlling and  FDR-controlling procedures. The results of $p$-value
correction can be more informative for choosing selection cut-offs. Here, we
are going to adjust the $p$-value using the procedure of Benjamini and
Hochberg for controlling the false discovery rate (FDR). Note that FDR is the
expected proportion of false positives among the genes that are called
differentially expressed. 

\begin{Exercise}
  Use the function \Rfunction{mt.rawp2adjp} and  Benjamini and Hochberg
  procedure to adjust the $p$-value. How many genes have adjusted $p$-values
  less than $0.05$? Print out the ten highest-ranking genes with respect to
  the adjusted $p$-value.
\end{Exercise}

\begin{Solution}
<<get multtest>>=
library(multtest)
mt <- mt.rawp2adjp(tt$p.value, proc="BH")
head(mt$adjp) ## print out the first 6 genes with their raw p and adjusted p
sum(mt$adjp[, "BH"] < 0.05)
mt$index[1:10]
featureNames(ALLfilt_bcrneg)[mt$index[1:10]]
@ 
\end{Solution}

\section{Linear model and moderate $t$-tests}
\Rpackage{limma} is a Bioconductor package that provides advanced statistical
methods for linear modeling of microarray data and for identifying 
differentially expressed genes. Our goal here is to get familiar with the
steps of fitting a linear model to the data and using moderate $t$-tests for
assessing differential expression.

\subsection{Example 1}
First, define the design matrix.
<<design matrix>>=
library(limma)
design <- model.matrix( ~mol.biol, ALLfilt_bcrneg)
head(design)
@ 
The above code chunk creates an $n \times 2$ binary matrix, where $n$ is the size of the
samples. This matrix is passed to the \Rfunction{lmFit} function to fit
linear model $y_i = \mu + \beta a_{ij}$ to the expression data. Note that
$a_{ij}=1$ if $i \in \{BCR/ABL\}$, $\mu$ is the mean expression of BCR/ABL
gene and $\beta$ represents the effect of NEG on expression level of gene $i$.
Next we pass this linear model to the function \Rfunction{eBayes} to calculate
the moderated $t$-statistics corresponding to the coefficients of the linear
model, i.e., $\mu$ and $\beta$.

<<linear models and ebayes>>=
fit1 <- lmFit(exprs(ALLfilt_bcrneg), design)
fit2 <- eBayes(fit1)
topTable(fit2, coef=2, adjust.method="BH", 
         number=5)
@     

\begin{Question}
  Why is the input argument \Robject{coef} of \Rfunction{topTable} set to $2$?
\end{Question}

\begin{Exercise}
  Let the cut-off of the adjusted $p$-value be $0.05$. How many probe sets are
  corresponding to differentially expressed genes?
\end{Exercise}
  
\begin{Solution}
<<expressed genes>>=
result <- topTable(fit2, coef=2, adjust.method="BH", 
                   p.value=0.05, number=nrow(fit2))
nrow(result)       

@ 
\end{Solution}

\subsection{Example 2}
Create another design matrix in the following way:
<<design matrix 2>>=
design <- model.matrix( ~0+mol.biol, ALLfilt_bcrneg)
colnames(design) <- c("BCR_ABL", "NEG")
head(design)
@ 
The linear model is then constructed in the following expression:
 \[ y_i = \beta_1 a_{ij} + \beta_2 b_{ij} + \varepsilon_i, \]
where $\beta_1$ is the mean expression of BCR/ABL and $\beta_2$ of NEG. In
this case, we need to let \Rfunction{eBayes} know which contrast
estimator to use. If we are interested in the effect of
BCR/ABL on the expression level of genes, relative to NEG, we can construct a
contrast matrix by the code below:
<<contrast matrix>>=
## note that BCR_ABL and NEG are the column names of design matrix
contr <- makeContrasts(BCR_ABL-NEG, levels=design)
@ 
<<simplier, eval=FALSE>>=
## or simply 
contr <- c(1, -1)
@ 
Now we are ready to perform further differential analysis.
<<example 2>>=
fit <- lmFit(exprs(ALLfilt_bcrneg), design)
fit1 <- contrasts.fit(fit, contr)
fit2 <- eBayes(fit1)
@ 
\begin{Exercise}
  Use the function \Rfunction{topTable} to print out the ten highest
  ranking genes, which are differentially expressed in the BCR/ABL samples
  relative to the NEG samples.
\end{Exercise}

\begin{Solution}
<<print out top 10 gene>>=
topTable(fit2,  adjust.method="BH",  number=10)
@ 
\end{Solution}


\end{document}
