%
% NOTE -- ONLY EDIT THE .Rnw FILE!!!  The .tex file is
% likely to be overwritten.
%
% \VignetteIndexEntry{Lab 5}
%\VignetteDepends{Biobase. golubEsets, ellipse,lattice, mva, MASS,cluster }
%\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 5: Cluster Analysis Using R and Bioconductor}

\begin{document}

\maketitle

\section*{Introduction}

In this lab we introduce you to various notions of distance and to
some of the clustering algorithms that are available in R. The lab
comes in two parts, in the first we consider different distance
measures while in the second part we consider the clustering
methods. It is worth mentioning that all clustering and machine
learning methods rely on some measure of distance. Almost all
implementations have some measure that they use as a default (this is
usually Euclidean). As a user of this software you need to decide what
distance measures are of interest or are relevant for your situation
and to subsequently ensure that the clustering methods use the
distance metric that \textbf{you} believe is appropriate.

The examples in this section will rely on the data reported in
\citet{Golub99}. The basic comparisons here are between B-cell derived
ALL, T-cell derived ALL and AML. In some cases we will just compare
ALL to AML while in others comparisons between the three groups will
be considered.

<<loadlibs>>=
 library(Biobase)
 library(annotate)
 library(golubEsets)
 library(genefilter)
 library(ellipse)
 library(lattice)
 library(cluster)
 library(mva)
@

We transform the data in much the same way that \citet{Golub99}
did. The sequence of commands needed are given below.

<<setup>>=
data(golubTrain)
X<-exprs(golubTrain)
X[X<100]<-100
X[X>16000]<-16000

mmfilt <- function(r=5, d=500, na.rm=TRUE) {
function(x) {
  minval <-  min(x, na.rm=na.rm)
  maxval <-  max(x, na.rm=na.rm)
  (maxval/minval > r) && (maxval-minval > d)
  }
}

mmfun <- mmfilt()

ffun <- filterfun(mmfun)
sub <- genefilter(X, ffun )
sum(sub) ## Should get 3051

X <- X[sub,]
X <- log2(X)
golubTrainSub<-golubTrain[sub,]
golubTrainSub@exprs <- X

Y <- golubTrainSub$ALL.AML

Y<-paste(golubTrain$ALL.AML,golubTrain$T.B.cell)
Y<-sub("NA","",Y)

@

Now that the data are set up, we are ready to compute distances between
tumor samples and apply clustering procedures to these samples.

\section{Distances}

We first compute the correlation distance between the 38 tumor samples
from the training set, using all 3,051 genes retained after the above
filtering procedure.

%%FIXME: other distances?

<<cordist>>=
r<-cor(X)
dimnames(r)<-list(as.vector(Y),as.vector(Y))
d<-1-r
@


We rely on the \texttt{plotcorr} function from the \textit{ellipse}
package and the \texttt{levelplot} function from \textit{lattice} for
displaying this 38 $\times 38$ distance matrix.

<<plotcorr1,fig=TRUE>>=
plotcorr(r,
main="Leukemia data: Correlation matrix for 38 mRNA samples\n All 3,051 genes")
@

<<plotcorr2,fig=TRUE>>=
plotcorr(r,numbers=TRUE,
main="Leukemia data: Correlation matrix for 38 mRNA samples\n All 3,051 genes")
@

<<levelplot>>=
levelplot(r,col.region=heat.colors(50),
main="Leukemia data: Correlation matrix for 38 mRNA samples\n All 3,051 genes")
@

We can see in this plots that there is not a strong grouping of the
ALL and AML samples. We now consider selecting genes that are good
differentiators between the two groups. We use a $t$-test filter and
require that the $p$-value be below 0.01 for selection (no adjustment
has been made).

%%FIXME: this takes a while and we might want to store the end result
<<getDiffGenes>>=

gtt <- ttest(golubTrainSub$ALL, p=0.01)
gf1 <- filterfun(gtt)
whT <- genefilter(golubTrainSub, gf1)

##subset the genes
gTrT <- golubTrainSub[whT,]

@

Now we can repeat the above computations and plots using the data in
\Robject{gTrT}. There were \Sexpr{sum(whT)} genes selected using the
$t$-test criterion above.

%%FIXME: I'm worried about this plot -- shouldn't their be a negative
%%correlation between the ALL and AML patients.
<<corrDsub>>=
rS <- cor(exprs(gTrT))
dimnames(rS)<-list(gTrT$ALL, gTrT$ALL)
dS<-1-rS

