%\VignetteIndexEntry{SOLUTIONS: Introduction to R}
%\VignettePackage{SeattleIntro2010}

\documentclass{article}
\usepackage[colorlinks=true,linkcolor=blue]{hyperref}
\usepackage{theorem}
\usepackage{underscore}
\usepackage{color}
\usepackage{url}
\usepackage{subfigure}

\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Rmethod}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textsf{#1}}}
\newcommand{\Rfunarg}[1]{{\texttt{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}

\newcommand{\code}[1]{{\texttt{#1}}}
\newcommand{\file}[1]{{\texttt{#1}}}

\newcommand{\software}[1]{\textsf{#1}}
\newcommand{\R}{\software{R}}
\newcommand{\bioc}{\software{Bioconductor}}

%% Excercises and Questions
\theoremstyle{break} \newtheorem{Ex}{Exercise}
\theoremstyle{break} \newtheorem{Q}{Question}
%% And solution or answer
\newenvironment{solution}{\color{blue}}{\bigskip}

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

\title{Introduction to R}
\author{Nishant Gopalakrishnan, Martin Morgan}
\date{9 December, 2010}
\maketitle
\tableofcontents
\listoffigures

<<setup, echo=FALSE>>=
options(width=60)
olocale = Sys.setlocale(locale="C")
vers <- R.Version()
vers <- paste(vers$major, vers$minor, sep = ".")
@

\section{Introduction}

Microarray expression data from 128 individuals with acute lymphoblastic
leukemia has been provided as a comma separated file \Robject{exprsMat.csv}. 
Several covariates such as age, sex, type, stage of the disease etc. have been 
provided as another comma separated file \Robject{pData.csv}.
Our goal in this exercise is to
\begin{itemize}
\item Read in the expression data and the covariates into \R{} objects. 
\item Perform data manipulations such as subsetting to arrive at a smaller 
expression dataset with samples that we are interested in.
\item Generate MA plots from our subsetted data using the \Rpackage{lattice} 
package.
\end{itemize}

\section{Load Microarray Expression data into an \R{}session}

We shall proceed to load the data for microarray expression values and the 
covariates from the supplied csv files into an \R{} session. The expression data 
is contained in the file \Robject{exprsMat.csv} and the covariates are in the file
\Robject{pData.csv}. These files are located in the \Robject{extdata} folder of the 
\Rpackage{SeattleIntro2010} package. 

A convenient way to get to the path of the files is to make use of the 
\Rfunction{system.file} function after loading the \Rpackage{SeattleIntro2010}
package using the \Rfunction{library} function. 

\R{} provides several functions such as \Rfunction{read.table} for reading in 
meta data files into a \Robject{data.frame} with appropriate column and row names 
from the header information provided in the file. Convenience functions such as
\Rfunction{count.fields} are also available for discovering problems in files
(such as certain rows in the file having different number of fields) when using 
the \Rfunction{read.table} function.

\begin{Ex}
\begin{itemize}
  \item Use the \Rfunction{system.file} to locate the path to the files 
  \Robject{exprsMat.csv} and \Robject{pData.csv}
  \item Make use of the \Rfunction{count.fields} function on the \Robject{pData.csv}
  file to ensure that the file has the same number of fields in each line of the 
  file.
  \item Use the \Rfunction{read.table} function to read in the experimental 
  meta data into an \R{} variable \Robject{pdOrig}.
  \item Observe the effect, the \Robject{check.names} argument of
  \Rfunction{read.table} set to FALSE has on reading in the column names of the 
  variable \Robject{pdOrig}. The function \Rfunction{colnames} can be used to
  extract the column names of the \Robject{pdOrig}. Hint: Observe the 
  column name \Robject{pmol biol}.
  \item Use the \Rfunction{read.table} function to read in the expression matrix
  into a variable \Robject{expOrig} with the argument
  \Robject{check.names=FALSE}. 
  \end{itemize}
\end{Ex}

\begin{solution}
<<readData>>=
library(SeattleIntro2010)
phenoPath <- system.file( "extdata", "pData.csv", package="SeattleIntro2010") 
pdOrig <- read.table(phenoPath, check.names = FALSE)
colnames(pdOrig)
expPath <- system.file("extdata", "exprsMat.csv", package="SeattleIntro2010") 
expOrig <- read.table(expPath, check.names = FALSE)
@
\end{solution}

\section{Subset data}

The \Robject{expOrig} variable is a \Robject{data.frame} that contains the expression
intensities for the experiment. The column names of the \Robject{data.frame}
represent the sample Id's and the row names represent the genes that are being
studied. A convenient function to look at just the first couple of entries in 
large \Robject{data.frame} is to use the \Rfunction{head} function.

The \Robject{pdOrig} variable is a \Robject{data.frame} with columns representing the
various covariates that describe the experiment (age, sex, etc.). The row names 
correspond to the sample Id's. The column BT is a factor indicating the tumour 
cell type (B1, B2, T1 ,T2 etc. with B indicating B-cell type and T indicating T-cell type). 
Similarly the column \Robject{mol biol} is a factor indicating the molecular 
biology of the cancer. (BCR/ABL, NEG, E2A/PBX1 etc.) The \Robject{mol biol}
column  can be accessed using \Rfunction{pdOrig[["mol biol"]]}.

In this section, we shall try to make use of some of concepts of indexing,
subsetting, factors etc. that we discussed earlier to split our \Robject{expOrig} 
and \Robject{pdOrig} variables to narrow down to the samples that we are 
interested in.  We are specifically interested in B-cell tumours with molecular
biology type "NEG" or "BCR/ABL". We proceed to subset our \Robject{expOrig} and
\Robject{pdOrig} variables to narrow down to our samples of interest.

\begin{Ex}
\begin{itemize}
  \item Use the \Rfunction{grep} function on the \Robject{BT} column in the
 \Robject{pdOrig} \Robject{data.frame} to identify the B cell tumours.
  
  \item Identify the samples that are "NEG" or "BCR/ABL" molecular
  biology type using the column \Robject{mol biol}, the \Rfunction{\%in\%}
  function and the \Rfunction{which} functions in \R{}.

  \item Identify samples that are both B cell tumours and are "BCR/ABL" or "NEG"
   using the \Rfunction{intersect} function on the indices that we have 
   previously computed.
  
  \item Subset the meta data variable \Robject{pdOrig} to create a new variable
  \Robject{psubData}. Note:The rows of the data.frame represent samples.(Subset 
   along rows)

  \item Subset the expression matrix \Robject{expOrig} to create a new
  \Robject{data.frame} \Robject{exprsMat} with only the samples 
  with the cell and the tumour types we interested in. Note:The columns of the 
  expression matrix represents samples.(Subset along columns)
  \end{itemize}
\end{Ex}

\begin{solution}
<<subsetPhenoData>>=
bcell <- grep("^B", as.character(pdOrig$BT))
types <- c("NEG", "BCR/ABL")
moltyp <- which(as.character(pdOrig[["mol biol"]]) %in% types)
indx <- intersect(bcell, moltyp)
psubData <- pdOrig[indx,]
exprsMat <- expOrig[, indx]
@
\end{solution}

\section{Recode the factor levels}

The covariate data for some variables, for example \Robject{BT} is represented
using a variable of type factor with distinct levels. These variables can only 
take a distinct number of categorical values. The levels can be obtained by
using the \Rfunction{levels} function on the factor variable. Operations such as
subsetting on a factor variable subsets the factor but leaves the levels of the
variable unchanged. In this exercise, we shall take a look at the levels of the
factor variables that we have just subsetted and attempt to correct this problem.

\begin{Ex}
\begin{itemize}
  \item Observe the levels for the \Robject{mol biol} and \Robject{moltyp} variables
   Do you notice any problem ?.
  \item Recode the factor levels for the \Robject{mol biol} and moltyp variables
  using the \Rfunction{factor} function.
   \end{itemize}
\end{Ex}

\begin{solution}
<<recodeFactors>>=
psubData$BT <-  factor(psubData$BT)
levels(psubData$BT)
psubData[["mol biol"]] <- factor(psubData[["mol biol"]])
levels(psubData[["mol biol"]])
@
\end{solution}

\section{Compute summary statistics}

\R{} includes several functions that allows you to do a lot while writing only few
lines of code. A good example is the \Rfunction{aggregate} function that splits
the data into subsets, computes summary statistics for each section and
returns the result in a convenient form. For more details regarding this
function, please type in \Rfunction{help("aggregate")} into an \R{} session. We
will be making use of the \Rfunction{formula} and the \Rfunction{data.frame}
methods for the \Rfunction{aggregate} function in this example. Another useful
function for creating contingency tables that we will be using is the 
\Rfunction{xtabs} function.

We proceed to create some summary statistics on the meta data variable
\Robject{psubData} using the \Rfunction{aggregate} and \Rfunction{xtabs}
functions.

\begin{Ex}
\begin{itemize}
  \item Find the average age of male and females in our subsetted metadata
  variable \Robject{psubData} for the NEG and BCR/ABL groups using the  
  \Robject{data.frame} interface for the \Rfunction{aggregate} function. The
  table generated by using the \Rfunction{aggregate} should look similar to the 
  one found below. Hint: Try passing \Robject{na.rm = TRUE} to the aggregate
  function
  
  \begin{table}[ht]
\begin{center}
\begin{tabular}{rllr}
  \hline
 & sex & molBiol & age \\ 
  \hline
1 & F & BCR/ABL & 39.94 \\ 
  2 & M & BCR/ABL & 40.50 \\ 
  3 & F & NEG & 29.75 \\ 
  4 & M & NEG & 24.86 \\ 
   \hline
\end{tabular}
\end{center}
\end{table}

  \item Recalculate the average age of male and females in our subsetted metadata
  variable \Robject{psubData} for the NEG and BCR/ABL groups, this time using the  
  \Robject{formula} interface for the \Rfunction{aggregate} function. Make sure
  that the results are identical to the one from the previous step.

  \item  The column \Robject{relapse} in \Robject{psubData} is a logical vector 
  indicating whether the patient had a relapse or not. Create a contingency table 
  of the number of subjects that have had a relapse for the samples included in
  \Robject{psubData} using the \Rfunction{xtabs} function and the covariates
  \Robject{relapse}, \Robject{mol biol} and \Robject{sex}. The table generated
  by using the \Rfunction{xtabs} function should look similar to the 
  one found below.


\begin{table}[ht]
\begin{center}
\begin{tabular}{rrr}
  \hline
 & BCR/ABL & NEG \\ 
  \hline
F & 7.00 & 3.00 \\ 
  M & 9.00 & 18.00 \\ 
   \hline
\end{tabular}
\end{center}
\end{table}

   \end{itemize}
\end{Ex}

\begin{solution}
<<contingencytables>>=

aggregate( psubData[, "age", drop = FALSE], by= list("sex"= psubData$sex,
            "molBiol"= psubData[["mol biol"]]),  FUN = mean,
        na.rm = TRUE )
aggregate(age ~ sex + `mol biol`, data = psubData, FUN = mean)

xtabs( relapse ~ sex + `mol biol`, data = psubData)
@
\end{solution}

\section{Data visualization}

The \Rpackage{lattice} package implements the Trellis graphics system and
provides high-level functions for visualization of multivariate data.
A typical call to the lattice function \Rfunction{xyplot} is shown below.
<<xyplotFormula, eval=FALSE>>=
xyplot(y ~ x | c, data, groups = g)
@

The arguments to a lattice function can be summarized in terms of
\begin{enumerate}
	\item{lattice function: A lattice plotting function such as \Rfunction{xyplot},
	\Rfunction{dotplot} etc.}
    \item{formula: The first argument to a lattice method is a formula. The
	formula for our example is \code{y} \verb+~+ \code{x | c}. If the lattice 
	method takes only a single vector as input, the formula can be expressed as 
	 \verb+~+ \code{ x | c}.}
	 \begin{itemize}
	\item{primary variables: Variables y (Y axis of the plot) and x (X axis of the
	plot) that defines the lattice display separated by the \verb+~+ character.
	}
	\item{conditioning variable: Variable c in the example separated from the
	primary variables by the character \code{|}. The conditioning variable divides the plot
	into separate panels.}
	\end{itemize}
	\item{grouping variable: The variable g in the example. The grouping variable
	segregates data into subgroups within each panel.}
	\item{data: A \Rclass{data.frame} with column names corresponding to the
	variables y, x, c and g.}
\end{enumerate}

We proceed to create M-A plots for the first 20 samples in our subsetted 
expression intensity data.frame \Robject{exprsMat}, using the formula provided below
and the \Rfunction{xyplot} function from the \Rpackage{lattice} package.

\begin{equation}
 A = rowMeans(m)\\
\end{equation}
\begin{equation}
 M_i = m_i - A\\
\end{equation}

where 
\begin{itemize} 
    \item m is the expression \Robject{data.frame}.
    \item i is from 1 to number of samples in the expression
    \Robject{data.frame}
    \item M is variable for the Y axis of the M-A plot
    \item Mi is the column for sample i
    \item A is the variable for the X axis of the M-A plot
\end{itemize}


\begin{Ex}

\begin{itemize}
  \item  Create a new variable \Robject{mat} by subsetting the \Robject{exprsMat}
  variable to include only the first 20 samples. Using the \Rfunction{rowMeans}
  function, create the variable \Robject{A} 

  \item Create a variable \Robject{M} using the equation provided above

  \item Create a \Robject{data.frame} \Robject{df} with variables M, A and the
  sample names in the expression matrix. The column names of the data.frame
  should be "M", "A" and "Sample"  respectively.

  \item Use the \Rfunction{xyplot} function from the lattice package to plot the variable
  M on the Y axis and the variable A on the x axis. Make use of the \Robject{Sample}
  variable as a conditioning variable to split the plot into sub plots one each
  for each sample.
  \end{itemize}
\end{Ex}

\begin{solution}
<<lattice plots>>=
library(lattice)
mat <- exprsMat[,1:20]
A <- rowMeans(log2(mat))
M <- log2(unlist(mat)) - A
Sample <- rep(colnames(mat), each=nrow(mat))
df <- data.frame(M, A, Sample, row.names=NULL, check.names = FALSE)
plt <-  xyplot(M ~ A | Sample, df, panel=panel.smoothScatter)

@
\end{solution}

<<xyPlotPrint,fig=TRUE,include=FALSE,prefix=FALSE,width=4,height=4,echo=FALSE>>=
print(plt)
@
\begin{figure}[p]
          \includegraphics[width=.8\textwidth]{xyPlotPrint}
          \caption{M-A plot for 20 Samples produced using \Rfunction{xyplot}}

\end{figure}


\section{Session information}

<<sessionInfo, results=tex, echo=FALSE>>=
toLatex(sessionInfo())
@

\end{document}

