BACKGROUND
peptide microarrays have emerged as a promising technique for the simultaneous high throughput analysis of peptide characteristics.
synthesized peptides are spotted in a grid-layout on glass slides which allow the screening of thousands of peptides within a single experiment with requiring only a small quantity of sample.
applications range from studying the humoral response to hiv  <cit>  or food allergens  <cit>  to the detection of cancer biomarkers  <cit>  and antibody signatures  <cit>  to the characterization of protein-protein interactions  <cit>  and of kinase substrates  <cit> .
while peptide microarrays offer enormous potential for a wide range of applications and while significant improvements have been made with regard to the reliable spotting of small amounts of peptides at closely neighboring, well-defined spatial positions  <cit> , a major bottle neck remains in the automated analysis of the acquired data.
numerous tools have been developed for the analysis of dna microarray data, see e.g.
<cit>  for reviews, but these cannot easily be transferred and require major adjustments.
for instance, rather than quantifying the impact of differential expression, peptide microarrays experiments commonly only result in a single wavelength measurement per peptide.
usually, only a small proportion of the spotted peptides is expected to show a signal and requires reliable identification.
further specific challenges of the analysis of peptide microarray data include the diverse sources of noise, ranging from peptide synthesis artifacts to unspecific binding effects to peptides.
existing tools specifically developed for the analysis of peptide microarray experiments can be categorized into three groups, graphical analysis tools, differential analysis tools, and general analysis tools.
as graphical tools, a combined approach of clustering and principal component analysis to visualize similarly behaving groups of peptides has been introduced  <cit>  as well as an integrated webserver for the storage of peptide microarray experiment data with several graphical analysis steps  <cit> .
for differential analysis, several approaches have adapted differential expression detection schemes to peptide microarrays, including a support vector machine driven webtool for distinguishing peptide binding intensities of two experimental groups  <cit>  and adapted statistical tests for differentiating measured intensities for two populations  <cit> .
the more general question of identifying signal carrying peptide spots and accounting for peptide microarray specific sources of noise has been addressed using a robust version of a z-score for the difference of the intensity of a specific peptide spot to empty spots to identify signal carrying peptide spots  <cit>  and a linear model fit on all peptide spot measurements to account for several systematic effects  <cit> .
this approach was extended by a signal calling step based on a t-test and a cutoff based removal of secondary binding spots  <cit> .
here, we introduce a novel method for the general analysis of peptide microarray data, that significantly extends the existing approaches.
using several classes of control peptides, we apply a linear model for normalization and the removal of systematic array effects.
further, we use a mixture-model to identify secondary antibody binding peptide spots and apply a probabilistic approach for signal calling which does not rely on arbitrary thresholds, but provides a slide specific estimate.
additionally, we provide a machine-learning driven quality control procedure to computationally exclude intensity measurements of low reliability.
after describing our methods in detail, we apply it to data from a cancer-biomarker detection study and demonstrate improvements relative to the existing general analysis tools, particularly with regards to sensitivity.
method
peptide microarrays
the layout of the microarray slides used in this study is based on a three level hierarchy : each array  consists of three subarrays  which are identical in terms of individual peptide placement but are printed consecutively; thus, each individual peptide is spotted as a triplicate.
each subarray has  <dig> blocks  which are arranged in a four by four pattern; each of these blocks is printed by a single print tip.
each block has  <dig> rows and  <dig> columns, resulting in  <dig> peptide spots per array.
printed slides are scanned at  <dig> nm wavelength for quality control, with peptide spots being visible due to scatterlight effects.
in brief, the incubation experiment and image acquisition consist of the following steps:  blocking of the array surface to reduce unspecific binding,  incubation of the slides with sample material,  incubation with a fluorescently labeled secondary antibody and  scanning at  <dig> nm wavelength and spot detection.
we refer to a peptide as being reactive if we observe fluorescence of the labeled secondary antibody indicating that primary antibody of the sample bound to the peptide.
for a detailed protocol of microarray incubation experiments including the motivation of acquisition conditions see  <cit> .
it is noteworthy that our computational methods described in the further sections accept generic data as input and are not restricted to the array layout and procedures mentioned above.
algorithmic overview
the analysis pipeline is implemented in r  <cit>  and available from http://www.tron-mz.de/compmed.
read out and transformation of raw data
to transform an image file from the scan of a peptide microarray to individual intensity values for each single peptide, we apply genepix pro, version  <dig> , using the mean foreground intensity  from the genepix result  file as a measure for intensity.
it has been previously described  <cit>  that peptide microarray measurements show strong heteroscedasticity with larger variances for high intensity measurements.
by applying a binary logarithm transformation to each peptide spot mean intensity, we stabilize the variance allowing the application of a linear model procedure.
control peptide based linear model fit
peptide microarray intensity measurements can be perturbed by a variety of factors, some of which have been proposed in the model of  <cit> .
these factors include array-to-array effect  resulting from the printing process or storage conditions.
further, depending on the spatial position of a peptide spot, it might be incubated with a different amount of sample.
this relates to the overall row  and column  position of a spot as well as to the subarray  to which a spot belongs.
further, variation might be introduced by the printing needle  used.
the effect of a peptide itself is then denoted by pn.
thus, an intensity measurement yijklmn can be modeled as
with the variance stabilizing logarithm transformation applied in the previous step, the residuals  can be assumed to follow a normal distribution with constant variance.
in most applications, the number of repeats per peptides are very limited and the placement of peptides is not in random order, but rather defined by practical considerations of the printing process.
thus, confounding of effects can arise, for instance the placement of predominantly strongly reacting peptides at positions spotted with one needle might result in the wrongful attribution of a strongly positive needle effect.
to avoid these confounding effects, we limit the fitting of the linear model of equation  <dig> to a set of control peptides, which are printed with high number of repeats for each needle and subarray to allow reliable estimation.
while this limitation renders the estimation of a peptide effect pn infeasible, all remaining coefficients can be reliably estimated and can be used to correct all peptide measurements.
this allows us to estimate the parameters ai, rj, ck, sl, and nm for all values of i, j, k, l, and m. we then correct all measurements, for which we could not estimate the peptide specific effect pn, by computing .
unreliable spot removal
when historic data is available for peptide arrays with identical slide layout, this data can be used to remove unreliable spots.
since this is not possible in all setups, we regard this as an optional step for improving and robustifying peptide microarray data analysis.
we apply a random forest procedure  <cit>  to identify potentially unreliable peptide spots  <cit> .
unreliable spots may arise from a variety of not explicitly modeled experimental effects, e.g.
the different decay of spotted peptides, dust effecting peptide spots or errors in the alignment of spots by software such as genepix.
these variation can drastically reduce the reproducibility of peptide array experiments and may result in misleading conclusions.
using historic data from previous experiments with slides of an identical layout, we train a random forest classifier to detect spots which are not reproducible between subarrays.
a random forest classifier is a machine learning procedure, which learns label categories from existing data and then applies them to new, previously unclassified data sets.
it shows good predictive power while being robust to parameter setting  <cit>  and favorable theoretical properties .
for the training of the classifier, we use the fact that each peptide is repeated on various subarrays .
we label all subarray spots which are not within a 95% linear regression confidence band for a single linear regression on either of the remaining two subarray repeats as unreliable and all other spots as reliable.
we use all columns of the result file besides the sequence as well as all columns of the quality control result file of a corresponding scatterlight scan for each spot as explanatory variable .
we use random sampling to balance the set of unreliable and reliable spots in their size and then use the prediction of the random forest to label all spots of a new data set.
for each peptide, we directly discard all those spots indicated as unreliable and only average over the remaining spots, weighted by their probability of being reliable.
if all spots of one peptide are marked as unreliable, we exclude this peptide from any further analysis.
mixture model for secondary antibody binding detection
in peptide microarray experiments, secondary antibody binding peptides pose a major complication.
when peptides show a direct reaction with the dye-coupled secondary antibody, two effects are confounded.
it cannot be distinguished to which extent a measured intensity is explained by this undesired and unspecific reaction rather than the specific reaction of interest.
thus, to avoid ambiguous results which hinder interpretation and may result in misleading conclusions, peptides showing significant unspecific reaction with the secondary antibody need to be identified and excluded from further analysis.
reactivities on empty slides without any sample have been identified as important indicators  <cit> .
in experimental settings not allowing empty slides, this step can be excluded from the overall data analysis pipeline with an increased risk of false positives.
however, this approach is tailored to identifying the unspecific binding of the secondary antibody to the peptide, not to the slide itself, which should affect all spots in a similar mode and thereby increase the noise level.
this effect is addressed by the subsequent signal calling step.
rather than relying on a fixed threshold, we propose a mixture model strategy to derive a probabilistic criterion for the identification of secondary antibody binding peptides from empty slides.
for this purpose, we assume that the normalized intensities of all peptides on the empty slides stem from two populations: one population for peptides, which do not show any reaction to the secondary antibody besides noise and a second population for peptides which show a reaction to the secondary antibody.
we model each population by a normal distribution and let the mean and variance of the first population be denoted by μ <dig> and  and for the second population by μ <dig> and .
a schematic display of a histogram of the measured intensities and the two underlying distributions in red  and green  is given in figure 1d.
the orange line indicates the resulting cutoff.
we apply an expectation maximization algorithm  <cit>  to iteratively adapt the estimates for μ <dig>  μ <dig>  , and , and for each measured normalized intensity, we update the likelihood that it belongs to either one of the two populations  <cit> .
for the initialization, we rely on the empirical mean and variance estimate for two groups of control peptides.
one group for which we do not expect any reaction with the secondary antibody and a second group of peptides which are known to show various degrees of reactivity.
we then use the resulting distribution of the peptides which do not show any reaction to the secondary antibody, to derive a criterion for secondary antibody binding.
since we expect most peptides not to show any reaction, we expect this distribution to be more reliably estimated and mostly driven by random variations.
we regard the upper 5% quantile  of this distribution to be the cutoff.
thus, for any measured normalized intensity on the empty slides , we compute
any resulting value of  below  <dig>  indicates only low likelihood for this intensity being this large by chance  <cit>  and thus the corresponding peptide is disregarded for the further analysis.
with this cutoff generation we rely on the assumption that at least the population of peptides not showing a secondary antibody binding follow a normal distribution.
however our historic slides at least approximately support this assumption .
instead of this cutoff generation, it is also possible to directly compute a false discovery rate from the likelihood of the mixture model components , however this relies more heavily on the correctness of the normality assumption for the group of the reactive peptide than the proposed approach, which depends more on the not reactive peptides, in which we are more confident.
mixture model for signal call
to distinguish those peptide spots which carry a signal from those spots which only display noise, we again apply a mixture model.
we assume that both, the signal-carrying peptide spots as well as the noise spots can be modeled by normal distributions with mean μsignal respectively μnoise and standard deviations  respectively .
a schematic display of a histogram of the measured intensities and the two underlying distributions in red  and green  is given in figure 1e; the orange line indicates the resulting cutoff.
we initialize these parameters again based on suitable control classes.
for one class of peptides expected to show a reaction, and for a second class of peptides not expected to show a reaction, we separately compute the empirical mean and variance and use these as starting values.
we then update these values as well as the likelihood for each data point using the expectation maximization algorithm as described for the secondary antibody binding detection.
we then use the distribution of the noise spots to identify the signal carrying spots which depart significantly from the noise distribution by identifying the upper 5% quantile  of the noise distribution as a cutoff.
thus, for any measured normalized intensity , we compute  analog to equation  <dig>  any resulting value of  below  <dig>  indicates only low likelihood for this intensity being this large by chance.
experiments
experiments are conducted on a total of  <dig>  peptide spots which were printed on six pepstar peptide microarrays  with an identical design and printed in two batches .
two separate batches were chosen to show the maximum possible inter-experiment variation.
on each slide 10% of the total  <dig>  spots are used as dedicated control spots.
5% of all peptides are meant as positive controls, expected to show strong specific reactions,  <dig> % are secondary antibody controls, which are expected to show strong reactions with the secondary antibody.
1% of all peptides are negative controls not expected to show any reaction and the remaining  <dig> % are process controls.
among the remaining  <dig>  peptide spots per array,  <dig> spots contain peptides of human plac- <dig>  and  <dig> spots contain peptides of nyeso, vegfa, and cd <dig> .
all proteins are represented by peptides with a length of  <dig> amino acids and an overlap of  <dig> amino acids.
the slides were incubated in a hs  <dig> pro hybridization station .
one slide per print batch was incubated with diluent only , one slide per batch with a pooled plasma of six patients without indication of cancer and low concentration  of a plac1-specific spike-in antibody  and the remaining slide per batch with the same sample material, but high concentration  of the plac1-specific spike-in antibody.
the experimental setup is summarized in table  <dig>  a cy5-conjugated affinipure mouse anti-human igg   was used as a secondary antibody.
all slides were scanned with a genepix  <dig> microarray scanner  and converted to a result file using genepix pro, version  <dig> .
all result files are available from the authors' website.
summary of the experimental setup for the six microarrays used in this study.
the analysis of results is focused on five aspects: we evaluate  the efficiency of the linear model in reducing the unexplained variance,  the unreliable spot finding procedure by comparing the reproducibility of the two print batches before and after quality control,  the removal of secondary antibody binding peptides by focusing on the nyeso, vegfa, and cd <dig> peptides, which are not expected to show any reaction with the given normal sera, with the exception of secondary binding effects, and  the efficiency of the signal call by analyzing the specificity and sensitivity of our approach based on the nyeso, vegfa, and cd <dig> peptides which are expected to be non-reactive and the plac1-peptides, which we expect to show a reaction.
we regard nyeso, vegfa, and cd <dig> peptides showing a reaction as false positives  and those not showing a reaction as true negatives .
similarly, we regard plac- <dig> peptides not showing a reaction as false negatives  and those showing a reaction as true positives .
sensitivity, specificity, and accuracy are then given by
in a final step , we compare our approach to the approaches introduced by  <cit>  and  <cit>  based on the accuracy.
to judge the significance of the results of experiments -, we computed 95% bootstrapped confidence intervals for all sensitivity, specificity and accuracy values based on  <dig> times resampling from the peptide intensities of each class.
in order not to confound the results of the unreliable spot finding step  with the subsequent steps and to allow the fair comparison with the existing methods which do not identify unreliable spots, experiments - were run on all peptide spots, even though some were marked as unreliable.
RESULTS
we analyze the effect of the major steps of our approach including the linear model, the quality control and the secondary antibody binding detection and signal call as well as comparing the overall procedure on the antibody reactivity data.
with regard to the linear model fit, a summary of results for print batch  <dig> and the high spike-in antibody concentration is given in table  <dig>  it is evident that the peptide effect itself is the major influence of the measured intensity.
still, all remaining explanatory variables, which model systematic effects, are highly significant as shown by the p-values for the corresponding f-test.
row and column effects, so the spatial position of a peptide, show the overall strongest influence of the systematic effects.
the needle and subarray effect are smaller, but still highly significant.
however, it should be noted that a confounding of effects cannot be excluded, e.g.
when considering a single needle only a subset of row and column positions are addressable on a subarray.
including all interaction effects requires an unfeasibly high number of degrees of freedom, while choosing subsets does not significantly improve results .
overall, the systematic effects contribute to reducing the unexplained variation in the measured peptide intensity data by 65%.
linear model fit summary.
the f-values and the corresponding probabilities ) clearly indicate that all explanatory variables used are highly significant and that each contributes to reducing the variance present in the arrays.
the peptide sequence itself shows by far the strongest effect while requiring only twelve degrees of freedom , it shows a sum of squares for its effect  of  <dig> , much stronger than the residuals sum of squares of  <dig>  while not as strong as the peptide effect, the remaining explanatory effects reduce the residuals by 65%.
column and row effects are strongest, but subarray and needle effects require fewer degrees of freedom, resulting in strong mean sum of squares  values.
after training on a total of  <dig> historic slides with identical layout, the unreliable spot finding procedure marks  <dig> % of all spots in print batch  <dig> and  <dig> % of all spots in print batch  <dig> as unreliable for the high concentration experiment.
since peptides are excluded only in cases when all three replicates for one peptide are marked as unreliable, this results in the removal of  <dig> % of all peptides in print batch  <dig> and  <dig> % of all peptides in print batch  <dig>  visual inspection  shows that these peptides primarily show large variation between the two data sets.
by removing these peptide spots, we see an improvement of 3% in the coefficient of variation from r <dig> =  <dig>  to r <dig> =  <dig> .
still, not all outliers are detected by the approach.
analysis of the missed outliers however indicates that these undetected outliers consistently affect replicates of only nine different peptide sequences, in contrast to the  <dig> different peptide sequences marked as unreliable by the unreliable spot finding procedure.
we evaluated our secondary antibody binding removal approach by computing the number of peptides which were excluded from the further analysis because of secondary antibody removal.
out of  <dig> non control peptides 354/ <dig>  were identified with reaction on the empty slides and excluded from the further analysis.
this affected none of the plac- <dig> peptides which were expected to show a specific reaction and 0/ <dig>  of the nyeso, vegfa and cd <dig> peptides not expected to show a reaction.
consequently, when regarding the sensitivity and specificity which were computed based on these peptides, we observe that with the secondary antibody binding removal the specificity is slightly increased while the sensitivity remains unchanged when compared to running our approach without this step.
the sensitivity and specificity for our approach are shown in figure  <dig> in comparison to the approaches of  <cit>  and  <cit> .
for two print batches, two slides were incubated with a low and a high concentration of spike-in antibody, respectively.
overall, we observe for all approaches that sensitivity is more affected by the reduction of the antibody concentration than specificity, which is rather stable.
in comparison, our approach shows best sensitivity across all settings.
it shows a sensitivity of  <dig>  to  <dig>  compared to  <dig>  to  <dig>  for the other two approaches for the high antibody concentration data.
this difference is more expressed for the low antibody concentration data with a sensitivity of  <dig>  to  <dig>  for our approach compared to  <dig>  to  <dig>  for the previous two approaches.
specificity is slightly worse for our approaches for print batch  <dig> with values of  <dig>  to  <dig>  for our approach and  <dig>  to  <dig>  for the existing approaches.
similarly for print batch  <dig>  our approach features specificity of  <dig>  to  <dig>  compared to  <dig>  to  <dig> .
as a consequence of the substantially improved sensitivity, our approach also displays best accuracy.
while this difference is small for the high antibody concentration with accuracy values of  <dig>  to  <dig>  for our approach compared to  <dig>  to  <dig>  for the previous approaches, it is much more substantial for low antibody concentration.
here, our approach shows an accuracy ranging from  <dig>  to  <dig> , compared to  <dig>  to  <dig>  for the two other approaches.
rapmad is based on the assumption of  normality, especially for the linear model fit and the mixture model.
for all arrays incubated with plasma, we analyzed these assumptions and find support for approximate normality .
discussion
removal of systematic bias
by applying a linear model which explains 65% of the noise variation in the data, rapmad significantly improves the reliability of peptide microarray experiments.
in addition to previous findings  <cit> , we establish that also the spatial position on the peptide microarray is of high impact for removing systematic effects.
this is important for distinguishing low intensity signals from noise.
since the effects are of a systematic nature , it also reduces the risk that systematic effects are mistakenly identified as peptide signal effects.
despite the strong reduction in unexplained variation, not all occurring systematic effects are necessarily included in the model.
for instance, we decided against including all interaction effects.
while it can be of high interest to model the interaction of a specific peptide with a specific needle, there is usually not a sufficient number of replicates to allow a reliable estimation in a standard experiment.
however, for large-scale experiments, it could be beneficial to design specific slides with sufficient replicates to quantify these effects.
the effects of the secondary antibody are accounted for within subsequent analysis steps and are thus not modeled explicitly.
improved data quality after removal of unreliable measurements and secondary binding peptides
the production of peptide arrays as well as the binding of the secondary antibody involve sophisticated chemical processes.
in spite of numerous experimental improvements  <cit>  and constant monitoring, artifacts resulting from errors in the spotting of peptides or the unspecific binding of the secondary antibody can strongly influence measurements and result in misleading interpretations.
we have proposed two computational statistics driven approaches to identify and exclude these artifacts and show that they can be successfully applied.
by marking approximately  <dig> % of the peptide spots as unreliable, we can improve the reproducibility and thus the data quality.
it should be noted however that still not all unreliable and irreproducible peptides are detected.
the small number of peptide sequences being connected to the undetected outlying observations show that this entails a sequence-specific phenomenon rather than a print process artefact which would affect all sequences more evenly.
while characteristics of the print process are captured within the available learning data, for instance in the scatterlight scans, sequence-specific effects are difficult to detect.
in addition, the unreliable peptide spot finding procedure is limited by the available training data.
we rely on labels based on the assumption that a specific replicate which departs from the remaining two replicates should be labeled as unreliable.
we are aware that this definition is not necessarily correct since there are cases when the two remaining replicates are unreliable, but the specific replicate itself is correct.
this can result in incorrect labels, which strongly reduce the quality of the random forest prediction.
while manually validation could improve the label-quality, this reduces the applicability of the procedure to other data than our own.
thus, we accept a possibly smaller impact of the quality control procedure to allow the fully automated adaption to other datasets.
historic data for training is not always available and thus, this is an optional step for our pipeline.
by identifying  <dig> % of the spots as secondary antibody binding, we can improve the specificity of the procedure without sacrifying its sensitivity.
while the impact on the data shown here is comparatively small, this step can be of key importance to avoid wrongful interpretation.
when dealing with peptides prone to printing or synthesis problems or the binding of secondary antibody to peptides, any data analysis tool is challenged by the fact that the relevant information was lost before data acquisition, rendering the complete recovery of the information impossible.
while not recovering the correct information, the magnitude of excluding wrongful results is significant.
if wrongful results were not excluded, research could result in wrongful conclusions.
with the removal of wrongful results no conclusions are possible for these doubtful peptide spots.
since in most cases, only subsets of peptide measurements are affected and designs usually contain several replicates, repeats of experiments are often not required to come to reliable conclusions.
strong sensitivity for low intensities
a particular strength of rapmad is the fact that it is especially suitable for low intensity measurements such as those occurring for low antibody concentrations.
while the problem of separating signal from noise is substantially more difficult when the signal intensity is reduced, we only see a comparatively minor reduction in accuracy when reducing spike-in antibody concentration.
sensitivity is affected more significantly than specificity, but still at or above  <dig> .
this robustness of our approach can be explained by the adaptivity of the signal calling procedure, which is based on correctly estimating the noise distribution and generating a noise-threshold based on this distribution.
thereby, it is less affected by changes in the signal intensity.
as a consequence, the noise-threshold is adaptive and even low intensity signals can be detected to be above the threshold, resulting in a strong sensitivity.
favorable comparison to existing approaches
the strength of rapmad to maintain a high sensitivity even for low intensity settings is also what sets it apart from existing algorithms.
thus, the sensitivity and accuracy of our approach in low antibody concentration experiment is competitive to the ones of the two other approaches for high antibody concentrations.
this is a result of the adaptive procedures which are better suited to distinguish signal from noise in low intensity situation than preset thresholds and demonstrates that our computational approach is capable of compensating for experimental limitations.
for high intensities as seen with high antibody concentration, it is generally easier to achieve high sensitivity due to a better separation of signal and noise in the raw data.
that is also why we see less of a benefit of rapmad in comparison to existing algorithms.
increasing the effectiveness of screening procedures
peptide microarrays have been primarily used within screening settings with follow-up experiments focusing on a small number of promising candidates to ensure specificity  <cit> .
unidentified reactive peptides from a peptide microarray experiment cannot be distinguished from the large number of rightfully excluded unreactive peptides and are thus not further regarded.
by substantially increasing the sensitivity for low intensity reactions, our approach avoids the wrongful exclusion of peptides and thereby increases the effectiveness of screening procedures.
while standard settings were used in all experiments, rapmad additionally allows the user to adapt the quantiles used in the mixture models to trade additional sensitivity for specificity.
CONCLUSIONS
within this contribution, we introduce rapmad, a tool for the robust and rapid analysis of peptide microarray data.
rapmad is an automated, multi-step approach that combines several computational statistics procedures to augment the data quality of peptide microarray and to allow a reliable analysis.
its steps include the preprocessing of the data by removing systematic effects, the exclusion of unreliable measurements and secondary antibody binding peptides, and a probabilistic signal call for reactive peptides.
in comparison with existing algorithms, it shows competitive performance for high antibody concentration.
for low antibody concentration, it shows a significant increase in accuracy over existing approaches.
this is mainly due to its substantially improved sensitivity at competitive specificity.
with its increased sensitivity and its automated data analysis, rapmad can thereby contribute to establish and broaden the usage of peptide microarrays as a standard tool for a wide-range of peptidomics applications.
authors' contributions
byr, ml, ur, jcc, and us designed the methods, byr and ml wrote code and analyzed data.
all authors contributed to the design of the experiments.
yk and ar set up and ran the experiments.
byr wrote the manuscript with contributions from ml and jcc.
all authors read and approved the manuscript.
supplementary material
additional file 1
supplementary figures.
pdf document containing supplementary figures.
click here for file
acknowledgements
the authors would like to thank christina fleig  for data acquisition, antonia masch  for critical reading of the manuscript and ganymed ag, mainz for providing antibodies.

 peptide microarrays have emerged as a promising technique for the simultaneous high throughput analysis of peptide characteristics. synthesized peptides are spotted in a grid-layout on glass slides which allow the screening of thousands of peptides within a single experiment with requiring only a small quantity of sample. applications range from studying the humoral response to hiv or food allergens to the detection of cancer biomarkers and antibody signatures to the characterization of protein-protein interactions and of kinase substrates. while peptide microarrays offer enormous potential for a wide range of applications and while significant improvements have been made with regard to the reliable spotting of small amounts of peptides at closely neighboring, well-defined spatial positions, a major bottle neck remains in the automated analysis of the acquired data. numerous tools have been developed for the analysis of dna microarray data, see e.g. for reviews, but these cannot easily be transferred and require major adjustments. for instance, rather than quantifying the impact of differential expression, peptide microarrays experiments commonly only result in a single wavelength measurement per peptide. usually, only a small proportion of the spotted peptides is expected to show a signal and requires reliable identification. further specific challenges of the analysis of peptide microarray data include the diverse sources of noise, ranging from peptide synthesis artifacts to unspecific binding effects to peptides. existing tools specifically developed for the analysis of peptide microarray experiments can be categorized into three groups, graphical analysis tools, differential analysis tools, and general analysis tools. as graphical tools, a combined approach of clustering and principal component analysis to visualize similarly behaving groups of peptides has been introduced as well as an integrated webserver for the storage of peptide microarray experiment data with several graphical analysis steps. for differential analysis, several approaches have adapted differential expression detection schemes to peptide microarrays, including a support vector machine driven webtool for distinguishing peptide binding intensities of two experimental groups and adapted statistical tests for differentiating measured intensities for two populations. the more general question of identifying signal carrying peptide spots and accounting for peptide microarray specific sources of noise has been addressed using a robust version of a z-score for the difference of the intensity of a specific peptide spot to empty spots to identify signal carrying peptide spots and a linear model fit on all peptide spot measurements to account for several systematic effects. this approach was extended by a signal calling step based on a t-test and a cutoff based removal of secondary binding spots. here, we introduce a novel method for the general analysis of peptide microarray data, that significantly extends the existing approaches. using several classes of control peptides, we apply a linear model for normalization and the removal of systematic array effects. further, we use a mixture-model to identify secondary antibody binding peptide spots and apply a probabilistic approach for signal calling which does not rely on arbitrary thresholds, but provides a slide specific estimate. additionally, we provide a machine-learning driven quality control procedure to computationally exclude intensity measurements of low reliability. after describing our methods in detail, we apply it to data from a cancer-biomarker detection study and demonstrate improvements relative to the existing general analysis tools, particularly with regards to sensitivity. method peptide microarrays the layout of the microarray slides used in this study is based on a three level hierarchy: each array consists of three subarrays which are identical in terms of individual peptide placement but are printed consecutively; thus, each individual peptide is spotted as a triplicate. each subarray has blocks which are arranged in a four by four pattern; each of these blocks is printed by a single print tip. each block has rows and columns, resulting in peptide spots per array. printed slides are scanned at nm wavelength for quality control, with peptide spots being visible due to scatterlight effects. in brief, the incubation experiment and image acquisition consist of the following steps: blocking of the array surface to reduce unspecific binding, incubation of the slides with sample material, incubation with a fluorescently labeled secondary antibody and scanning at nm wavelength and spot detection. we refer to a peptide as being reactive if we observe fluorescence of the labeled secondary antibody indicating that primary antibody of the sample bound to the peptide. for a detailed protocol of microarray incubation experiments including the motivation of acquisition conditions see. it is noteworthy that our computational methods described in the further sections accept generic data as input and are not restricted to the array layout and procedures mentioned above. algorithmic overview the analysis pipeline is implemented in r and available from http://www.tron-mz.de/compmed. read out and transformation of raw data to transform an image file from the scan of a peptide microarray to individual intensity values for each single peptide, we apply genepix pro, version, using the mean foreground intensity from the genepix result file as a measure for intensity. it has been previously described that peptide microarray measurements show strong heteroscedasticity with larger variances for high intensity measurements. by applying a binary logarithm transformation to each peptide spot mean intensity, we stabilize the variance allowing the application of a linear model procedure. control peptide based linear model fit peptide microarray intensity measurements can be perturbed by a variety of factors, some of which have been proposed in the model of. these factors include array-to-array effect resulting from the printing process or storage conditions. further, depending on the spatial position of a peptide spot, it might be incubated with a different amount of sample. this relates to the overall row and column position of a spot as well as to the subarray to which a spot belongs. further, variation might be introduced by the printing needle used. the effect of a peptide itself is then denoted by pn. thus, an intensity measurement yijklmn can be modeled as with the variance stabilizing logarithm transformation applied in the previous step, the residuals can be assumed to follow a normal distribution with constant variance. in most applications, the number of repeats per peptides are very limited and the placement of peptides is not in random order, but rather defined by practical considerations of the printing process. thus, confounding of effects can arise, for instance the placement of predominantly strongly reacting peptides at positions spotted with one needle might result in the wrongful attribution of a strongly positive needle effect. to avoid these confounding effects, we limit the fitting of the linear model of equation to a set of control peptides, which are printed with high number of repeats for each needle and subarray to allow reliable estimation. while this limitation renders the estimation of a peptide effect pn infeasible, all remaining coefficients can be reliably estimated and can be used to correct all peptide measurements. this allows us to estimate the parameters ai, rj, ck, sl, and nm for all values of i, j, k, l, and m. we then correct all measurements, for which we could not estimate the peptide specific effect pn, by computing. unreliable spot removal when historic data is available for peptide arrays with identical slide layout, this data can be used to remove unreliable spots. since this is not possible in all setups, we regard this as an optional step for improving and robustifying peptide microarray data analysis. we apply a random forest procedure to identify potentially unreliable peptide spots. unreliable spots may arise from a variety of not explicitly modeled experimental effects, e.g. the different decay of spotted peptides, dust effecting peptide spots or errors in the alignment of spots by software such as genepix. these variation can drastically reduce the reproducibility of peptide array experiments and may result in misleading conclusions. using historic data from previous experiments with slides of an identical layout, we train a random forest classifier to detect spots which are not reproducible between subarrays. a random forest classifier is a machine learning procedure, which learns label categories from existing data and then applies them to new, previously unclassified data sets. it shows good predictive power while being robust to parameter setting and favorable theoretical properties. for the training of the classifier, we use the fact that each peptide is repeated on various subarrays. we label all subarray spots which are not within a 95% linear regression confidence band for a single linear regression on either of the remaining two subarray repeats as unreliable and all other spots as reliable. we use all columns of the result file besides the sequence as well as all columns of the quality control result file of a corresponding scatterlight scan for each spot as explanatory variable. we use random sampling to balance the set of unreliable and reliable spots in their size and then use the prediction of the random forest to label all spots of a new data set. for each peptide, we directly discard all those spots indicated as unreliable and only average over the remaining spots, weighted by their probability of being reliable. if all spots of one peptide are marked as unreliable, we exclude this peptide from any further analysis. mixture model for secondary antibody binding detection in peptide microarray experiments, secondary antibody binding peptides pose a major complication. when peptides show a direct reaction with the dye-coupled secondary antibody, two effects are confounded. it cannot be distinguished to which extent a measured intensity is explained by this undesired and unspecific reaction rather than the specific reaction of interest. thus, to avoid ambiguous results which hinder interpretation and may result in misleading conclusions, peptides showing significant unspecific reaction with the secondary antibody need to be identified and excluded from further analysis. reactivities on empty slides without any sample have been identified as important indicators. in experimental settings not allowing empty slides, this step can be excluded from the overall data analysis pipeline with an increased risk of false positives. however, this approach is tailored to identifying the unspecific binding of the secondary antibody to the peptide, not to the slide itself, which should affect all spots in a similar mode and thereby increase the noise level. this effect is addressed by the subsequent signal calling step. rather than relying on a fixed threshold, we propose a mixture model strategy to derive a probabilistic criterion for the identification of secondary antibody binding peptides from empty slides. for this purpose, we assume that the normalized intensities of all peptides on the empty slides stem from two populations: one population for peptides, which do not show any reaction to the secondary antibody besides noise and a second population for peptides which show a reaction to the secondary antibody. we model each population by a normal distribution and let the mean and variance of the first population be denoted by μ and and for the second population by μ and. a schematic display of a histogram of the measured intensities and the two underlying distributions in red and green is given in figure 1d. the orange line indicates the resulting cutoff. we apply an expectation maximization algorithm to iteratively adapt the estimates for μ μ, and, and for each measured normalized intensity, we update the likelihood that it belongs to either one of the two populations. for the initialization, we rely on the empirical mean and variance estimate for two groups of control peptides. one group for which we do not expect any reaction with the secondary antibody and a second group of peptides which are known to show various degrees of reactivity. we then use the resulting distribution of the peptides which do not show any reaction to the secondary antibody, to derive a criterion for secondary antibody binding. since we expect most peptides not to show any reaction, we expect this distribution to be more reliably estimated and mostly driven by random variations. we regard the upper 5% quantile of this distribution to be the cutoff. thus, for any measured normalized intensity on the empty slides, we compute any resulting value of below indicates only low likelihood for this intensity being this large by chance and thus the corresponding peptide is disregarded for the further analysis. with this cutoff generation we rely on the assumption that at least the population of peptides not showing a secondary antibody binding follow a normal distribution. however our historic slides at least approximately support this assumption. instead of this cutoff generation, it is also possible to directly compute a false discovery rate from the likelihood of the mixture model components, however this relies more heavily on the correctness of the normality assumption for the group of the reactive peptide than the proposed approach, which depends more on the not reactive peptides, in which we are more confident. mixture model for signal call to distinguish those peptide spots which carry a signal from those spots which only display noise, we again apply a mixture model. we assume that both, the signal-carrying peptide spots as well as the noise spots can be modeled by normal distributions with mean μsignal respectively μnoise and standard deviations respectively. a schematic display of a histogram of the measured intensities and the two underlying distributions in red and green is given in figure 1e; the orange line indicates the resulting cutoff. we initialize these parameters again based on suitable control classes. for one class of peptides expected to show a reaction, and for a second class of peptides not expected to show a reaction, we separately compute the empirical mean and variance and use these as starting values. we then update these values as well as the likelihood for each data point using the expectation maximization algorithm as described for the secondary antibody binding detection. we then use the distribution of the noise spots to identify the signal carrying spots which depart significantly from the noise distribution by identifying the upper 5% quantile of the noise distribution as a cutoff. thus, for any measured normalized intensity, we compute analog to equation any resulting value of below indicates only low likelihood for this intensity being this large by chance. experiments experiments are conducted on a total of peptide spots which were printed on six pepstar peptide microarrays with an identical design and printed in two batches. two separate batches were chosen to show the maximum possible inter-experiment variation. on each slide 10% of the total spots are used as dedicated control spots. 5% of all peptides are meant as positive controls, expected to show strong specific reactions, % are secondary antibody controls, which are expected to show strong reactions with the secondary antibody. 1% of all peptides are negative controls not expected to show any reaction and the remaining % are process controls. among the remaining peptide spots per array, spots contain peptides of human plac- and spots contain peptides of nyeso, vegfa, and cd. all proteins are represented by peptides with a length of amino acids and an overlap of amino acids. the slides were incubated in a hs pro hybridization station. one slide per print batch was incubated with diluent only, one slide per batch with a pooled plasma of six patients without indication of cancer and low concentration of a plac1-specific spike-in antibody and the remaining slide per batch with the same sample material, but high concentration of the plac1-specific spike-in antibody. the experimental setup is summarized in table a cy5-conjugated affinipure mouse anti-human igg was used as a secondary antibody. all slides were scanned with a genepix microarray scanner and converted to a result file using genepix pro, version. all result files are available from the authors' website. summary of the experimental setup for the six microarrays used in this study. the analysis of results is focused on five aspects: we evaluate the efficiency of the linear model in reducing the unexplained variance, the unreliable spot finding procedure by comparing the reproducibility of the two print batches before and after quality control, the removal of secondary antibody binding peptides by focusing on the nyeso, vegfa, and cd peptides, which are not expected to show any reaction with the given normal sera, with the exception of secondary binding effects, and the efficiency of the signal call by analyzing the specificity and sensitivity of our approach based on the nyeso, vegfa, and cd peptides which are expected to be non-reactive and the plac1-peptides, which we expect to show a reaction. we regard nyeso, vegfa, and cd peptides showing a reaction as false positives and those not showing a reaction as true negatives. similarly, we regard plac- peptides not showing a reaction as false negatives and those showing a reaction as true positives. sensitivity, specificity, and accuracy are then given by in a final step, we compare our approach to the approaches introduced by and based on the accuracy. to judge the significance of the results of experiments -, we computed 95% bootstrapped confidence intervals for all sensitivity, specificity and accuracy values based on times resampling from the peptide intensities of each class. in order not to confound the results of the unreliable spot finding step with the subsequent steps and to allow the fair comparison with the existing methods which do not identify unreliable spots, experiments - were run on all peptide spots, even though some were marked as unreliable.  we analyze the effect of the major steps of our approach including the linear model, the quality control and the secondary antibody binding detection and signal call as well as comparing the overall procedure on the antibody reactivity data. with regard to the linear model fit, a summary of results for print batch and the high spike-in antibody concentration is given in table it is evident that the peptide effect itself is the major influence of the measured intensity. still, all remaining explanatory variables, which model systematic effects, are highly significant as shown by the p-values for the corresponding f-test. row and column effects, so the spatial position of a peptide, show the overall strongest influence of the systematic effects. the needle and subarray effect are smaller, but still highly significant. however, it should be noted that a confounding of effects cannot be excluded, e.g. when considering a single needle only a subset of row and column positions are addressable on a subarray. including all interaction effects requires an unfeasibly high number of degrees of freedom, while choosing subsets does not significantly improve results. overall, the systematic effects contribute to reducing the unexplained variation in the measured peptide intensity data by 65%. linear model fit summary. the f-values and the corresponding probabilities ) clearly indicate that all explanatory variables used are highly significant and that each contributes to reducing the variance present in the arrays. the peptide sequence itself shows by far the strongest effect while requiring only twelve degrees of freedom, it shows a sum of squares for its effect of, much stronger than the residuals sum of squares of while not as strong as the peptide effect, the remaining explanatory effects reduce the residuals by 65%. column and row effects are strongest, but subarray and needle effects require fewer degrees of freedom, resulting in strong mean sum of squares values. after training on a total of historic slides with identical layout, the unreliable spot finding procedure marks % of all spots in print batch and % of all spots in print batch as unreliable for the high concentration experiment. since peptides are excluded only in cases when all three replicates for one peptide are marked as unreliable, this results in the removal of % of all peptides in print batch and % of all peptides in print batch visual inspection shows that these peptides primarily show large variation between the two data sets. by removing these peptide spots, we see an improvement of 3% in the coefficient of variation from r = to r =. still, not all outliers are detected by the approach. analysis of the missed outliers however indicates that these undetected outliers consistently affect replicates of only nine different peptide sequences, in contrast to the different peptide sequences marked as unreliable by the unreliable spot finding procedure. we evaluated our secondary antibody binding removal approach by computing the number of peptides which were excluded from the further analysis because of secondary antibody removal. out of non control peptides 354/ were identified with reaction on the empty slides and excluded from the further analysis. this affected none of the plac- peptides which were expected to show a specific reaction and 0/ of the nyeso, vegfa and cd peptides not expected to show a reaction. consequently, when regarding the sensitivity and specificity which were computed based on these peptides, we observe that with the secondary antibody binding removal the specificity is slightly increased while the sensitivity remains unchanged when compared to running our approach without this step. the sensitivity and specificity for our approach are shown in figure in comparison to the approaches of and. for two print batches, two slides were incubated with a low and a high concentration of spike-in antibody, respectively. overall, we observe for all approaches that sensitivity is more affected by the reduction of the antibody concentration than specificity, which is rather stable. in comparison, our approach shows best sensitivity across all settings. it shows a sensitivity of to compared to to for the other two approaches for the high antibody concentration data. this difference is more expressed for the low antibody concentration data with a sensitivity of to for our approach compared to to for the previous two approaches. specificity is slightly worse for our approaches for print batch with values of to for our approach and to for the existing approaches. similarly for print batch our approach features specificity of to compared to to. as a consequence of the substantially improved sensitivity, our approach also displays best accuracy. while this difference is small for the high antibody concentration with accuracy values of to for our approach compared to to for the previous approaches, it is much more substantial for low antibody concentration. here, our approach shows an accuracy ranging from to, compared to to for the two other approaches. rapmad is based on the assumption of normality, especially for the linear model fit and the mixture model. for all arrays incubated with plasma, we analyzed these assumptions and find support for approximate normality. discussion removal of systematic bias by applying a linear model which explains 65% of the noise variation in the data, rapmad significantly improves the reliability of peptide microarray experiments. in addition to previous findings, we establish that also the spatial position on the peptide microarray is of high impact for removing systematic effects. this is important for distinguishing low intensity signals from noise. since the effects are of a systematic nature, it also reduces the risk that systematic effects are mistakenly identified as peptide signal effects. despite the strong reduction in unexplained variation, not all occurring systematic effects are necessarily included in the model. for instance, we decided against including all interaction effects. while it can be of high interest to model the interaction of a specific peptide with a specific needle, there is usually not a sufficient number of replicates to allow a reliable estimation in a standard experiment. however, for large-scale experiments, it could be beneficial to design specific slides with sufficient replicates to quantify these effects. the effects of the secondary antibody are accounted for within subsequent analysis steps and are thus not modeled explicitly. improved data quality after removal of unreliable measurements and secondary binding peptides the production of peptide arrays as well as the binding of the secondary antibody involve sophisticated chemical processes. in spite of numerous experimental improvements and constant monitoring, artifacts resulting from errors in the spotting of peptides or the unspecific binding of the secondary antibody can strongly influence measurements and result in misleading interpretations. we have proposed two computational statistics driven approaches to identify and exclude these artifacts and show that they can be successfully applied. by marking approximately % of the peptide spots as unreliable, we can improve the reproducibility and thus the data quality. it should be noted however that still not all unreliable and irreproducible peptides are detected. the small number of peptide sequences being connected to the undetected outlying observations show that this entails a sequence-specific phenomenon rather than a print process artefact which would affect all sequences more evenly. while characteristics of the print process are captured within the available learning data, for instance in the scatterlight scans, sequence-specific effects are difficult to detect. in addition, the unreliable peptide spot finding procedure is limited by the available training data. we rely on labels based on the assumption that a specific replicate which departs from the remaining two replicates should be labeled as unreliable. we are aware that this definition is not necessarily correct since there are cases when the two remaining replicates are unreliable, but the specific replicate itself is correct. this can result in incorrect labels, which strongly reduce the quality of the random forest prediction. while manually validation could improve the label-quality, this reduces the applicability of the procedure to other data than our own. thus, we accept a possibly smaller impact of the quality control procedure to allow the fully automated adaption to other datasets. historic data for training is not always available and thus, this is an optional step for our pipeline. by identifying % of the spots as secondary antibody binding, we can improve the specificity of the procedure without sacrifying its sensitivity. while the impact on the data shown here is comparatively small, this step can be of key importance to avoid wrongful interpretation. when dealing with peptides prone to printing or synthesis problems or the binding of secondary antibody to peptides, any data analysis tool is challenged by the fact that the relevant information was lost before data acquisition, rendering the complete recovery of the information impossible. while not recovering the correct information, the magnitude of excluding wrongful results is significant. if wrongful results were not excluded, research could result in wrongful conclusions. with the removal of wrongful results no conclusions are possible for these doubtful peptide spots. since in most cases, only subsets of peptide measurements are affected and designs usually contain several replicates, repeats of experiments are often not required to come to reliable conclusions. strong sensitivity for low intensities a particular strength of rapmad is the fact that it is especially suitable for low intensity measurements such as those occurring for low antibody concentrations. while the problem of separating signal from noise is substantially more difficult when the signal intensity is reduced, we only see a comparatively minor reduction in accuracy when reducing spike-in antibody concentration. sensitivity is affected more significantly than specificity, but still at or above. this robustness of our approach can be explained by the adaptivity of the signal calling procedure, which is based on correctly estimating the noise distribution and generating a noise-threshold based on this distribution. thereby, it is less affected by changes in the signal intensity. as a consequence, the noise-threshold is adaptive and even low intensity signals can be detected to be above the threshold, resulting in a strong sensitivity. favorable comparison to existing approaches the strength of rapmad to maintain a high sensitivity even for low intensity settings is also what sets it apart from existing algorithms. thus, the sensitivity and accuracy of our approach in low antibody concentration experiment is competitive to the ones of the two other approaches for high antibody concentrations. this is a result of the adaptive procedures which are better suited to distinguish signal from noise in low intensity situation than preset thresholds and demonstrates that our computational approach is capable of compensating for experimental limitations. for high intensities as seen with high antibody concentration, it is generally easier to achieve high sensitivity due to a better separation of signal and noise in the raw data. that is also why we see less of a benefit of rapmad in comparison to existing algorithms. increasing the effectiveness of screening procedures peptide microarrays have been primarily used within screening settings with follow-up experiments focusing on a small number of promising candidates to ensure specificity. unidentified reactive peptides from a peptide microarray experiment cannot be distinguished from the large number of rightfully excluded unreactive peptides and are thus not further regarded. by substantially increasing the sensitivity for low intensity reactions, our approach avoids the wrongful exclusion of peptides and thereby increases the effectiveness of screening procedures. while standard settings were used in all experiments, rapmad additionally allows the user to adapt the quantiles used in the mixture models to trade additional sensitivity for specificity.  within this contribution, we introduce rapmad, a tool for the robust and rapid analysis of peptide microarray data. rapmad is an automated, multi-step approach that combines several computational statistics procedures to augment the data quality of peptide microarray and to allow a reliable analysis. its steps include the preprocessing of the data by removing systematic effects, the exclusion of unreliable measurements and secondary antibody binding peptides, and a probabilistic signal call for reactive peptides. in comparison with existing algorithms, it shows competitive performance for high antibody concentration. for low antibody concentration, it shows a significant increase in accuracy over existing approaches. this is mainly due to its substantially improved sensitivity at competitive specificity. with its increased sensitivity and its automated data analysis, rapmad can thereby contribute to establish and broaden the usage of peptide microarrays as a standard tool for a wide-range of peptidomics applications. authors' contributions byr, ml, ur, jcc, and us designed the methods, byr and ml wrote code and analyzed data. all authors contributed to the design of the experiments. yk and ar set up and ran the experiments. byr wrote the manuscript with contributions from ml and jcc. all authors read and approved the manuscript. supplementary material additional file 1 supplementary figures. pdf document containing supplementary figures. click here for file acknowledgements the authors would like to thank christina fleig for data acquisition, antonia masch for critical reading of the manuscript and ganymed ag, mainz for providing antibodies.
