%\VignetteIndexEntry{Lymphoma: Exploring cDNA data}
%\VignetteDepends{Biobase,multtest,vsn,marray,limma,arrayMagic}
%\VignetteKeywords{Course in Practical Microarray Analysis}
%\VignettePackage{}

\documentclass[a4paper, twoside, 11pt]{article}

% packages %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{latexsym}
\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage[dvips]{graphicx}
\usepackage[round]{natbib}
\usepackage{alltt}


% section style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{sectsty}
\usepackage{caption}
\renewcommand{\captionfont}{\it\sffamily}
\renewcommand{\captionlabelfont}{\bf\sffamily}
\allsectionsfont{\sffamily}

% color %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{color}
\definecolor{gray}{rgb}{0.9,0.9,0.9}

% page style %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[a4paper,left=2cm,right=2cm,top=2cm,bottom=2cm,nohead]{geometry}
\usepackage{bioCexercises}

\setlength{\parskip}{1.5ex}
\setlength{\parindent}{0cm}
\pagestyle{empty}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% begin document %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}
\sffamily


% header %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{--- Exploring cDNA Data---\\}
\author{Achim Tresch, Andreas Buness, Tim Bei\ss barth, Florian Hahne,
  Wolfgang Huber}
\maketitle
\thispagestyle{empty}

\textbf{Please note: the package lymphoma and this vignette are deprecated. They will not
be supported after Bioconductor release 2.0.}

The following exercise will guide you through the first steps of a
spotted cDNA microarray analysis.  These steps comprise loading data
into R/Bioconductor, quality control of the measurements, and
preprocessing of the raw data via normalization up to the generation of
very simple gene lists.  Make extensive use of the
\code{help(<object>)} command to find information about particular
objects.  Use the \code{vignette(<package>)} command to get
introductory material for a certain package.\\

<<no.nonsense,eval=FALSE,echo=FALSE,results=hide>>=
q(save="no")
@

\begin{exercises}
%%--------------------Ex.1------------------------------ 
\item {\bf Preliminaries.}  To go through this exercise, you need to
  have installed R >= 2.1.0, the release 1.6 versions of the
  Bioconductor libraries Biobase, marray, multtest, limma, vsn, and
  arrayMagic and the library lymphoma, which contains the excercises
  and part of the lymphoma data set.
