%\VignetteIndexEntry{protViz primer}

\documentclass[12pt]{article}

\usepackage[utf8]{inputenc}
\usepackage{hyperref}
\usepackage{color}
\renewcommand*\familydefault{\sfdefault}

\definecolor{NoteGrey}{rgb}{0.96,0.97,0.98}

\textwidth=6.2in
\textheight=9.5in
@ \parskip=.3cm
\oddsidemargin=.1in
\evensidemargin=.1in
\headheight=-1in

\newcommand{\scscst}{\scriptscriptstyle}
\newcommand{\scst}{\scriptstyle}
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\code}[1]{{\texttt{#1}}}

\author{Christian Panse \and Jonas Grossmann}
\begin{document}
\title{\code{protViz}: Visualizing and Analyzing Mass Spectrometry Related Data in Proteomics}

\maketitle

\fcolorbox{black}{NoteGrey} {
\begin{minipage}{13.5cm}
\vspace{0.5cm}
\begin{center}
\textbf{ Vignette for package version $\ge v.0.1.48$}\\
\vspace{0.5cm}
Recent changes, updates, and  new features:
\begin{itemize}
\item 0.2.9 added \code{sigmamix} data
\item 0.2.5 export aa2mass
\item 0.2.0 peakplot margin bug fix
\end{itemize}
\end{center}
\vspace{0.5cm}
\end{minipage}
}

\tableofcontents
\newpage

\section{Preliminary Note}

\code{protViz} is an R package to do quality checks, visualizations and analysis of mass spectrometry data, coming from proteomics experiments.
The package is developed, tested and used at the Functional Genomics Center Zurich. 
We use this package mainly for prototyping, teaching, and having \code{fun} with proteomics data.
But it can also be used to do data analysis for small scale data sets. Nevertheless, if one is patient, it also handles large data sets.

\section{Related Work}
{\em The methode of choice in proteomics is mass spectrometry.}
There are already packages in R which deal with mass spec related data.
Some of them are listed here:
\begin{itemize}
\item \code{MSnbase} package (basic functions for mass spec data including quant aspect with iTRAQ data)\\
\url{http://www.bioconductor.org/packages/release/bioc/html/MSnbase.html}
\item \code{plgem} -- spectral counting quantification, applicable to MudPIT experiments\\
\url{http://www.bioconductor.org/packages/release/bioc/html/plgem.html}
\item \code{synapter} -- MSe (Hi3 = Top3 Quantification) for Waters Q-tof data aquired in MSe mode\\ 
\url{http://www.bioconductor.org/packages/release/bioc/html/synapter.html}
\item \code{mzR} \\ 
\url{http://www.bioconductor.org/packages/release/bioc/html/mzR.html}
\item isobar iTRAQ/TMT quantification package\\
\url{http://www.bioconductor.org/packages/release/bioc/html/isobar.html}
\item \code{readMzXmlData}\\ 
\url{http://cran.r-project.org/web/packages/readMzXmlData/}
\item \code{msQC}\\
\url{http://bioconductor.org/packages/3.0/bioc/html/msQC.html}
\end{itemize}

\section{Get Data In -- Preprocessing}
The most time consuming and challenging part for data analysis and visualization is shaping the data that they can easily be processed. In this package, we intentionally left this part away because it is very infrastructure dependent. Moreover we use also commercial tools to analyze data and export the data into R accessible formats. We provide different kind of importers if these formats are available, but with very little effort, one can bring other exports in a similar format which will make it easy to use our package for a variety of tools.

\subsection{Identification - In-silico from Proteins to Peptides}
For demonstration we use a sequence of peptides derived from a tryptics digest using the Swissprot \code{FETUA\_BOVIN Alpha-2-HS-glycoprotein precursor (Fetuin-A) (Asialofetuin)} protein.

\code{fcat} and \code{tryptic-digest} are commandline programs which are included in the package.
\code{fcat} removes the lines starting with \code{>} and all 'new line' character within the protein sequence while \code{tryptic-digest} is doing the triptic digest of a protein sequence applying the rule: cleave after arginine (R) and lysine (K) except followed by proline(P).

\begin{verbatim}
$ cat Fetuin.fasta
MKSFVLLFCLAQLWGCHSIPLDPVAGYKEPACDDPDTEQAALAAVDYINKHLPRGYKHTL
NQIDSVKVWPRRPTGEVYDIEIDTLETTCHVLDPTPLANCSVRQQTQHAVEGDCDIHVLK
QDGQFSVLFTKCDSSPDSAEDVRKLCPDCPLLAPLNDSRVVHAVEVALATFNAESNGSYL
QLVEISRAQFVPLPVSVSVEFAVAATDCIAKEVVDPTKCNLLAEKQYGFCKGSVIQKALG
GEDVRVTCTLFQTQPVIPQPQPDGAEAEAPSAVPDAAGPTPSAAGPPVASVVVGPSVVAV
PLPLHRAHYDLRHTFSGVASVESSSGEAFHVGKTPIVGQPSIPGGPVRLCPGRIRYFKI
\end{verbatim}
%\code{cat Fetuin.fasta | fcat | tryptic-digest}
\begin{verbatim}
$ cat Fetuin.fasta | fcat | tryptic-digest
MK
SFVLLFCLAQLWGCHSIPLDPVAGYK
EPACDDPDTEQAALAAVDYINK
HLPR
GYK
HTLNQIDSVK
VWPR
RPTGEVYDIEIDTLETTCHVLDPTPLANCSVR
QQTQHAVEGDCDIHVLK
QDGQFSVLFTK
CDSSPDSAEDVR
K
LCPDCPLLAPLNDSR
VVHAVEVALATFNAESNGSYLQLVEISR
AQFVPLPVSVSVEFAVAATDCIAK
EVVDPTK
CNLLAEK
QYGFCK
GSVIQK
ALGGEDVR
VTCTLFQTQPVIPQPQPDGAEAEAPSAVPDAAGPTPSAAGPPVASVVVGPSVVAVPLPLHR
AHYDLR
HTFSGVASVESSSGEAFHVGK
TPIVGQPSIPGGPVR
LCPGR
IR
YFK
I
\end{verbatim}


\section{Peptide Identification}
{\em The currency in proteomics are the peptides.}
In proteomics, proteins are digested to so-called peptides since peptides are much easier to handle biochemically than proteins. Proteins are very different in nature some are very sticky while others are soluble in aqueous solutions while again are only sitting in membranes. Therefore, proteins are chopped up into peptides because it is fair to assume, that for each protein, there will be a number of peptides behaving well, so that they can actually be measured with the mass spectrometer. This step introduces another problem, the so-called protein inference problem. In this package here, we do not at all touch upon the protein inference.

\subsection{Computing the Parent Ion Mass}
<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
library(protViz)
op<-par(mfrow=c(1,1))
fetuin<-c('MK', 'SFVLLFCLAQLWGCHSIPLDPVAGYK',
'EPACDDPDTEQAALAAVDYINK',
'HLPR', 'GYK', 'HTLNQIDSVK', 'VWPR', 
'RPTGEVYDIEIDTLETTCHVLDPTPLANCSVR',
'QQTQHAVEGDCDIHVLK', 'QDGQFSVLFTK',
'CDSSPDSAEDVR', 'K', 'LCPDCPLLAPLNDSR',
'VVHAVEVALATFNAESNGSYLQLVEISR',
'AQFVPLPVSVSVEFAVAATDCIAK',
'EVVDPTK', 'CNLLAEK', 'QYGFCK',
'GSVIQK', 'ALGGEDVR', 
'VTCTLFQTQPVIPQPQPDGAEAEAPSAVPDAAGPTPSAAGPPVASVVVGPSVVAVPLPLHR',
'AHYDLR', 'HTFSGVASVESSSGEAFHVGK',
'TPIVGQPSIPGGPVR', 'LCPGR', 'IR', 'YFK', 'I')
(pm<-parentIonMass(fetuin))
op<-par(mfrow=c(2,1))
plot(pm, ylab="peptide mass [in Da]", 
    main="Fetuin Peptide tryptic digested.")
hist(pm, xlab="peptide mass [in Da]")
@

\subsection{In-silico Peptide Fragmentation}
The fragment ions of a peptide can be computed following the rules proposed in \cite{pmid6525415}.
Beside the \code{b} and \code{y} ions the \code{FUN} argument of \code{fragmentIon} defines 
which ions are computed.
the default ions beeing computed are defined in the function \code{defaultIon}.
The are no limits for defining other forms of fragment ions for ETD (c and z ions) CID (b and y ions).
<<keep.source = TRUE, eval = TRUE, fig = FALSE>>=
defaultIon
@

<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
peptides<-c('HTLNQIDSVK', 'ALGGEDVR', 'TPIVGQPSIPGGPVR')

pim<-parentIonMass(peptides)
fi<-fragmentIon(peptides)

par(mfrow=c(3,1)); 
for (i in 1:length(peptides)){
    plot(0,0,
        xlab='m/Z',
        ylab='',
        xlim=range(c(fi[i][[1]]$b,fi[i][[1]]$y)),
        ylim=c(0,1),
        type='n',
        axes=FALSE,
        sub=paste( pim[i], "Da"));
    box()
    axis(1,fi[i][[1]]$b,round(fi[i][[1]]$b,2))
    pepSeq<-strsplit(peptides[i],"")
    axis(3,fi[i][[1]]$b,pepSeq[[1]])

    abline(v=fi[i][[1]]$b, col='red',lwd=2) 
    abline(v=fi[i][[1]]$c, col='orange') 
    abline(v=fi[i][[1]]$y, col='blue',lwd=2)
    abline(v=fi[i][[1]]$z, col='cyan')
}
@

The next lines compute the singly and doubly charged fragment ions of the \code{HTLNQIDSVK} peptide. Which are usually the ones that can be used to make an identification.
<<keep.source = TRUE, eval = TRUE, fig = FALSE>>=
Hydrogen<-1.007825
(fi.HTLNQIDSVK.1<-fragmentIon('HTLNQIDSVK'))[[1]]
(fi.HTLNQIDSVK.2<-(fi.HTLNQIDSVK.1[[1]] + Hydrogen) / 2)
@


\subsection{Peptide Sequence -- Fragment Ion Matching}
Given a peptide sequence and a tandem mass spectrum.
For the assignment of a canditate peptide an in-silico fragment ion
spectra \code{fi}  is computed.
The function \code{findNN} determines for each fragment ion the closesed 
peak in the MS2.
If the difference between the in-silico mass and the measured mass  
is inside the 'accuracy' mass window of the mass spec device the 
in-silico fragment ion is considered as potential hit.

<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
    peptideSequence<-'HTLNQIDSVK'
    spec<-list(scans=1138,
        title="178: (rt=22.3807) [20080816_23_fetuin_160.RAW]",
        rtinseconds=1342.8402,
        charge=2,
        mZ=c(195.139940, 221.211970, 239.251780, 290.221750, 
    316.300770, 333.300050, 352.258420, 448.384360, 466.348830, 
    496.207570, 509.565910, 538.458310, 547.253380, 556.173940, 
    560.358050, 569.122080, 594.435500, 689.536940, 707.624790, 
    803.509240, 804.528220, 822.528020, 891.631250, 909.544400, 
    916.631600, 973.702160, 990.594520, 999.430580, 1008.583600, 
    1017.692500, 1027.605900),
        intensity=c(931.8, 322.5, 5045, 733.9, 588.8, 9186, 604.6,
    1593, 531.8, 520.4, 976.4, 410.5, 2756, 2279, 5819, 2.679e+05,
    1267, 1542, 979.2, 9577, 3283, 9441, 1520, 1310, 1.8e+04,
    587.5, 2685, 671.7, 3734, 8266, 3309))

    fi<-fragmentIon(peptideSequence)
    n<-nchar(peptideSequence)

    by.mZ<-c(fi[[1]]$b, fi[[1]]$y)
    by.label<-c(paste("b",1:n,sep=''), paste("y",n:1,sep=''))

    # should be a R-core function as findInterval!
    idx<-findNN(by.mZ, spec$mZ) 

    mZ.error<-abs(spec$mZ[idx]-by.mZ)

    plot(mZ.error[mZ.error.idx<-order(mZ.error)],
        main="Error Plot",
        pch='o',
        cex=0.5,
        sub='The error cut-off is 0.6Da (grey line).',
        log='y')
    abline(h=0.6,col='grey')
    text(1:length(by.label), 
        mZ.error[mZ.error.idx],  
        by.label[mZ.error.idx],
        cex=0.75,pos=3) 
@

The graphic above is showing the mass error of the assingment between the MS2 \code{spec} and the singly charged fragment ions of \code{HTLNQIDSVK}.
The function \code{psm} is doing the peptide sequence matching. Of course, the more theoretical ions match (up to a small error tolerance, given by the system) the actually measured ion series, the more likely it is, that the measured spectrum indeed is from the inferred peptide (and therefore the protein is identified)
\subsection{Modifications}
<<keep.source = TRUE, eval = TRUE, fig = FALSE>>=
library(protViz)

ptm.0<-cbind(AA="-", 
    mono=0.0, avg=0.0, desc="unmodified", unimodAccID=NA)

ptm.616<-cbind(AA='S',
    mono=-27.010899, avg=NA, desc="Substituition", 
    unimodAccID=616)

ptm.651<-cbind(AA='N',
    mono=27.010899, avg=NA, desc="Substituition", 
    unimodAccID=651)


m<-as.data.frame(rbind(ptm.0, ptm.616, ptm.651))

genMod(c('TAFDEAIAELDTLNEESYK','TAFDEAIAELDTLSEESYK'), m$AA)

fi<-fragmentIon(c('TAFDEAIAELDTLSEESYK', 
    'TAFDEAIAELDTLNEESYK', 'TAFDEAIAELDTLSEESYK', 
    'TAFDEAIAELDTLNEESYK'), 
        modified=c('0000000000000200000', 
        '0000000000000100000', '0000000000000000000', 
        '0000000000000000000'), 
    modification=m$mono)

#bh<-c('TAFDEAIAELDTLNEESYK', 'TAFDEAIAELDTLSEESYK')
#fi<-fragmentIon(rep('HTLNQIDSVK',2), 
#    modified=c('0000000100','0000000000'), 
#    modification=m[,2])
@


\subsection{Labeling Peaklists}
The labeling of the spectra can be done with the \code{peakplot} function.
%\SweaveOpts{width=9.5, height=6}
%\setkeys{Gin}{width=0.95\columnwidth}
<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
data(msms)
op<-par(mfrow=c(2,1))
peakplot("TAFDEAIAELDTLNEESYK", msms[[1]])
peakplot("TAFDEAIAELDTLSEESYK", msms[[2]])
par(op)
@
%\subsection{Finding Glyco Patterns}

\section{Quantification}
For an overview on Quantitative Proteomics read \cite{pmid22772140, pmid22821268}.
The authors are aware that meaningful statistics usually require much higher number
of biological replicates. 
In almost all cases there are not more than three to six repitions.
For the moment there are limited options due to the availabilty of machine time
and the limits of the technologies.

\subsection{Label-free methods on protein level}
The data set \code{fetuinLFQ} contains a subset of our results descriped in \cite{pmid20576481}.
The example below shows a visualization using trellis plots. It graphs 
the abundance of four protein in dependency from the 
fetuin concentration spiked into the sample.

<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
library(lattice)
data(fetuinLFQ)

cv<-1-1:7/10
t<-trellis.par.get("strip.background")
t$col<-(rgb(cv,cv,cv))
trellis.par.set("strip.background",t)

print(xyplot(abundance~conc|prot*method,
    groups=prot,
    xlab="Fetuin concentration spiked into experiment [fmol]",
    ylab="Abundance",
    aspect=1,
    data=fetuinLFQ$t3pq[fetuinLFQ$t3pq$prot 
        %in% c('Fetuin', 'P15891', 'P32324', 'P34730'),],
    panel = function(x, y, subscripts, groups) {
        if (groups[subscripts][1] == "Fetuin")  {
            panel.fill(col="#ffcccc")
        }
        panel.grid(h=-1,v=-1)
        panel.xyplot(x, y)
        panel.loess(x,y, span=1)
        if (groups[subscripts][1] == "Fetuin")  {
            panel.text(min(fetuinLFQ$t3pq$conc),
                max(fetuinLFQ$t3pq$abundance),
                paste("R-squared:", 
                round(summary(lm(x~y))$r.squared,2)),
                cex=0.75,
                pos=4)
        }
    }
))
@

The plot shows the estimated concentration of the four proteins using the top three most intense peptides. The Fetuin peptides are spiked in with increasing concentration while the three other yeast proteins are kept stable in the background. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{pgLFQ -- LCMS based label-free quantification} 

LCMS based label-free quantification is a very popular method to extract relative quantitative information from mass spectrometry experiments.
At the FGCZ we use the software ProgenesisLCMS for this workflow \url{http://www.nonlinear.com/products/progenesis/lc-ms/overview/}.
Progenesis is a graphical software which does the aligning between several LCMS experiments, extracts signal intensities from LCMS maps and annotates the mastermap with peptide and protein labels.

<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
data(pgLFQfeature)
data(pgLFQprot)

featureDensityPlot<-function(data, n=ncol(data), nbins=30){
    my.col<-rainbow(n);
    mids<-numeric()
    density<-numeric()
    for (i in 1:n) { 
        h<-hist(data[,i],nbins, plot=FALSE)
        mids<-c(mids, h$mids)
        density<-c(density, h$density)
    }
    plot(mids,density, type='n')
    for (i in 1:n) { 
        h<-hist(data[,i],nbins, plot=FALSE)
        lines(h$mids,h$density, col=my.col[i])
    }
    legend("topleft", names(data), cex=0.5,
        text.col=my.col
    )
}

par(mfrow=c(1,1)); 
featureDensityPlot(asinh(pgLFQfeature$"Normalized abundance"),
    nbins=25)
@

The \code{featureDensityPlot} shows the normalized signal intensity 
distribution (asinh transformed) over 24 LCMS runs which are aligned in this experiment.

<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
op<-par(mfrow=c(1,1),mar=c(18,18,4,1),cex=0.5)
samples<-names(pgLFQfeature$"Normalized abundance")
image(cor(asinh(pgLFQfeature$"Normalized abundance")),
    col=gray(seq(0,1,length=20)),
    main='pgLFQfeature correlation',
    axes=FALSE)

axis(1,at=seq(from=0, to=1, 
    length.out=length(samples)), 
    labels=samples, las=2)

axis(2,at=seq(from=0, to=1, 
    length.out=length(samples)), labels=samples, las=2)
par(op)
@

This image plot shows the correlation between runs on feature level (values are asinh transformed). White is perfect correlation while black indicates a poor correlation.

<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
op<-par(mfrow=c(1,1),mar=c(18,18,4,1),cex=0.5)
image(cor(asinh(pgLFQprot$"Normalized abundance")),
    main='pgLFQprot correlation',
    axes=FALSE,
    col=gray(seq(0,1,length=20)))
axis(1,at=seq(from=0, to=1, 
    length.out=length(samples)), labels=samples, las=2)
axis(2,at=seq(from=0, to=1, 
    length.out=length(samples)), labels=samples, las=2)
par(op)
@

This figure shows the correlation between runs on protein level (values are asinh transformed). 
White is perfect correlation while black indicates a poor correlation.
Stricking is the fact that the six biological replicates for each condition cluster very well.

<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
par(mfrow=c(2,2),mar=c(6,3,4,1))
ANOVA<-pgLFQaov(pgLFQprot$"Normalized abundance", 
    groups=as.factor(pgLFQprot$grouping), 
    names=pgLFQprot$output$Accession,
    idx=c(15,16,196,107),
    plot=TRUE)
@

This figure shows the result for four proteins which either differ significantly in expression accross conditions (green boxplots) using an analysis of variance test, or non differing protein expression (red boxplot).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{iTRAQ -- Two Group Analysis} 

The data for the next section is an iTRAQ-8-plex experiment where two conditions are compared (each condition has 4 biological replicates)

\subsubsection{Sanity Check}
<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
data(iTRAQ)
x<-rnorm(100)
par(mfrow=c(3,3),mar=c(6,4,3,0.5));
for (i in 3:10){
    qqnorm(asinh(iTRAQ[,i]), 
        main=names(iTRAQ)[i])
    qqline(asinh(iTRAQ[,i]), col='grey')
}
b<-boxplot(asinh(iTRAQ[,c(3:10)]), main='boxplot iTRAQ')
@

A first quality check to see if all reporter ion channels are having the same distributions. Shown in the figure are Q-Q plots of the individual reporter channels against a normal distribution.
The last is a boxplot for all individual channels.

\subsubsection{On Protein Level}
<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
data(iTRAQ)
group1Protein<-numeric()
group2Protein<-numeric()

for (i in c(3,4,5,6))
    group1Protein<-cbind(group1Protein,
        asinh(tapply(iTRAQ[,i], paste(iTRAQ$prot), sum, na.rm=TRUE)))
         
for (i in 7:10)
    group2Protein<-cbind(group2Protein,
        asinh(tapply(iTRAQ[,i], paste(iTRAQ$prot), sum, na.rm=TRUE)))
                  
                  
par(mfrow=c(2,3),mar=c(6,3,4,1))
for (i in 1:nrow(group1Protein)){
    boxplot.color="#ffcccc"
    tt.p_value<-t.test(as.numeric(group1Protein[i,]), 
        as.numeric(group2Protein[i,]))$p.value       

    if (tt.p_value < 0.05)
        boxplot.color='lightgreen'

    b<-boxplot(as.numeric(group1Protein[i,]), 
        as.numeric(group2Protein[i,]),
        main=row.names(group1Protein)[i],
        sub=paste("t-Test: p-value =", round(tt.p_value,2)),
        col=boxplot.color,
        axes=FALSE)
    axis(1, 1:2, c('group_1','group_2')); axis(2); box()

    points(rep(1,b$n[1]), as.numeric(group1Protein[i,]), col='blue')
    points(rep(2,b$n[2]), as.numeric(group2Protein[i,]), col='blue')
}
@

This figure shows five proteins which are tested if they differ accross conditions using the four biological replicates with a t-test.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{On Peptide Level}

The same can be done on peptide level using the \code{protViz} function \code{iTRAQ2GroupAnalysis}. 
<<keep.source = TRUE, eval = TRUE, fig = FALSE>>=
data(iTRAQ)
q<-iTRAQ2GroupAnalysis(data=iTRAQ, 
    group1=c(3,4,5,6), 
    group2=7:10, 
    INDEX=paste(iTRAQ$prot,iTRAQ$peptide), 
    plot=FALSE)
q[1:10,]
@


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Pressure Profiles QC}

A common problem with mass spec setup is the pure reliability of the high pressure pump.
The following graphics provide visualizations for quality control.

On overview of the pressure profile data can be seen by using the \code{ppp} function.
<<keep.source = TRUE, eval = TRUE, fig = FALSE>>=
data(pressureProfile)
ppp(pressureProfile)
@

The lines plots the pressure profiles data on a scatter plot 'Pc' versus
'time' grouped by time range (no figure because of too many data items).

The Trellis \code{xyplot} shows the Pc development over each instrument 
run to a specified relative run time (25,30,...).

<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
pp.data<-pps(pressureProfile, time=seq(25,40,by=5))
print(xyplot(Pc ~ as.factor(file) | paste("time =", 
    as.character(time), "minutes"),
    panel = function(x, y){
        m<-sum(y)/length(y)
        m5<-(max(y)-min(y))*0.05
        panel.abline(h=c(m-m5,m,m+m5), 
            col=rep("#ffcccc",3),lwd=c(1,2,1))
        panel.grid(h=-1, v=0)
        panel.xyplot(x, y)
    },
    ylab='Pc [psi]',
    layout=c(1,4),
    sub='The three red lines indicate the average plus min 5%.',
    scales = list(x = list(rot = 45)),
    data=pp.data))
@

While each panel in the \code{xyplot} above shows the data to a given point in time, we try to use the levelplot 
to get an overview of the whole pressure profile data.
<<keep.source = TRUE, eval = TRUE, fig = TRUE>>=
pp.data<-pps(pressureProfile, time=seq(0,140,length=128))
print(levelplot(Pc ~ time * as.factor(file),
    main='Pc(psi)',
    data=pp.data,
    col.regions=rainbow(100)[1:80]))
@

\label{lab:datapreparation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Session information}\label{sec:sessionInfo}

<<keep.source = TRUE, eval = TRUE, fig = FALSE>>=
sessionInfo()
packageDescription('protViz')
@

\bibliography{protViz}
\bibliographystyle{plain}

\end{document}
