%\VignetteDepends{ALL, affy, arrayQualityMetrics, genefilter, lattice, limma}
%\VignetteIndexEntry{Microarray Analysis}
%\VignettePackage{SeattleIntro2010}

\documentclass{beamer}

\usepackage{BioconductorSlides}
\SweaveOpts{keep.source=TRUE}

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

\begin{document}

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

\frame{\titlepage}

\section[Outline]{}
\frame{\tableofcontents}

\section{Basics}

\begin{frame}
  \frametitle{Introduction}
  \begin{itemize}
    \item Identify differentially expressed genes associated with biological
      or experimental conditions.
    \item Primarily concerned with two-class problems.
    \item Data with $n$ samples and $p$ probes ($p >> n$).
      \begin{table}
        \begin{tabular}{ccccc|ccccc}
          A & A & A & A & A &  B & B & B & B & B \\ \hline
          $x_{1,1}$ & $x_{1,2}$ & $x_{1,3}$ & $x_{1,4}$ & $x_{1,5}$ & $x_{1,6}$ &
          $x_{1,7}$ & $x_{1,8}$ & $x_{1,9}$ & $x_{1,10}$  \\
          
          $x_{2,1}$ & $x_{2,2}$ & $x_{2,3}$ & $x_{2,4}$ & $x_{2,5}$ & $x_{2,6}$ &
          $x_{2,7}$ & $x_{2,8}$ & $x_{2,9}$ & $x_{2,10}$ \\
          
          $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ &
          $\vdots$ & 
          $\vdots$ & $\vdots$ & $\vdots$ \\
          
          $x_{p,1}$ & $x_{p,2}$ & $x_{p,3}$ & $x_{p,4}$ & $x_{p,5}$ & $x_{p,6}$ &
          $x_{p,7}$ & $x_{p,8}$ & $x_{p,9}$ & $x_{p,10}$ 
          
        \end{tabular}
      \end{table}
  \end{itemize}
  
\end{frame}


\begin{frame}
  \frametitle{Approaches}
  \begin{itemize}
    \item Gene-by-gene hypothesis testing
      \begin{itemize}
        \item Treating each gene independently of others.
   %     \item Many properties may be shared among genes.
        \item Goal: find statistically significant associations of biological
          conditions. 
        \item Genes are deemed to be interesting if the $p$-value is small.
        \item Method: $t$-tests, moderated $t$-tests, ROC, $F$-test.
      \end{itemize}
    \item Machine learning
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{$t$-tests}
  
  \[t_g = \frac{\mu_{x}-\mu_y}{\sqrt{\sigma_x^2 - \sigma_y^2}}\]
  
Drawback: 
  \begin{itemize}
    \item Parametric assumptions hard to justify with few arrays.
    \item The variance in small samples might be noisy. 
    \item Genes with small fold-change might be significant from statistical,
      not biological point of view. 
  \end{itemize}
\end{frame}

\section{Moderated $t$-tests}
\begin{frame}
  \frametitle{Moderated $t$-statistics}
  \begin{itemize}
    \item Rather than estimating within-group variability for each gene, pool
      the global information from all other genes.
    \item Advantage: eliminate occurrence of accidentally large $t$-statistics
      due to accidentally small within-group variance.
  \end{itemize}
\end{frame}


% moderated t and limma
\begin{frame}[fragile]
  \frametitle{Moderated $t$-statistics}
  Using empirical Bayesian approach to estimate:
  \begin{itemize}
    \item Overall estimate variation $s_0^2$.
    \item Per-gene deviation variation $s_g^2$. 
    \item Shrinkage variation
      \[
        \tilde{s}_g^2 = \frac{d_0s_0^2+d_g s^2_g}{d_0+ d_g}
      \]
      
    \item Contrast estimator $\hat{\beta}_g$ -- the difference in means
      between two classes.
    \item Moderated $t$-statistics:
      \[
        \tilde{t_{g}} = \frac{\hat{\beta}_{g}}{\tilde{s}_g \sqrt{\nu_{g}}}
      \]
  \end{itemize}
\end{frame}


\begin{frame}[fragile]
  \frametitle{Using limma}
  \begin{enumerate}
    \item Define a design matrix to establish parameters of linear model
      \Rfunction{model.matrix}. 
    %\item Define contrast model if needed \structure{makeContrasts}.
    \item Fit a linear model for each gene based on the given design matrix
      (and a contrast matrix): \Rfunction{lmFit()}.
    %\item Fit the contrast model if needed \structure{contrasts.fit}
    \item Use function \Rfunction{eBayes} to get moderated $t$-statistics
      and relevant statistics.
  \end{enumerate}
  