plotcorr(r,
main="Leukemia data: Correlation matrix for 38 mRNA samples\n 609 genes")
@

\section{Multidimensional Scaling}

We now turn to {\em multidimensional scaling} (MDS) for representing
the tumor sample distance matrix. MDS is a data reduction method that
is appropriate for general distance matrices. It is much like
principal component analysis (PCA), but it is not the same -- except
for Euclidean distances.  Given an $n \times n$ distance matrix, MDS
is concerned with identifying $n$ points in Euclidean space with a
{\em similar} distance structure.  The purpose is to provide a lower
dimensional representation of the distances which can better convey
the information on the relationships between objects, such as the
existence of clusters or one-dimensional structure in the data (e.g.,
seriation). In most cases we would like to find a two or three
dimensional reduction that could easily be visualized.

<<mds2,fig=TRUE>>=
mds<- cmdscale(d, k=2, eig=TRUE)
plot(mds$points, type="n", xlab="", ylab="",
main="MDS for ALL AML data, correlation matrix, G=3,051 genes, k=2")
text(mds$points[,1],mds$points[,2],Y, col=as.integer(factor(Y))+1,
     cex=0.8)
@

<<mds3,fig=TRUE>>=
mds3 <- cmdscale(d, k=3, eig=TRUE)
pairs(mds3$points,
main="MDS for ALL AML data, correlation matrix, G=3,051 genes, k=3",
pch=c("B","T","M")[as.integer(factor(Y))], col = as.integer(factor(Y))+1)
@

%%FIXME: put in a call to Rgl --
<<mdsGL>>=

if( require(rgl) && interactive() ) {
  rgl.spheres(x=mds3$points[,1], mds3$points[,2],
        mds3$points[,3], col= as.integer(factor(Y))+1, radius=.01)
 }
@

Note, if we are measuring distances between samples on the basis of
$G=3,051$ probes then we are essentially looking at points in 3,051
dimensional space.  As with PCA, the quality of the representation
will depend on the magnitude of the first $k$ eigenvalues.  To assess
how many components there are, a \textit{scree} plot similar to that
used for principal components can be created.  This plot of the
eigenvalues suggests that much of the information is contained in the
first component. One might consider using either three or four
components as well.

<<scree,fig=TRUE>>=
 mdsScree <- cmdscale(d, k=8, eig=TRUE)
 plot(mdsScree$eig, pch=18, col="blue")
@


\section*{EDD and ROC}

In this section we consider a couple of other tools for visualizing
genomic data. The first is \Rpackage{edd} which stands for expression
density diagnostics. The idea behind edd is to look for genes whose
patterns of expression are quite different in two groups of
interest.

To assess whether two genes have different patterns of expression edd
transforms all expression values to the range 0,1 and then performs
density estimation on them. Thus it is looking for differences in
shape and ignores differences in scale or location. In some situations
location and scale are more important (however there are many
different tools that can be used in those situations).

In this case we will use \Rpackage{edd} to compare those patients with
AML to those with B-cell derived ALL.
We will rely on the merged data set so that we
have a number of samples from each group (38 ALLs and 25 AMLs).

Once we have subset the data to the appropriate cases we next perform
some non-specific filtering to remove any genes that show little
variation across samples. This is an important step when using
\Rpackage{edd}. Since \Rpackage{edd} standardizes the genes (across
samples) to give approximately common centering and scaling
non-informative genes will not be distinguishable from informative
ones. Therefore, we must explicitly remove the non-informative ones
before applying \Rpackage{edd}.

<<edd>>=

 gMs <- golubMerge[,(golubMerge$T.B != "T-cell"
            | is.na(golubMerge$T.B))]

 fF1 <- gapFilter(100, 200, .1)
 ff <- filterfun(fF1)
 whgMs <- genefilter(gMs, ff)
 sum(whgMs)  ##3082
 gMs <- gMs[whgMs,]

 library(edd)

@

We next split the data into two expression sets, one with the AML and
the other with the B-cell ALL samples. For each set we use
\Rfunction{edd} to classify samples into one of a set of
predetermined shapes.

We will use \Rfunction{edd} to perform the classification
and will use the default values of most of the parameters. That means
that we will be using the \Robject{multiCand} for candidate comparison
and \Rfunction{knn} for classification.

