**RProteomics Package Overview**

**Rich Haney**

**19 May 2006**

**Introduction**

This package involves using simple spectrum models and pipelines to
integrate low-level processing, peak finding, alignment and other
feature detection tasks, together with methods for high-level
classification and and other tasks. In the approach, we decompose
spectra in different ways. For example, we might want to decompose a
spectrum ( I ) into a set of background features, noise and a (
remaining ) signal according the following idea:

> I \~ runningQuantile(I,1001,percentile=70) +
>
> waveletNoise (I) +
>
> ISignal

Within the above, simple additive model, we solve for ISignal. In terms
of handling that, we define a pipeline of tasks to compute over a set of
just the first two spectra using the following:

p("

> tmp1=runningQuantile(I\[,1:2\],1001,percentile=70)
>
> tmp1a=minus(I\[,1:2\],tmp1\[,1:2\])
>
> tmp2=wavelet(tmp1a\[,1:2\])
>
> ISignal=minus(tmp1a\[,1:2\],tmp2\[,1:2\])
>
> plot2D(I\[,1:2\],ISignal\[,1:2\])

")

In the above expressions, note that "I" is for us a matrix of
intensities. In that matrix, columns of I comprise different scans, rows
of I comprise intensities, and a common vector of m/z values defines the
common m/z values of all rows.

In order to obtain such a matrix, we often interpolate among the
\<m/z,I\> data. For example, the following does spline-based
interpolation:

p("

I = interpolate(rawData,method='spline',stepSize=.5)

")

On the other hand, this does interpolation along a quadratic curve of
m/z intensity values that has been supplied:

p("

I = interpolate(rawData,method='linear',quadraticCurve)

")

Here "I" is apt to be too big too fit into memory, so it is apt to
represent a matrix of intensities stored on disk. On the other hand, the
quadratic curve can typically fit in memory. So, in a model such as
that, "rawData" and "I" are apt to both represent sets of files on disk,
while "quadraticCurve" represents data in memory.

To be sure, here are many methods of interpolation and approximation. So
we also allow for basic linear interpolation, splines, loess-based
interpolation, and a "sliding-loess" based method.

It is after the interpolation step that basic modeling steps are used.
In the first example above we used a nonparametric model for the
background of the spectrum. That example comprised use of a "running
Quantile" of the 70^th^ percentile of the spectrum. However, we might
like to segment that baseline, so that it is not modeled in that single,
non-parametric way. In this context, in the following case, we define
the background as comprising three segments altogether, including a
logistic fit of the spectrum from m/z values 2000 -- 2149, a quadratic
fit of values from 2150 -- 5000, and a running 70^th^ percentile of the
rest of data from 5001 to 8000:

I \~ logistic (I\[2000:2149,\]) +

fit1 (I\[2150:5000,\]) +

runningQuantile(I\[5001:80000,\],1001,percentile=70) +

waveletNoise (I\[2000:80000,\]) +

ISignal

Using a slightly more compressed syntax than the earlier example, we
might solve this using the following:

colStart \<- 1; colEnd \<- 50

p("

tmp1=logistic(I\[2000:2149,colStart:colEnd\], subtract=TRUE)

> tmp2=fit1(tmp1\[2150:5000, colStart:colEnd\], subtract=TRUE)
>
> tmp3=runningQuantile(tmp2\[5001:80000,colStart:colEnd\],1001,percentile=70,subtract=TRUE)
>
> ISignal=wavelet(tmp3\[2000:80000,colStart:colEnd\],subtract=TRUE)

plot2D(I\[,1:2\],ISignal\[,1:2\])

")

Having obtained a signal, we are then apt to find, align and quantify
peaks ( or handle peaks in some other order of processing.

This comprises a brief overall look as to what we are trying to do.

Next is a section on the overall ( relatively short ) list of functions
actually "exposed" to users in the application sense.

**Overall Routines directly "Explosed" at the user level.**

As my be gleaned from the above, the main exported command is:

1.  p ( for run pipeline)

Minor and supporting ones are:

2.  dataset( to define a dataset ), variable ( to define a variable ),
    and variables() to list or reset variables

3.  toMemory and toDisk ( to handle and covert data in terms of files)

4.  setOption -- to set one of 3 overall system defaults

5.  A set of routines for testing the installation ( e.g., runBGTests )

For the most part, the system only requires the use of two commands,
that is, "p" and "dataset".

The following summarizes the various ( exposed ) routines:

1.  **p ( for run pipeline)**

To help it communicate with Java, allow the package to make use of data
files stored on disk, and in general to run pipelines and so on, this
package makes use of a thin "services" or "pipeline" layer above R.

Usage:

p ( command )

Here "command" is a character string that defines a thin layer over R
functions that can then ( and with the help of some Java components ) be
gridified a little more easily.

Example, with I and I2 stored as sets of binary files on disk:

colStart \<- 1; colEnd \<- 10

p( "

I2 = highPassW(I\[,colStart:colEnd\],2048,numCoeffsToDrop=500)

plot2D(I2\[,colStart:colEnd\]')

" )

Notice that what is passed in as a string. I this particular case, the
variables "I" and "I2" are stored in the form of binary files on disk.

One can define a string of commands separated by semicolons(";"), or, as
in the above case as commands that are simply on separate lines.

The current order of processing is to look for data files on disk; if
they are not found, the system checks to see if there are R variables
known to the procedure. In this case, the variables colStart, colEnd and
the overall expression colStart:colEnd are variables in memory.

For the most part, it is not necessary to use other commands aside from
the dataset command.

2.  **dataset ( to define explicit location and default storage types of
    datasets ); variable and variables() to define storage types for
    individual variables**

Usage:

dataset(name, rootDirectory, storageType)

Examples:

dataset('patz2',rootDirectory='/data')

dataset('patz2',rootDirectory='c:/dataFiles', storageType='binaryFiles')

As noted again, in order for each step to be part of a pipeline, and
communicate with a gridified services layer, variables by default are
stored externally to the application. They can be in memory, take the
form of sets of binary files, single binary files ( primitive in this
version), single CSV files with multiple columns, or sets of CSV columns
with one or two columns each. ( Lamentably, although the system can read
in mzXML files and use databases, the use of terms with default storage
types of mzXML or databases is not supported in this version.)

Some times it is helpful to define variables in models to have a
specific storage type.

**variable** -- defines a variable and its storage type

Usage:

variable('I','patz2',storageType='csvFiles'

variable('I','storageType='memory')

**variables** -- list ( or clear out ) a current list of variables in
use by the system.

Usage:

variables()

variables(reset=TRUE)

Examples:

variables() simply returns the list of all variables in the system,
together with their associated datasets, storage types, location on the
system, and other attributes.

3.  **toMemory and toDisk -- moving data back and forth between disk and
    memory**

Usage:

variableInMemory \<- toMemory('variable')

toDisk('variableNameOnDisk', variableInMemory )

Examples ( together with other examples of other conversion functions )

colStart \<- 1; colEnd \<- 10

variable(\'variableStoredInOneFile\',storageType=\'csvFile\');

p(\'oneFile = I\[,colStart:colEnd\]\' )

variableStoredInMemory \<- toMemory(\'oneFile\[,1:2\]\');

toDisk('variableStoredInbinaryFiles1\',variableStoredInMemory)

variable('variableStoredInCsvFiles\',storageType=\'csvFiles\');

p('variableStoredInCsvFiles = oneFile\[,colStart:colEnd\]\')

variable('variableStoredInCsvFiles1\',storageType=\'csvFiles\');

p('variableStoredInCsvFiles12 = variableStoredInCsvFiles \[,1:2\]\')

4.  **setOption -- set a list of options**

There are just 3 current options in the system, each of which can take
the values of TRUE or FALSE

-   batchMode -- whether a plot command is expecting a user to zoom and
    pan a given plot.

-   nonNegative -- whether a given denoising command returns only
    positive values

-   maxIsSignal -- whether a given denoising command clips each
    potential value as being the minimum of the denoised or smoothed
    value, and the original signal value

Usage:

setOption(list)

Examples:

setOption(list(batchMode=TRUE))

#The following is the current system default:

setOption(list(batchMode=FALSE, nonnegative=TRUE, maxIsSignal=TRUE))

5.  **Set of procedures to test installations.**

The following tests instatllation:

runInterpolateTests();runBGTests();runNormalizationTests();

runDenoiseTests();runHighLevelTests();

**Examples of Simple Pipelines**

We can think of procedures as involving

-   Interpolation and approximation

-   Data normalization

-   Background removal

-   Denoising/smoothing

-   Peak finding and peak alignment

-   High-level analysis tasks ( classification and quantitation )

-   Evaluating the results

In this discussion, we give examples of each of the areas.

**Interpolation/approximation**

The following represent some examples of data normalization

colStart \<- 1; colEnd \<- 2

setOption(list(batchMode=TRUE))

#Do basic interpolation of all values, using the default of cubic spline
interpolation, and then plot the results;

p( \' IN= interpolate(twoColumns\[,colStart:colEnd\]);

plot2D(IN\[,colStart:colEnd\])\' )

#In this case explicitly give the method as spline-based interpolation.

p(\" IN =interpolate(twoColumns\[,colStart:colEnd\],method=\'spline\')

plot2D(IN\[,colStart:colEnd\])\" )

#In this case specify the method of linear approximation

p( \" IN = interpolate(twoColumns\[,colStart:colEnd\],method=\'linear\'
");

#In this case define an explicit grid of values to be supplied.

#Note that the grid could be a quadratic curve or some other
representation.

#Grids can be supplied for all except the "spline" method. That is, you
can pass a vector ( or grid ) of values.

grid \<- 2000:80000.

p ( \'IN = interpolate(twoColumns\[,colStart:colEnd\],grid)');

#In this case define a step size of .5, instead of the default of 1.

#A step size can be provided for all methods ( spline, loess,
slidingLoess, linear). ( At this point there is a limitation in that you
cannot pass a grid of values into the spline-based function.

p ( 'IN = interpolate(twoColumns\[,colStart:colEnd\],stepSize=.5)');

#In this case, define sliding-loess based interpolation, with various
default values. Here 301 is the window width.

p(
\'IL2=interpolate(twoColumns\[,1:1\],method=\"slidingLoess\",width=301,span=.92,degree=2)\');

#In this case, just grab the second column of two-column files.

#This makes sense if all mz values of all spectra are the same.

p ( \'justI = secondColumn(twoColumns\[,colStart:colEnd\])');

#In this case, we are picking out the mz values, that is, the first
columns

#of all of the data.

p( \'justMz = firstColumn(twoColumns\[,colStart:colEnd\])');

setOption(list(batchMode=FALSE))

**Examples of Normalization**

The following are some ( hopefully useful ) examples of data
normalization

colStart \<- 1; colEnd \<- 10; setOption(list(batchMode=TRUE))

#Find the cuberoot of values

p(\' IcubeRoot = cubeRoot(I\[,colStart:colEnd\]);
plot2D(I\[,colStart:colEnd\],IcubeRoot\[,colStart:colEnd\]\')

#Find the square of values in the spectrum

p(\' Isquare = square(I\[,colStart:colEnd\]);

plot2D(I\[,colStart:colEnd\],Isquare\[,colStart:colEnd\]\')

#Find square root

p(\' Isqrt = sqrt(I\[,colStart:colEnd\];

plot2D(I\[,colStart:colEnd\],Isqrt\[,colStart:colEnd\]\')

#Find the natural log

p(\' Ilog = log(I\[,colStart:colEnd\]);

plot2D(I\[,colStart:colEnd\],Ilog\[,colStart:colEnd\]\')

#Find the log ( base 2 )

p(\' Ilog2 = log2(I\[,colStart:colEnd\]);
plot2D(I\[,colStart:colEnd\],Ilog2\[,colStart:colEnd\]\')

p(\'Ilog10 = log10(I\[,colStart:colEnd\]);
plot2D(I\[,colStart:colEnd\],Ilog10\[,colStart:colEnd\]\')

#Find the exponent ( power

p(" Ipower = power(I\[,colStart:colEnd\],exponent=.4);

plot2D(I\[,colStart:colEnd\],Ipower\[,colStart:colEnd\]")

#Find the square root as the power = .5

p(" Ipower = power(I\[,colStart:colEnd\],exponent=1.5);
plot2D(I\[,colStart:colEnd\])")

#Find the exponent

exponent \<- 1.5

p(" I2 = power(I\[, colStart:colEnd\],exponent=exponent)

plot2D(I\[,colStart:colEnd\],I2\[,colStart:colEnd\]\")

#normalize by the maximum value

p(" normalizeByMax = normalizeByMax(I\[,colStart:colEnd\]);
plot2D(I\[,colStart:colEnd\],normalizeByMax\[,colStart:colEnd\])")

#Normalize by the maximum value

p("

normalizeByMax=NrmalizeByMax(I\[,colStart:colEnd\],rowStart=2500,rowEnd=10000)

plot2D(normalizeByMax\[,colStart:colEnd\])")

p(" normalizeUsingMinAndMax =
normalizeUsingMinAndMax(I\[,colStart:colEnd\]);
plot2D(normalizeUsingMinAndMax\[,colStart:colEnd\]")

p(" normalizeUsingMinAndMax =
normalizeUsingMinAndMax(I\[,colStart:colEnd\], rowStart=2000 , rowEnd =
6000)

plot2D(normalizeUsingMinAndMax\[,colStart:colEnd\]")

p(" normalizeByMax = normalizeByMax(I\[,colStart:colEnd\], rowStart=2000
, rowEnd = 6000))

plot2D(normalizeByMax\[,colStart:colEnd\] ")

p(" IOCNormalize = IOCNormalize(I\[,colStart:colEnd\]);
plot2D(I\[,colStart:colEnd\],IOCNormalize\[,colStart:colEnd\]) ")

p(" IOCNormalize =
IOCNormalize(I\[,colStart:colEnd\],rowStart=2500,rowEnd=8000);
plot2D(I\[,colStart:colEnd\],IOCNormalize\[,colStart:colEnd\]) ")

p(" iocnorm = IOCNormalize(I\[,colStart:colEnd\], rowStart=2000 , rowEnd
= 6000)

plot2D(Iiocnorm\[,colStart:colEnd\])\")

**Examples of Background Determination and Background Correction**

setOption(list(batchMode=TRUE))

colStart \<- 1; colEnd \<- 10

#Determine the running 70^th^ percentile of the data over a window of
length 1025 ( which we may then want to subtract from the data )

p(' I2 = runningQuantile(I\[,colStart:colEnd\],1025,percentile=70)

plot2D(I\[,colStart:colEnd\],I2\[,colStart:colEnd\])

I3 = minus(I2\[,colStart:colEnd\],I\[,colStart:colEnd\])

plot2D(I\[,colStart:colEnd\],I3\[,colStart:colEnd\])

\')

#Determine the running 10th percentile of the data over a window of
length 1025 ( which we may then want to subtract from the data )

p('

> I2=runningQuantile(I\[,colStart:colEnd\],1025,percentile=10,subtract=TRUE)

plot2D(I\[,colStart:colEnd\],I2\[,colStart:colEnd\])

\')

#Determine the running 80^th^ percentile, over a ( short ) window of
length 81 and subtract it from the data.

p(' I2=
runningQuantile(I\[,colStart:colEnd\],81,percentile=80,subtract=TRUE)

plot2D(I\[,colStart:colEnd\],I2\[,colStart:colEnd\])

\')

#Determine a linear fit to a very short part of the spectrum ( mz 2000
-- 2419 ) and just plot that small part of the spectrum

p(' I4 = fit1(I\[2000:2149,colStart:colEnd\])

plot2D(I4\[2000:2149,colStart:colEnd\])\')

#Determine a quadratic fit to a very short part of the spectrum ( mz
2000 -- 2419 ) and just plot that small part of the spectrum

p(' I4 = quadfit(I\[2000:2149,colStart:colEnd\]

plot2D(I\[,colStart:colEnd\],I4\[,colStart:colEnd\])\')

#Determine a quadratic fit to a very short part of the spectrum ( mz
2000 -- 2419 ) and just plot that small part of the spectrum using a
back-handed syntax

p(' I4=
fit1(I\[2000:2149,colStart:colEnd\],quadratic=TRUE,exponential=FALSE)

plot2D(I\[,colStart:colEnd\],I4\[,colStart:colEnd\])\')

#Determine an exponential fit to another section of the spectrum ( mz
5000 -- 20000 )

p(' I4 = expfit(I\[5000:20000,colStart:colEnd\])
plot2D(I\[,colStart:colEnd\],I4\[,colStart:colEnd\])\')

#Determine an exponential fit to another section of the spectrum ( mz
5000 -- 20000 ) with some additional ( redundant ) options )

p('
I4=expfit(I\[5000:20000,colStart:colEnd\],quadratic=FALSE,logistic=FALSE)\')

#Determine a logistic fit to the first part of the spectrum, a linear
fit to a successive piece, the 70^th^ ( running ) quantile for the next
portion; then add up the individual pieces and plot them

p(\' I1 = logistic(I\[2000:2149,colStart:colEnd\],subtract=TRUE)

I2 = fit1(I\[2150:5000,colStart:colEnd\],subtract=TRUE)

I3=runningQuantile(I\[5001:80000,colStart:colEnd\],1001,percentile=70,subtract=TRUE)

I4=add(I1\[2000:2149,colStart:colEnd\],I2\[2150:5000,colStart:colEnd\],I3\[5001:80000,colStart:colEnd\])

plot2D(I\[2000:80000,colStart:colEnd\],I4\[2000:80000,colStart:colEnd\])

')

#Determine a logistic fit to the first part of the spectrum, a linear
fit to a successive piece, the 70^th^ ( running ) quantile for the next
portion; then add up the individual pieces and plot them

p(\' I \~ logistic(I\[2000:2149,colStart:colEnd\]) +
fit1(I\[2150:5000,colStart:colEnd\]) +
runningQuantile(I\[5001:80000,colStart:colEnd\],1001,percentile=70) +
rWaveletNoise(I\[,colStart:colEnd\],501) +
runningQuantile(I\[2000:80000,colStart:colEnd\],1001,percentile=20) +
Signal\'

plot2D(I\[2000:80000,colStart:colEnd\])\')

#p(\'I \~ logistic(I\[2000:2149,colStart:colEnd\]) +
fit1(I\[2150:5000,colStart:colEnd\]) +
runningQuantile(I\[5001:80000,colStart:colEnd\],1001,percentile=70) +
rWaveletNoise(I\[,colStart:colEnd\],501) +
runningQuantile(I\[,colStart:colEnd\],1001,percentile=20) + Signal\')

#titles \<- \'logisticModel(I) + runningQuantile(I,1001,quantile=.7) +
rWaveletNoise(I) + signal(I)\'

titles \<- \'logisticModel(I) + rq(I,1001,quantile=.7) +
rWaveletNoise(I) + signal(I)\'

command \<-
paste(\'plot2D(I\[2000:80000,colStart:colEnd\],Signal\[2000:80000,colStart:colEnd\]\',\',
titles=\',\"\'\",titles,\"\'\", \")\", sep=\"\")

p(command)

p(\'I \~
runningQuantile(I\[5001:80000,colStart:colEnd\],1001,percentile=70) +

> rWaveletNoise(I\[,colStart:colEnd\],501) +
> runningQuantile(I\[,colStart:colEnd\],1001,percentile=20) + Signal\')

titles \<- \'runningQuantile(I,1001,quantile=.7) +

rWaveletNoise(I) + signal(I)\'

command \<-
paste(\'plot2D(I\[2000:80000,colStart:colEnd\],Signal\[2000:80000,colStart:colEnd\]\',\',
titles=\',\"\'\",titles,\"\'\", \")\", sep=\"\")

p(command)

p(\'I \~ logistic(I\[2000:2149,colStart:colEnd\]) +

> fit1(I\[2150:5000,colStart:colEnd\]) +
> runningQuantile(I\[5001:80000,colStart:colEnd\],1001,percentile=70) +
> rWaveletNoise(I\[,colStart:colEnd\],501) +
> runningQuantile(I\[,colStart:colEnd\],1001,percentile=20) + Signal\')

titles \<- \'logisticModel(I) + runningQuantile(I,1001,quantile=.7) +
WaveletNoise(I) + signal(I)\'

> command \<-
> paste(\'plot2D(I\[2000:80000,colStart:colEnd\],Signal\[2000:80000,colStart:colEnd\]\',\',
> titles=\',\"\'\",titles,\"\'\", \")\", sep=\"\")

p(command)

p('
I4=logistic(I\[2000:2149,colStart:colEnd\],logistic=TRUE,quadratic=FALSE,exponential=FALSE)\')

p(' I3=minus(I\[,colStart:colEnd\],I2\[,colStart:colEnd\])\')

p(' I \~fit1(I\[2000:2149,colStart:colEnd\]) +
fit1(I\[2150:5000,colStart:colEnd\]) +

expfit(I\[5001:20000,colStart:colEnd\]) +

Signal

plot2D(I\[,colStart:colEnd\],Signal\[,colStart:colEnd\])\')

**Examples of Denoising**

#Use a high-pass filter with a sliding window of length 2048, dropping
about 1/4th ( 500 of 2048 ) of the ( Fourier ) coefficients

p( \'I2 = highPassW(I\[,colStart:colEnd\],2048,numCoeffsToDrop=500)\')

#Use a high-pass filter with a sliding window length that is much
shorter, and drop out about half of the information.

p( \'I2 = highPassW(I\[,colStart:colEnd\],1024,numCoeffsToDrop=500)\')

#use a sliding-window, high pass filter that is ( again ) shorter than
the previous one, dropping out about 1/3 of the( frequency-domain)
information.

p( \'I2 = highPassW(I\[,colStart:colEnd\],512,numCoeffsToDrop=200)\')

#Use a high pass filter that is covers the entire 78K spectrum, just
ignoring the 500 ( highest ) frequency-domain coefficients.

p( \'I2 = highPass(I\[,colStart:colEnd\],numCoeffsToDrop=500)\')

#Use an undecimated wavelet transform over a set of windows. Here we
could

#also use the term "udwtW" for undecimated wavelet transform ( windowed
)

p( \'IW =
rWavelet(I\[,colStart:colEnd\],1024,threshold.multiplier=10)\')

#Use a decimated wavelet transform over the entire spectrum. Here we
could also use the term "udwt"( for undecimated wavelet transform )
instead of wavelet.

p( \'IW = wavelet(I\[,colStart:colEnd\],threshold.multiplier=10)\')

#Use a sliding-window, loess-based smoother, with an explicit step of
subtracting this from the spectrum

p( \"I2 = lo(I\[,colStart:colEnd\],1024,span=0.1);

I2rem = minus(I\[,colStart:colEnd \], I2\[,colStart:colEnd\])\")

#Use a sliding-loess filter, together with an implicit step of
subtracting this from the spectrum.

p( "I2a = lo(I2Rem\[,colStart:colEnd\],101,span=0.25,subtract=TRUE)

plot2D(I\[,colStart:colEnd\],I2Rem\[,colStart:colEnd\],I2a\[,colStart:colEnd\])

\")

#A set of miscellaneous examples

p( \"I2 = runningQuantile(I\[,colStart:colEnd\],101,percentile=10)\")

p( \'I2 = highPassW(I\[,colStart:colEnd\],512,numCoeffsToDrop=200)\')

p( \'I2 = highPass(I\[,colStart:colEnd\],numCoeffsToDrop=500)\')

p( \'IW = udwtW(I\[,colStart:colEnd\],1024,threshold.multiplier=10)\')

p( \'IW = udwt(I\[,colStart:colEnd\],threshold.multiplier=10)\')

#An example using a loess-based ( moving window ) method, in we use a
specific, manual step to do the subtraction steps

p(\" I2 = lo(I\[,colStart:colEnd\],1024,span=0.1);

I2rem = minus(I\[,colStart:colEnd \], I2\[,colStart:colEnd\])

I2a = lo(I2Rem\[,colStart:colEnd\],101,span=0.25)

I2aRemA = minus(I2rem\[,colStart:colEnd \], I2a\[,colStart:colEnd\]

plot2D(I\[,colStart:colEnd\],I2Rem\[,colStart:colEnd\],I2aRemA\[,colStart:colEnd\])

\")

**An example involving Peak Alignment**

sampleColumns \<- getSampleColumns(\"I\",25)

variable(c(\'alignedData\',\'alignedDataNormal\',\'alignedDataCancerous\'),storageType=\'binaryFiles\');variable(c(\'reducedData\',\'colMeans\',\'eigenVectors\',\'allPeaks1\',\'peakMatrix\'),storageType=\'csvFile\');variable(c(\'discriminant\',\'predicted\'),storageType=\'Rdata\');

#Do the basic alignment step

p(\"

\[allPeaks1,alignedData,mzValues\]=align(I\[,sampleColumns\$allColumns\],iterations=12,tolerance=8,minFractionPeaksInCommon=0.1)

#Determine peak quantities of aligned peaks

peakMatrix=peakQuantities(alignedData\[,sampleColumns\$allColumns\],allPeaks1)

\[reducedData,eigenValues,eigenVectors,colMeans\]=q5PCAStep(peakMatrix,numCoeffsToDrop=3)

#Using the discriminant found, do the classification step

discriminant = lda(reducedData,sampleColumns\$numericLabels)

#predicted =
q5.classifyData(peakMatrix,discriminant,colMeans,eigenVectors)

#Also just plot out the aligned data.

alignedDataNormal = alignedData\[,sampleColumns\$normalColumns\]

alignedDataCancerous = alignedData\[,sampleColumns\$advancedColumns\]

plot2D(alignedDataNormal,alignedDataCancerous)

\")

**Example of Classification**

A sample example involving a PCA-based method using results of a peak
finder is as follows:

sampleColumns \<- getSampleColumns(\"I\",25);

variable(c(\'allPeaks1\',\'peakMatrix\'),storageType=\'csvFile\');

variable(\'q5Results\',storageType=\'Rdata\');

p(\"

> \[allPeaks1,alignedData,mzValues\]=align(I\[,sampleColumns\$allColumns\],iterations=12,tolerance=8,minFractionPeaksInCommon=0.1)

peakMatrix=peakQuantities(alignedData\[,sampleColumns\$allColumns\],allPeaks1)

q5Results = q5PCAAllStep(peakMatrix,numCoeffsToDrop=3)

\")

**Evaluating the Results**

In general, the approach for decomposing the spectra involves taking
existing ( or new ) R functions, and then automatically generating
pipeline and grid wrappers for them. Doing that, we can make use of a
variety of methods. At each step of a processing pipeline, we aim to
include a parametric method, a semi-parametric method, and a
non-parametric method. We judge the results in three ways:

1.  Beauty-contest methods ( i.e., using plots )

2.  Non-parametric or other measures of model fitness ( e.g.,
    Chi-square, KL, other distance measures)

3.  How helpful ( or not ) the various features are when doing
    high-level tasks such as classification or peptide quantitation.

    **Overall Technical Aspects of this System**

    In the rest of this document, we give an overview of the general
    technical aspects of the system.

    **Summary of System Components**

We use the following components so far:

1.  Methods of interpolation, including linear, spline, and loess-based
    methods.

2.  Arithmetic methods of normalization, such as involving dividing by
    maximum values, areas under curves, and quantile-based
    normalization.

3.  Parametric methods for fitting background curves and peaks:
    exponential, bi-exponential, quadratic and linear fits of whole
    spectra or parts of the spectra.

4.  Methods for extracting noise, using Hilbert space projections,
    primarily signal processing methods using FFTs and Local FFTs.

5.  Methods using Besov spaces: basic wavelet methods using the Rice
    Wavelet Toolbox ( RWT ) library.

6.  Methods for alignment of spectra, at this point using linear, FDA
    and a recursive, linear method ( essentially dynamic time warping )

7.  A few methods of classification and quantitation, namely linear
    discriminant analysis, quadratic discriminant analysis, and (very
    shortly) random forests.

Note that there are other possible methods to take at each and every
step. For example, it is also common to do clustering-based approaches
to spectrum alignment.

Having mentioned the above, would want to note the following:

**Possible more advanced approaches that could be used over time:**

To be sure, it makes a lot of sense to incorporate a set of three types
of more advanced models to this procedure, that is:

1.  Hierarchical Bayesian models

2.  Nonlinear mixed effects models

3.  Models based upon convex programming or convex set models.

    However, our existing models are at this time more modest and
    primarily "one step" in nature.

**The two overall ways used to separate out "statistician" ( or
expert-level options ) from "biologist" ( or application-level options)\
**

The steps ( or terms ) that can be used within a pipeline ( "p" command
) are defined in a file called "functionInterpretations.csv". That file
also defines the function signatures and defaults for the various
arguments. This file is used to define what options are apt to be shown
to a biologist or other user. That is, within the Java UI and other
components, not all arguments are necessarily exposed to the user. Or,
there may be defauls for a call that are specific to a given application
area, and appropriate for use in analysis pipelines in that area. So,
they are not necessarily the same as the underlying defaults for the "R"
command

One other aspect of the file "functionInterpretations.csv", and the
system as as whole, is that it also separates arguments in terms of
important dataset, matrix and vector values, and handles those
separately from strings of scalar arguments. In this way, although it is
certainly an option to do so, the Java UI does not necessarily require a
user to specify each and every scalar value.

**Use by this system of a "services layer" over R functions.**

One practical application is to be able to easily handle 3-10 gigabytes
of LC/MS files in that domain, as well as to gridify existing R
routines. For this and other reason we implemented the various commands
using a "services layer" over existing ( and some new ) R functions.
This services layer makes use of notions of "pipelines" and "types".

The basic aspects of the types/pipeline layer in RProteomics are as
follows:

1.  For use in automated programming and other purposes, we define
    specifications for functions, including information on parameter
    types. Using that type information in different forms, we generate
    simple pipeline/grid wrappers and other code.

2.  We define a high-level workflow/pipeline language in which the
    "interpretations" of terms are different R functions used with
    specified list of defaults

3.  We allow different basic data types (e.g., matrices, vectors ) to be
    stored in different ways, such as in different kinds of files, so
    that a pipeline can integrate different functions in as natural a
    way as possible.

4.  We make use of the above through the use of a simple pipeline (
    workflow) engine within R, which can run through the datasets
    according to a couple of ( relatively simple ) approaches.

The section below considers each of these parts in turn.

**1. Defining function signatures.**

Information on function signatures is stored in several different ways.
For example, we make some use of an "AnalysisML", that is, an analysis
markup language; this is discussed at the end of this note. We also do
some writing and also reading of "Rd" files. However, over time, the
starting point has involved simply annotating R functions to define
their signature information and other details for use in the overall
pipeline. For example, the annotation of an external peak finder in
caMassClass is this:

#msc.peaks.find=function(X, SNR=2, span1=81, span2=11, zerothresh=0.9)

#Signature
X:matrix,SNR:real,span1:integer,span2:integer,zerothresh:real=\>matrix

#Keywords:Peaks

#Title: caMassClass peak finder

#Description: Finds peaks using running quantiles and other approaches.

#Example: variable(\'peaks\',storageType=\'csvFile\');

\# pipelineRun \<- p(\'peaks = msc.peaks.find(I,SNR=2.5)\')

The second line of the above is for the signature, which includes
information on the various data types. The last line in the above gives
an example of possible usage within the pipeline. (In the example, the
\"p\" is a synonym for \"runPipeline\", in this case involving running a
pipeline of one step). In this particular case, instead of the data
being created within an .RData file or some other binary form, data is
stored in the form of a single CSV file with multiple columns.

**2. Generating various wrappers, such as "pipeline wrappers".**

Given the basic data types and signatures of functions (that is, the
initial algebra or abstract syntax overall of the system), of course, it
is then possible to generate all sorts of wrappers. It certainly also
makes sense to take a ( complementary ) approach of generating R code to
handle such things as WSDL or CORBA wrappers as well, all in R. However,
we did not happen to do that; instead, we generated simple pipeline
wrappers. With the help of those wrappers, other aspects of the system
are then done not in R, but using Java and Java-based web services.

Within a pipeline, the "options" are collected together and passed
together as a list of R options. That list of options as well as other
information at a given pipeline step is also available within a
general-purpose \"callingParameters\" object. As used within the
pipeline, and in addition to the "callingParameters", the other
parameters are the main pieces of data, such as matrices.

As implied, the pipeline calls individual R functions with the help of
an automatically-generated "wrapper" function. Over time, I have changed
the wrappers in different ways. An example what one currently looks like
follows ( one line has been changed by hand.)

\###############################################################################

msc.peaks.findWrapper \<- function(callingParameters,X)

{

#Do the basic conversions from the list of options with defaults
populated

options \<- callingParameters\$options

SNRString \<- options\$SNR

if ( is.null(SNRString) )

{

SNRString \<- 2

}

SNR \<- as.double(SNRString)

span1String \<- options\$span1

if ( is.null(span1String) )

{

span1String \<- 81

}

span1 \<- as.integer(span1String)

\...\[deleted here\]

zerothreshString \<- options\$zerothresh

if ( is.null(zerothreshString) )

{

zerothreshString \<- 0.9

}

zerothresh \<- as.double(zerothreshString)

#rph: next line tweaked by hand

result \<- try ( msc.peaks.find (X=X, SNR= SNR, span = c(span1, span2),
zerothresh= zerothresh))

if (inherits(result, \'try-error\'))

{

logException(\"msc.peaks.find\", result)

return ( NULL )

}

return ( result)

}

\###############################################################################

We have found it helpful to generate the types-related information in
other forms, such as in \"good old spreadsheet\" form. For example, the
following is a subset of one of my spreadsheets that contains sets of
signatures as well as other information.

![](media/image1.png)

The above includes columns for:

-   key words

-   the name of the function as used within a pipeline ( \"model\" )
    syntax

-   the name of the actual R function used

-   the \"main\" signature as it involves vectors, matrices and other
    crucial pieces

-   (somewhat cropped ) a separate list of "options" and their types

-   (not shown ) default parameters assigned to the R function for this
    particular model name as used in the pipeline.

As noted above, the system attempts to the "main data" ( "main
parameters" ) from "user options".

Using the above spreadsheet, a file that describes and gives examples of
various basic types, and

files of "R" examples, we also generate various other files. In this
context, it is sometimes useful to have "Rd" files that contain content
that pertains to parameter data types. So, for a simple set of
nonparametric distance functions we sometimes generate sections of an Rd
file as follows:

...

\\usage{

results = distance(I, I2,type:STRING,exponent:REAL) #as a term within a
pipeline, or,

pipelineRun \<- runPipeline(\'results = distance(I,
I2,type:STRING,exponent:REAL)\')}

\\arguments{

\\item{I}{ Matrix of intensities, e.g., I or I\[10:20,510:1010\] }

\\item{I2}{ Matrix of intensities, e.g., I2 or I2\[10:20,510:1010\]

...

}

}

...

\\examples{

testName \<- \'distance: determine nonparametric quadratic distance
measure in a vector by vector way\'; dataset(\'patz\')

pipelineRun \<- p( \"distanceMinkowski =
distance(I\[,5:6\],I2\[5:6\],type=\'Minkowski\')\")

testName \<- \'distance_1: determine nonparametric distance measure in a
vector by vector way\'; dataset(\'patz\')

pipelineRun \<- p( \"distanceChiSquare =
distance(I\[,5:6\],I2\[5:6\],type=\'chiSquare\')\")

testName \<- \'distance_2: determine nonparametric distance measure in a
vector by vector way\'; dataset(\'patz\')

pipelineRun \<- p( \"distanceKL =
distance(I\[,5:6\],I2\[5:6\],type=\'KL\')\")

testName \<- \'distance_3: determine nonparametric distance measure in a
vector by vector way\'; dataset(\'patz\')

pipelineRun \<- p( \"distanceGeometric =
distance(I\[,5:6\],I2\[5:6\],type=\'Bhattacharyya\')\")

\[ At this point we have not yet integrated with Sweave.\]

The language for doing the automated wrapper generation and these other
components happens to use PERL. In particular, if types were available
within R in a native way, then pure R/SWIG itself might be better than
PERL to use.

**3. Allowing different data types to be stored in different ways ( that
is, have different \"interpretations\" )**

Much of the data handling in our system ( such as code to handle the XML
) is done using C/C++ or Java code. However, it has always made sense to
allow code within our R code to at least be able to handle some number
of data formats by itself. With that in mind, for example, we defined an
approach called StatML for wrapping datasets using "generic XML". ( At
this point, as noted, it still can only wrap matrices, not true data
frames. ) For general pipeline use, however, basic binary files have
been convenient and fine; two variations include having a single file
containing the matrix data in binary form; another variation involves
matrices with columns or rows split out as individual files. This does
not work for arbitrary R classes, of course. The idea was to master the
basic types of matrices and data frames, and only then move to use of
arbitrary R classes and more complex data formats such as netCDF.

For doing this, we used a very simple data structure. An example is:

fileInfo \<- list ( dataFormat = \"double\",

inputFilePrefix = \"cancerous\",

inputFileSuffix = \".dat\",

outputFileSuffix = \".dat\" )

semanticModel \<- list( semanticModelType = \"files\",

inputFunction = \"getSetOfVectorsInBinaryFiles\",

outputFunction = \"putSetOfVectorsInBinaryFiles\",

fileInfo = fileInfo )

As is done in some other systems, all of the input functions have the
same signature. The same holds for all of the output functions[^1]. In
this way, the system can store the data and/or communicate with its Java
modules in any desired way.

In this context, a sequence of commands at the application level such
as:

variable(\'J\',type=\"matrix\", dataset=\"patz\",

storageType=\'binaryFile\',rootDirectory=\"l:/\");

p(\'J = runStep(I)\' )

then puts the results of J as a single binary file within the
rootDirectory \"l:/patz/J\".

**4. The basic pipeline model**

As I believe has been done elsewhere, it is possible to just define
matrix indexing \[i,j\] operators to make use of data in databases, on
disk, or provided by grid services. The final decision to go with a
pipeline was in part to be able to define pipelines to do hierarchical
decompositions of spectra, and do that in an iterative manner. In other
words, a pipeline involves just a type of formula associated with the
solution of some iterative method. In this case, the locus of control
needs to be at the pipeline ( model ) layer. As implied, there are both
advantages and disadvantages of this approach.

The basic pipeline needed to be able to have slightly different steps
for handling functions of very different kinds of signatures. It also
needed to move data through the system ( that is, process matrices and
datasets ) in different ways, such as:

-   By handling matrices as individual sets of columns or rows ( one by
    one ) in order to reduce memory needs.

-   By handling matrices all at once, with all data in memory at one
    time.

-   By handling matrices as sets of vectors, but where the final result
    is a single vector.

Different processing methods involve the use of slightly different
pipeline procedures. For example:

if ( callingParameters\$ProcessType == \'vectorsInMatrixOut\' )

{

pipelineRun \<- applyFunctionToSetsOfVectors(callingParameters,

> rowAttributes,semanticModel,debug)

}

The above call acts on matrices as sets of vectors, that is, runs data
through the pipeline by moving just one column of data through each
procedure at a time. Other approaches apply to other ways of moving the
data through the system.

**Conclusion**

The above summarizes some basic aspects of the RProteomics package,
together with ways by which we automatically add new pipeline/grid
wrappers for use within the RProteomics pipeline.

**Footnote on an \"AnalysisML\" we developed**

We also developed an \"AnalysisML\". That included terms for
pre-conditions, post-conditions and other invariants, as well as
exceptions. The names of these terms ( e.g., \"contract\" ) were the
typical computer programming versions of those terms, rather than using,
say, the vocabulary of algebraic types.

However, editing the actual XML text was actually not as simple as
editing simple \"Rd\" files in R, or just defining the signatures as
annotations in short \"R\" files that just defined signatures. Another
minor issue is that there were some additional fields in our
\"AnalysisML\" we probably still might want to add. Finally, the overlap
with Rd files was in the nature of an overlap, rather than a strict
superset or subset. Having said that, we still hope to make use of our (
or someone else\'s ) AnalysisML to a greater extent in the future.

[^1]: This also applies to a set of routines for getting data in and out
    of database tables. However, since our data comes indirectly through
    XML-based databases, and through our Java layer, we have never
    actually integrated those in practice.