\end{frame}

\section{Using limma}
% linear model
\begin{frame}[fragile]
  \frametitle{Deriving linear models}
   Suppose we define a design matrix as the following:
   {\small
   \begin{center} 
     \begin{tabular}{c c c}
         \texttt{sample} $i$ & \texttt{(intercept)} & \texttt{mol.biolNEG} \\
         \hline
         \texttt{NEG}    & 1 & 1 \\
         \texttt{BCR/ABL}        & 1 & 0 \\
         \texttt{NEG}    & 1 & 1 \\
         \vdots              & \vdots & \vdots \\
     \end{tabular}
   \end{center}
   }
   Each gene $Y_{j}$ for all sample $ i$, the expression level can be
   expressed by 
   \[
     \left[
     \begin{array}{c}
       Y_{NEG_i, j} \\ Y_{BCR/ABL_i, j}
     \end{array}
     \right] = 
     \left[
     \begin{array}{cc}
       1 & 1 \\
       1 & 0 
     \end{array}
     \right]
     \left[
     \begin{array}{c}
       \beta_{intercept} \\ \beta_{mol.biolNEG}
     \end{array}
     \right] + \epsilon
   \]
   
   \[ \Rightarrow \beta_{mol.biolNEG} = Y_{BCR/ABL_i, j} - Y_{NEG_i, j} + \epsilon \]
   \[ y_j = \beta_{intercept} + \beta_{mol.biolNEG} a_{ij} + \epsilon \]
   \[ \Rightarrow y_j = \beta_0 + \beta_1 a_{ij} + \epsilon \]
   
\end{frame}

\begin{frame}[fragile]
  \frametitle{Using limma}

  Step 1: 
   \begin{block}{Code: define design matrix and contrast model}
{\small
<<design matrix, eval=FALSE>>=
library(limma)
design <- model.matrix( ~mol.biol, ALLfilt_bcrneg)

@ 
}
  \end{block}
  Step 2: 
  \begin{block}{Code: linear models and eBayes}
{\small
<<linear models and ebayes, eval=FALSE>>=
fit1 <- lmFit(exprs(ALLfilt_bcrneg), design)
fit2 <- eBayes(fit1)
topTable(fit2, coef=2, adjust.method="BH", 
         number=5)
@     
}
  \end{block}
\end{frame}

\begin{frame}[fragile]
  \frametitle{Deriving linear models}
   Suppose we define a design matrix as the following:
   {\small
   \begin{center} 
     \begin{tabular}{c c c}
         \texttt{sample} $i$ & \texttt{mol.biolBCR} & \texttt{mol.biolNEG} \\
         \hline
         \texttt{BCR/ABL}        & 1 & 0 \\
         \texttt{BCR/ABL}        & 1 & 0 \\
         \texttt{BCR/ABL}        & 1 & 0 \\
         \vdots              & \vdots & \vdots \\
         \texttt{NEG}    & 0 & 1 \\
         \texttt{NEG}    & 0 & 1 \\
         \texttt{NEG}    & 0 & 1 \\
         \vdots              & \vdots & \vdots \\         
     \end{tabular}
   \end{center}
   }
    
   \[ y_i = \beta_1 a_{ij} + \beta_2 b_{ij} + \varepsilon_i \]
   
\end{frame}

\begin{frame}[fragile]
  \frametitle{Using limma}

  Step 1: 
   \begin{block}{Code: define design matrix and contrast model}
{\small
<<design matrix, eval=FALSE>>=
library(limma)
design <- model.matrix( ~0+mol.biol, ALLfilt_bcrneg)
colnames(design) <- c("BCR_ABL", "NEG")
contr <- makeContrasts(BCR_ABL-NEG, levels=designs)
# contr <- c(1, -1)
@ 
}
  \end{block}
  Step 2: 
  \begin{block}{Code: linear models and eBayes}
{\small
<<linear models and ebayes, eval=FALSE>>=
fit <- lmFit(exprs(ALLfilt_bcrneg), design)
fit1 <- contrasts.fit(fit, contr)
fit2 <- eBayes(fit1)
topTable(fit2,  adjust.method="BH",  number=5)
@     
}
  \end{block}