When \Robject{multicand} is selected edd selects 100 samples from each
of the eight candidate densities. Each of sample is scaled to have
median 0 and mad 1. These constitute the reference samples and for
each gene we will perform the same transformation and then compare the
resulting data to the reference samples. Any machine learning
technique can in principle be used to make this comparison. In this
example we will use the default classifier which is \Rfunction{knn}.

The density, across samples, for each gene is
found and scaled to have median 0 and mad 1 (as for the reference
samples).

<<eddGM>>=

  gMsALL <- gMs[,gMs$ALL=="ALL"]
  gMsAML <- gMs[,gMs$ALL=="AML"]

  ##set the seed for the rng for reproducibility
  set.seed(1234)

  edd.ALL <- edd(gMsALL, method="knn", k=4, l=2)
  ##doubt
  sum(is.na(edd.ALL))
  edd.AML <- edd(gMsAML, method="knn", k=4, l=2)
  sum(is.na(edd.AML))
  table(edd.ALL, edd.AML)
@

For \Sexpr{sum(is.na(edd.ALL))} of the ALL patients no classification
was made. So properly they should be classified as \textit{doubt}. For
the AML group \Sexpr{sum(is.na(edd.AML))} were classified as
\textit{doubt}.
We can see from the resulting table that there are a number of genes
for which the ALL classification is a mixture (\texttt{mix1}) while
the AML classification is as a Normal random variable.
To further explore these will need to make use of some plots.

<<geneComp>>=

 diff1 <- edd.AML == "N(0,1)" & edd.ALL == ".75N(0,1)+.25N(4,1)"
 table(diff1)  ##should be 47 TRUE
 diff1[is.na(diff1)] <- FALSE
 AMLexprs <- exprs(gMsAML)[diff1,]
 ALLexprs <- exprs(gMsALL)[diff1,]
@

Now we have the two sets of expression data and we can see what the
differences are.
%%FIXME: need to figure out how to get the output of fq.matrows onto
%%the histograms in some sensible fashion

<<transF>>=
 tAML  <- fq.matrows(t(apply(AMLexprs, 1, centerScale)))
 tALL  <- fq.matrows(t(apply(ALLexprs, 1, centerScale)))

 par(mfrow=c(2,2))
 hist(AMLexprs[1,])
 hist(ALLexprs[1,])
 hist(AMLexprs[2,])
 hist(ALLexprs[2,])

 par(mfrow=c(1,1))
@

\subsection*{ROC}

Receiver operator curves are a tool that can be used to assess the
classificaton capabilities of a covariate. In situations where there are
two classes (for example, AML versus ALL) then one can consider
splitting the data into two groups depending on the value of the
covariate. The relative proportions of AML and ALL in the two
resulting groups can be used to assess the capability of that
covariate for splitting the data.


\section*{Hierarchical Clustering}

We can use the \texttt{hclust} function in the \textit{mva} package for
{\em agglomerative hierarchical clustering}. We consider the following
agglomeration methods, i.e., methods for computing {\em between
cluster distances}:

\begin{description}
\item[Single linkage] The distance between two clusters is the minimum
  distance between any two objects, one from each cluster.
\item[Average linkage] The distance between two clusters is the
  average of all pairwise distances between the members of both
  clusters.
\item[Complete linkage] The distance between two clusters is the
maximum distance between two objects, one from each cluster.
\end{description}

The nested sequence of clusters resulting from hierarchical clustering
methods can be represented graphically using a dendrogram. A
\textit{dendrogram} is a binary tree diagram in which the terminal
nodes, or leaves, represent individual observations and in which the
height of the nodes reflects the dissimilarity between the two
clusters that are joined.  While dendrograms are quite appealing
because of their apparent ease of interpretation, they can be
misleading.

First, the dendrogram corresponding to a given hierarchical clustering
is not unique, since for each merge one needs to specify which subtree
should go on the left and which on the right.  For $n$ observations,
there are $n-1$ merges and hence $2^{(n-1)}$ possible orderings of the
terminal nodes. Therefore, $2^{(n-1)}$ different dendrograms are
consistent with any given sequence of hierarchical clustering
operations. The ordering is of practical importance because it will
result in different genes being placed next to each other in the
heat-maps and possibly different interpretations of the same data. A
number of heuristics have been suggested for ordering the terminal
nodes in a dendrogram.  The default in the R function \texttt{hclust}
from the \textit{cluster} package is to order the subtrees so that the
tighter cluster is on the left (i.e., the most recent merge of the
left subtree is at a lower value than the last merge of the right
subtree).