since data characteristics of peptide microarrays differ from dna oligonucleotide microarrays, computational methods need to be tailored to these specifications to allow a robust and automated data analysis. providing sensitivity is thus a primary goal of data analysis procedures. peptide microarrays offer an enormous potential as a screening tool for peptidomics experiments and have recently seen an increased field of application ranging from immunological studies to systems biology. we evaluated rapmad in antibody reactivity experiments for several thousand peptide spots and compared it to two existing algorithms for the analysis of peptide microarrays. while follow-up experiments can ensure the specificity of results, sensitivity cannot be recovered in later steps. particularly, it shows substantially improved sensitivity for low intensity settings without sacrificing specificity.
 peptide microarrays offer an enormous potential as a screening tool for peptidomics experiments and have recently seen an increased field of application ranging from immunological studies to systems biology. by allowing the parallel analysis of thousands of peptides in a single run they are suitable for high-throughput settings. since data characteristics of peptide microarrays differ from dna oligonucleotide microarrays, computational methods need to be tailored to these specifications to allow a robust and automated data analysis. while follow-up experiments can ensure the specificity of results, sensitivity cannot be recovered in later steps. providing sensitivity is thus a primary goal of data analysis procedures. to this end we created rapmad, a novel computational tool implemented in r.  we evaluated rapmad in antibody reactivity experiments for several thousand peptide spots and compared it to two existing algorithms for the analysis of peptide microarrays. rapmad displays competitive and superior behavior to existing software solutions. particularly, it shows substantially improved sensitivity for low intensity settings without sacrificing specificity. it thereby contributes to increasing the effectiveness of high throughput screening experiments.  rapmad allows the robust and sensitive, automated analysis of high-throughput peptide array data. the rapmad r-package as well as the data sets are available from http://www.tron-mz.de/compmed. 