\end{frame}

\begin{frame}
  \frametitle{$t$-tests vs. moderated $t$-tests}
  \begin{itemize}
    \item In larger sample size, there is not big difference between the
      ordinary and the moderated tests.
    \item For smaller sample size the difference will be larger.
  \end{itemize}
  
  The empirical Bayes moderation is more useful in cases with fewer replicates.
\end{frame}

\begin{frame}
  \frametitle{$t$-tests vs. moderated $t$-tests}
<<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
@ 
<<ttest, echo=FALSE>>=
tt <- rowttests(ALLfilt_bcrneg, "mol.biol")
@ 
<<limma, echo=FALSE>>=
library(limma)
design <- model.matrix( ~mol.biol, ALLfilt_bcrneg)
fit1 <- lmFit(exprs(ALLfilt_bcrneg), design)
fit2 <- eBayes(fit1)
@ 
\begin{figure}
  \centering
<<plotComparison1, echo=FALSE, fig=TRUE, include=TRUE>>=
library(lattice)
print(xyplot(-fit2$t[,2] ~ tt$statistic, 
       main="79 samples",
       xlab="t-statistics from two-sample t-test", 
       ylab="t-statistics  from moderated t-test (limma)",
       panel=function(...) { 
         panel.xyplot(...); panel.grid(); 
         panel.lmline(..., col="red"); panel.abline(a=0,b=1, col="blue") }))
#plot(-log10(tt$p.value), -log10(fit2$p.value[, 2]),
#     xlab= "-log(10)(p) from two-sample t-test",
#     ylab="-log(10)(p) from moderated t-test (limma)",
#     pch=20, col="blue", main="79 samples")
#abline(c(0,1), col="red")
@ 
\end{figure}
\end{frame}



\begin{frame}
  \frametitle{$t$-tests vs. moderated $t$-tests}
<<subs, echo=FALSE>>=
subs = c(35, 65, 75, 1, 69, 71)
ALLsubs <- ALLfilt_bcrneg[, subs]
tt2 <- rowttests(ALLsubs, "mol.biol")
fit3 <- eBayes(lmFit(exprs(ALLsubs), design=design[subs,]))
@ 
  \begin{figure}
    \centering
<<plotComparison2, echo=FALSE, fig=TRUE, include=TRUE>>=
print(xyplot(-fit3$t[,2] ~ tt2$statistic, 
       main="6 samples -- 3 for each group",
       xlab="t-statisitcs from two-sample t-test", 
       ylab="t from moderated t-test (limma)",
       panel=function(...) { 
         panel.xyplot(...); panel.grid(); 
         panel.lmline(..., col="red"); panel.abline(a=0,b=1, col="blue") }))
##plot(-log10(tt2$p.value), -log10(fit3$p.value[, 2]),
##     xlab= "-log(10)(p) from two-sample t-test",
 #    ylab="-log(10)(p) from moderated t-test (limma)",
 #    pch=20, col="blue", main="6 samples -- 3 for each group")
#abline(c(0,1), col="red")
@ 
  \end{figure}
\end{frame}

\section{$p$-value Correction}
\begin{frame}
  \frametitle{$p$-value corrections}
  \begin{itemize}
    \item Basic idea: reduce critical value used to reject.
      
    \item Trade-off between sensitivity and specificity.
      
    \item Approaches implemented in the \Rpackage{multtest} package:
      \begin{itemize}
        \item criteria for error rate control include family-wise error rate
          (FWER) and false discovery rate (FDR).
        \item Permutation-based maxT methods.
      \end{itemize}
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Lab activity}
  \begin{itemize}
  \item Chapter 6 and 7 in {\it Bioconductor Case Studies}.
  \item Goals: get familiar with functions provided by \Bioconductor{}
    packages to perform differential expression analysis. 

  \end{itemize}
\end{frame}

\section{Resources}

\begin{frame}{Resources}
   \begin{itemize}
    \item G.K. Smyth, Linear models and empirical Bayes methods for
        assessing differential expression in microarray experiments,
        {\it Statistical Applications in Genetics and Molecular Biology},
        3(1), 2004. 
     \item G. K. Smyth, {\it limma: Linear Models for Microarray Data},  
       Bioconductor package vignette, 2005.  
        
     \item Florian Hahne et. al., {\it Bioconductor Case Studies}, Springer,
       2007. 

  \end{itemize} 
\end{frame}

\end{document}