A second, and perhaps less recognized shortcoming of dendrograms, is
that they \textit{impose} structure on that data, instead of
\textit{revealing} structure in these data. Indeed, dendrograms are
often viewed as graphical summaries of the data, rather than of the
results of the clustering algorithm. Such a representation will be
valid only to the extent that the pairwise dissimilarities possess the
hierarchical structure imposed by the clustering algorithm.  Note that
in particular, the same matrix of pairwise distances between
observations will be represented by a different dendrogram depending
on the distance function that is used to compute between-cluster
distances (e.g., complete or average linkage).  The {\em cophenetic
correlation coefficient} can be used to measure how well the
hierarchical structure from the dendrogram represents the actual
distances. This measure is defined as the correlation between the
$n(n-1)/2$ pairwise dissimilarities between observations and their
{\em cophenetic} dissimilarities from the dendrogram, i.e., the
between cluster dissimilarities at which two observations are first
joined together in the same cluster. \\

Next, we apply three agglomerative hierarchical clustering methods to
cluster tumor samples. For each methods, we plot a dendrogram, compute
the corresponding cophenetic correlation, and compare
the three clusters obtained from cutting the tree into three branches
to the actual tumor labels (ALL B-cell, ALL T-cell, and AML).

<<average,fig=TRUE>>=
hc1 <- hclust(as.dist(d), method="average")
coph1<-cor(cophenetic(hc1),as.dist(d))

plot(hc1,
main=paste("Dendrogram for ALL AML data: Coph = ", round(coph1,2)),
sub="Average linkage, correlation matrix, G=3,051 genes")

cthc1 <- cutree(hc1, 3)
table(Y, cthc1)

@

<<single,fig=TRUE>>=
hc2 <- hclust(as.dist(d), method="single")
coph2<-cor(cophenetic(hc2),as.dist(d))

plot(hc2,
main=paste("Dendrogram for ALL AML data: Coph = ", round(coph2,2)),
sub="Single linkage, correlation matrix, G= 3,051 genes")

cthc2 <- cutree(hc2, 3)
table(Y, cthc2)

@

<<complete,fig=TRUE>>=
hc3 <- hclust(as.dist(d), method="complete")
coph3<-cor(cophenetic(hc3),as.dist(d))

plot(hc3,
main=paste("Dendrogram for ALL AML data: Coph = ", round(coph3,2)),
sub="Complete linkage, correlation matrix, G= 3,051 genes")

cthc3 <- cutree(hc3, 3)
table(Y, cthc3)
@

Divisive hierarchical clustering is available using \texttt{diana} from the
\textit{cluster} package.

<<diana,fig=TRUE>>=
di1 <- diana(as.dist(d))
cophdi<- cor(cophenetic(di1), as.dist(d))

plot(di1, which.plots=2, #$
main=paste("Dendrogram for ALL AML data: Coph = ", round(cophdi,2)),
sub="Divisive algorithm, correlation matrix, G= 3,051 genes")

ct.di <- cutree(di1, 3)
table(Y, ct.di)
@

\section{Partitioning Methods}

Here we apply the {\em partitioning around medoids} (PAM) method to cluster
tumor mRNA samples.

<<PAM3,fig=TRUE>>=
##PAM
set.seed(12345)

pm3 <- pam(as.dist(d), k=3,diss=TRUE)
table(Y, pm3$clustering)

clusplot(d, pm3$clustering, diss=TRUE, labels=3, #$
col.p=1, col.txt=as.integer(factor(Y))+1,
main="Bivariate cluster plot for ALL AML data\n Correlation matrix, K=3, G=3,051 genes")

@
%<<PAM3plot>>=
%plot(pm3,which.plots=1)
%
%@

<<PAM3plot2>>=
plot(pm3,which.plots=2, main="Silhouette plot for ALL-AML Data")
@

One of the more difficult questions that arises when clustering data
is the assessment of how many clusters there are in the data.  We now
consider using PAM with four clusters.

<<PAM4,fig=TRUE>>=
##what happens if we try four?
set.seed(12345)
pm4 <- pam(as.dist(d), k=4, diss=TRUE)

clusplot(d, pm4$clustering, diss=TRUE, labels=3, #$
col.p=1, col.txt=as.integer(factor(Y))+1,
main="Bivariate cluster plot for ALL AML data\n Correlation matrix, K=4, G=3,051 genes")
@


\end{document}