%(see http://llmpp.nih.gov/lymphoma/).

<<libraries,echo=TRUE,results=hide>>= 
  library("vsn")
  library("multtest")
  library("marray")
  library("arrayMagic")
  library("RColorBrewer")
@

<<loadData, echo=FALSE>>=
data(lymphoma, package="lymphoma")
@ 
%%--------------------read.table one file------------
\item {\bf Reading data files.}  
Considerable attention should be payed to data import. Experience
tells that this is one of the most error prone steps.

\begin{exercise}

  \item Your data has to be stored in one folder, with one file
  corresponding to one sample. Our example folder is located at
  \code{<R library path>/lymphoma/extdata/}. You can find the path to
  that location using \code{system.file}.
  
<<setlibpath,echo=TRUE>>= 
 path = system.file("extdata", package="lymphoma")
 path
@
  The file names are \code{lc7b047rex.DAT, lc7b048rex.DAT, ...}.

<<dirlibpath,echo=TRUE>>= 
 dir(path, pattern=".DAT$")
@
  \item Open the file \file{lc7b048rex.DAT} in a text editor.  This is
  the typical file format for the results from the image analysis on a
  cDNA slide. Different image analysis programs use slightly different
  conventions and column headings, but we will describe an import
  method which is suited to the most common software (Genepix, Spot,
  ...??).\\


  \item For a both easy and flexible import of the data, there has to
  be a description file. The description file contains a table with
  all the hybridization data file names in one column and possibly
  additional sample information in further columns.  Create a
  tab-delimited text file of this kind. We have done this for you, the
  description file is named \code{phenoData.txt}. You may examine its
  structure with any text editor. There is a convenient method for
  converting such a table into a \code{phenoData} object.

<<phenodata.load>>=  
lymphenoData = read.phenoData(file.path(path,"phenoData.txt"))
@

\item The \code{phenoData} object is used to simulatneously import all
data files.

<<lymph.load,eval=T,results=hide, eval=FALSE>>=
lymphRaw = readIntensities(pData(lymphenoData), loadPath=path,
           fileNameColumn="fileName", slideNameColumn="slideNumber",
           type="ScanAlyze")
@

<<lymph.normalize,eval=T,results=hide, eval=FALSE>>=
lymphNormvsn = normalise(lymphRaw, subtractBackground = T, method="vsn",
spotIdentifier = "SPOT")
@

<<lymphoma.set,eval=FALSE>>=
lymphoma=as.exprSet(lymphNormvsn)
@

\end{exercise}



%%-------------------description of class exprSet -----------------
\item {\bf The Bioconductor class \code{exprSet}. }
\begin{exercise}
\item The object \code{lymphoma} is of class \code{exprSet}.  This
class is the standard representation of a microarray experiment in
Bioconductor.  It consists of the objects ("slots")\\

\begin{tabular}{lcl}
      exprs &:& A spots $\times$ samples matrix containing the
      expression levels\\ se.exprs &:& A spots $\times$ samples matrix
      containing an estimate of the standard error\\ & & of each
      single spot measurement\\ phenoData &:& An object of class
      \code{phenoData}, essentially a data frame containing\\ & &
      phenotypical information about the samples that were
      hybridized\\ annotation &:& Textual annotation\\ description &:&
      Object of class \code{MIAME} which incorporates those
      MIAME-entries\\ & & that are not covered by other objects of the
      \code{exprSet} class\\ notes &:& Text containing additional
      remarks\\
\end{tabular}


Slots can be accessed directly with "@"
(e.g. \code{lymphoma@phenoData}), but one should use the accessor
methods for the class \code{exprSets}. See \code{help(exprSets)} for
details. The most interesting objects to us are the expression matrix,
given by \code{exprs(lymphoma)} and the data frame with the phenotype
data, \code{pData(lymphoma)}.  Have a look at them.\\

<<pData1,echo=FALSE,eval=TRUE>>=
spotnames = paste("gene",1:nrow(exprs(lymphoma)),sep="")
rownames(exprs(lymphoma)) <- spotnames
@

\code{> dim(exprs(lymphoma))  \# genes $\times$ samples}
<<pDatashow,echo=FALSE>>=
dim(exprs(lymphoma))
@
<<pDatashow>>=
exprs(lymphoma)[1:3,1:6] # show upper left corner of the expression matrix
@
\code{> dim(pData(lymphoma))  \# samples $\times$ descriptors}
<<pDatashow,echo=FALSE>>=
dim(pData(lymphoma))
@
<<pDatashow,echo=TRUE,eval=FALSE>>=
pData(lymphoma)  # original phenoData for the \code{lymphoma} experiment
@
<<pDatashow,echo=FALSE,eval=TRUE>>=
pData(lymphoma)[1:3,]
@
$\vdots$\\

\item You might want to add a column containing the hybridization colour.
<<pDatashow,echo=TRUE,eval=TRUE,results=hide>>=
colour = c(rep("red",8),rep("green",8))
pData(lymphoma) = cbind(pData(lymphoma),colour)
@

<<pDatashow,echo=FALSE,eval=TRUE>>=
pData(lymphoma)[1:3,]
@
$\vdots$\\
\end{exercise}



%%-------------------histogram and scatterplot -----------------
\item {\bf Simple plots.}
\begin{exercise}
\item 
We will perform some elementary diagnostic plots for quality control.
Most analyses are carried out on the log transformed data, so \code{lymphoma} 
contains (generalized) log transformed expression values. For convenience, we 
extract these values into another variable.

<<logtransformation>>=
logexpr = exprs(lymphoma) #log(exprs(lymphoma))
#logexpr[is.infinite(logexpr)]= NaN  # remove infinite values
@

It is possible to examine each single channel or to focus on log
ratios or difference between them. Notice that in our expSet the data
for the red channel are in columns 1 to 8 and data for the green
channel in columns 9 to 16. For the following plots we will first use the
two channels of slide number 5. Later we will also look at the log ratios.

Now produce a histogram and a density plot
of the log intensities in channel 1 of slide 5.\\

<<logplots,echo=T>>=
slidenr=5
ch1=logexpr[,slidenr]; ch2=logexpr[,slidenr+8] 
@
<<logplots2,eval=F,results=hide>>=
plot(hist(ch1),main="Histogram", col="blue")
plot(density(ch1),main="Densityplot")
@

<<grafik1,echo=FALSE,results=hide>>=
pdf("histogramdensity.pdf",height=5,width=8.5)
par(mfrow=c(1,2))
hist(ch1,col="blue",main="Histogram")
plot(density(ch1),main="Densityplot")
dev.off()
@

\begin{center}
\includegraphics[width=0.7\textwidth]{histogramdensity.pdf}
\end{center}

Here are some plots that help to detect bad hybridizations. Compare the 
log intensity boxplots of the slides

<<compare,eval=T,results=hide>>=
boxplot(split(t(logexpr),1:ncol(logexpr)),col=as.vector(pData(lymphoma)[,"colour"]),
 main="boxplot")
@

<<grafik2,echo=F,results=hide>>=
pdf("boxplot.pdf",height=5,width=8.5)
par(mfrow=c(1,2))
boxplot(split(t(logexpr),1:ncol(logexpr)),col=as.vector(pData(lymphoma)[,"colour"]),
 main="Boxplot")
qqplot(ch1,ch2,main="q-qplot")
dev.off()
@

\begin{center}
\includegraphics{boxplot.pdf}
\end{center}

A convenient way to compare the expression distributions between two
samples is a quantile-quantile plot. For equal distributions it should
be a diagonal line.

<<qqplot,results=hide>>=
qqplot(ch1,ch2,main="q-q plot")
@

Do you think the q-qplot of slide number 5 looks OK?\\

\item Save one of the plots as a PDF. Copy and paste it into an MS-Office application.

<<saveplot,eval=F>>=
pdf(file="savedplot")
qqplot(ch1,ch2,main="q-q plot")
dev.off()
@

\item A more detailed view is provided by a scatter plot and its corresponding M-A plot

\code{
> plot(ch1,ch2,pch=".",main="Scatterplot")\\
> abline(a=0,b=1,col="blue"); abline(a=1,b=1,col="red"); abline(a=-1,b=1,col="red")\\
> plot((ch1+ch2)/2,ch2-ch1,pch=".",xlab="A",ylab="M",main="M-A plot")\\
> abline(h=0,col="blue"); abline(h=1,col="red"); abline(h=-1,col="red")
}

<<scatterplot2,echo=F,results=hide>>=
pdf(file="grafik3.pdf",height=5,width=8.5)
par(mfrow=c(1,2))
plot(ch1,ch2,pch=".",main="Scatterplot"); abline(a=0,b=1,col="blue")
abline(a=1,b=1,col="red"); abline(a=-1,b=1,col="red")
plot((ch1+ch2)/2,ch2-ch1,pch=".",xlab="A",ylab="M",main="M-A plot"); 
abline(h=0,col="blue")
abline(h=1,col="red"); abline(h=-1,col="red")
dev.off()
@
\vspace{2mm}

\begin{center}
\includegraphics[width=0.8\textwidth]{grafik3.pdf}
\end{center}

\item We hid one "outlier slide" among the original lymphoma data. 
Find it using diagnostic plots!

\end{exercise}



\item {\bf Normalization}
\begin{exercise}

\item Before we started analysis, we tacitly normalized our data (cf.
\code{normalise(lymphraw...)}). Try two other commonly used
normalization methods:

<<lymph.normalize,eval=T,results=hide>>=
lymphNormloess = normalise(lymphRaw, subtractBackground = T, method="loess", spotIdentifier = "SPOT")
logexpr2=as.exprSet(lymphNormloess)
lymphNormquantile = normalise(lymphRaw, subtractBackground = T, method="quantile", spotIdentifier = "SPOT")
logexpr3=as.exprSet(lymphNormquantile)
@

\item These commands take their time! You can save the results into a file with
the \func{save} function, and later restore them with the \func{load}
function. In MS-Windows, you can use the GUI for the latter.

\item Compare the results of the variance stabilization method to the loess method! 
Which Plots are appropriate for that?

\end{exercise}



\item {\bf Testing for differential expression}
\begin{exercise}

\item  Let us now get the log-ratios of the 2 channels
for all 8 slides in our data set. This can be done quite easily using
function \code{getExprSetLogRatio}. Look into its manual and also find
out about \code{exprSetRG} objects, which extend regular
\code{exprSets} to include information about the two different
channels.

<<log-ratios>>=
lymphRatio = getExprSetLogRatio(lymphNormvsn)
logRatios = exprs(lymphRatio)
@

Finally we are ready to calculate test statistics and to select
genes.  {\em Note:} The number of replicates (4 versus 4) that we are
considering here is too small to derive significant conclusions about
individual genes.  The full data set contains many more chips.  Here
we restrict ourselves to a few of them in order to keep things simple
for the purpose of this course.

\begin{exercise}
\item Look at the built-in function \func{t.test}, and at
  \func{mt.teststat} from the package \code{multtest}. Here, we use
  \func{mt.teststat} to calculate the $t$-test statistic for the
  comparison. The package \code{multtest} provides extensive
  functionality to calculate multiple-testing adjustments.

<<t,width=5,height=5,prefix=FALSE,fig=TRUE,include=FALSE>>=
classlabel = c(0,0,0,0,1,1,1,1)
tStat = mt.teststat(logRatios, classlabel)
summary(tStat)
hist(tStat, breaks=100, col="#fb6090")
@
%
\begin{figure}[htb]
\begin{center}  
  \includegraphics[width=0.4\textwidth]{t.pdf}
  \includegraphics[width=0.4\textwidth]{fdc.pdf}
  \caption{\label{tfdc} a) \textit{left:} histogram of $t$-statistic, 
  b) \textit{right:} false disocvery count.}
 \end{center}
\end{figure}

\item Similar to the FDR (false discovery rate) we estimate
      the significance of the most extreme
      t-values. The false discovery count is calculated for
      several lists of various length of the highest scored genes
      (Fig.~\ref{tfdc} right).

<<fdcCalc>>=
fc = fdc(logRatios, factor(classlabel), teststatfun="rowttests")
@ 

<<fdc,prefix=FALSE,fig=TRUE,include=FALSE>>=
plot(fc$nrgenesel,fc$fdc,main="False Discovery Count", 
     xlab="number of selected genes", ylab="false discovery count", 
     type="b")
@

\item Now we load the spot (gene) description table with function
 \code{read.delim}. You can find it in the same folder as the raw data
 in the file annotationData.txt. 

What does \code{read.delim} do?

<<spotlabels.do>>=
spotDescr = read.delim(system.file("extdata", "annotationData.txt", package="lymphoma"))
@

Let's print the 5 genes with the lowest values of the $t$-statistic
<<low>>=
selection = order(tStat)[1:5]
selection
@
as well as the 5 genes with the highest values of the $t$-statistic
<<high>>=
selection = order(tStat, decreasing=TRUE)[1:5]
selection
@

In a following step we extract the annotation information
for the 5 selected genes and generate a html-report
(Fig.~\ref{candidates}).

<<writeTable>>=
spotIDs = getSpotAttr(lymphRaw)$SPOT[selection] 
geneAnno = spotDescr[spotDescr[,"spotID"] %in% spotIDs,]
write.htmltable(geneAnno,filename="candidateGenes",title="Candidate Genes")
@

\begin{figure}[h]
\begin{center}  
  \includegraphics[width=0.5\textwidth]{candidateGenes}
  \caption{\label{candidates}}
 \end{center}
\end{figure}

\end{exercise}


\item Due to their pervasive power, heatmaps enjoy high popularity 
(although they hardly prove anything). We can produce one that shows
 the expression levels of the 20 genes with highest values of the
 t-statistc in a few lines of R code.

<<draw.heatmap>>=
selection = order(tStat, decreasing=TRUE)[1:20]
heatmap(logRatios[selection,],col=brewer.pal(10,"RdBu"))
@

<<draw.heatmap,echo=F,results=hide>>=
pdf("heatmap.pdf")
heatmap(logRatios[selection,],col=brewer.pal(10,"RdBu"))
dev.off()
@


\begin{center}
\includegraphics[width=0.9\textwidth]{heatmap.pdf}
\end{center}

\end{exercise}



\item {\bf Further quality assessment}
\begin{exercise}

The package \code{arrayMagic} provides additional measures for quality
control. It produces a bunch of graphics which are saved in your
current working directory. Take the time to examine some of them. Have
a look at the vignette \code{arrrayMagicVignette} for details before
proceeding.

<<vignetteArrayMagic,eval=FALSE>>=
vignette("arrayMagicVignette")
@ 

<<qualityArrayMagic,eval=TRUE>>=
qP <- qualityParameters(lymphRaw, lymphNormvsn, resultFileName="qP.txt", spotIdentifier="SPOT", slideNameColumn="fileName")
qualityDiagnostics(lymphRaw, lymphNormvsn, qP)
visualiseHybridisations(lymphRaw[,1], mappingColumns=list(Block="GRID", Column="COL", Row="ROW"))
@

here are some samples of the output:

\begin{center}
\includegraphics[width=0.8\textwidth]{slideDistancesHeatmap.pdf}
\end{center}
\begin{center}
\includegraphics{hybridisation_raw_1.pdf}
\end{center}
\end{exercise}
\end{exercises}
\end{document}
