BACKGROUND
in cancer, the path from normal to malignant cell involves multiple genomic alterations including losses and gains of genomic dna.
a long series of studies have demonstrated the biological and clinical relevance of studying such genomic alterations .
genome-wide scans of copy number alterations may be obtained with array-based comparative genomic hybridization , snp arrays and high-throughput sequencing .
after proper normalization and transformation of the raw signal intensities obtained from such technologies, the next step is usually to perform segmentation to identify regions of constant copy number.
many segmentation algorithms are designed to analyse samples individually , while most studies involve multiple samples, multiple tracks, or both.
joint handling of multiple samples is computationally and conceptually challenging, see e.g.
<cit> .
most systematic approaches for this problem are based on individual segmentation of each sample followed by post-processing to combine results across samples , while some recent publications propose strategies for joint segmentation of all samples  <cit> .
recently, the emergence of new technologies have pushed the limit of genomic resolution, opening new vistas for studying very short aberrations, including aberrations affecting only part of a gene or gene regulatory sites in the dna.
a major challenge raised by these novel technologies is the steadily growing length of the data tracks, which drastically increases the demand for computationally efficient algorithms.
the occurrence of extreme observations  of biological or technical origin pose an additional challenge, as most segmentation methods are substantially affected by such observations.
picard et al.
<cit>  propose a least squares based segmentation method that results in a piecewise constant fit to the copy number data.
their approach assumes that the user either supplies the desired number of segments or leaves to the method to automatically determine this number.
in this paper, we describe a related approach.
in particular, the proposed method utilizes penalized least squares regression to determine a piecewise constant fit to the data.
introducing a fixed penalty γ> <dig> for any difference in the fitted values of two neighboring observations induces an optimal solution of particular relevance to copy number data: a piecewise constant curve fully determined by the breakpoints and the average copy number values on each segment.
the user defined penalty γessentially controls the level of empirical evidence required to introduce a breakpoint.
given the number of breakpoints, the solution will be optimal in terms of least squares error.
to achieve high processing efficiency, dynamic programming is used .
to further increase computational efficiency, a novel vector based algorithm is proposed, and even further speed optimization is obtained through heuristics.
a central aim of the present work has been to provide methodology and high-performance algorithms for solving single- and multiple-track problems within a statistically and computationally unified framework.
all proposed algorithms are embedded in a comprehensive software suite for copy number segmentation and visualization, available as the bioconductor package copynumber.
main features of the package include:
• independent as well as joint segmentation of multiple samples
• segmentation of allele-specific snp array data
• preprocessing tools for outlier detection and handling, and missing value imputation.
• visualization tools
implementation
systems overview
the copynumber package provides functionality for many of the tasks typically encountered in copy number analysis: data preprocessing tools, segmentation methods for various analysis scenarios, and visualization tools.
figure  <dig> shows an overview of the typical work flow.
input is normalized and log2-transformed copy number measurements from one or more acgh, snp-array or hts experiments.
allele-frequencies may also be specified for the segmentation of snp-array data.
it is strongly recommended to detect and appropriately modify extreme observations  prior to segmentation, as these can have a substantial negative effect on the analysis.
for this purpose, a specially designed winsorization method is included in the software package.
a missing-value imputation method appropriate for copy number data is also available.
segmentation methods for three different scenarios  are implemented in the package.
all these methods are referred to as piecewise constant fitting  algorithms and seek to minimize a penalized least squares criterion.
in single sample pcf, individual segmentation curves are fitted to each sample.
in multi-sample pcf, segmentation curves with common segment borders are simultaneously fitted to all samples.
in allele-specific pcf, the segmentation curves are fitted to bivariate snp-array data, providing identical segment borders for both data tracks.
a set of graphical tools are also available in the package to visualize data and segmentation results, and to plot aberration frequencies and heatmaps.
also included are diagnostics to explore different trade-offs between goodness-of-fit and parsimony in terms of the number of segments.
in the remaining part of this section, a formal description of the algorithms is given.
however, note that these details are not a prerequisite for reading later sections or for using the copynumber package.
preprocessing: outlier handling
a challenging factor in copy number analysis is the frequent occurrence of outliers - single probe values that differ markedly from their neighbors.
such extreme observations can be due to the presence of very short segments of dna with deviant copy numbers, to technical aberrations, or a combination.
when identification of cnvs is a purpose of the study, the multi-sample method described below may be applied for such detection.
however, when the focus is on detection of broader aberrations, the potentially harmful effect of extreme observations on aberration detection methods induces a need for outlier handling procedures .
since the copynumber package is based on least squares, an extreme observation will tend to cause the detection of a short segment.
when searching for broader segments, such short  segments will represent noise and may also affect the identification of other segments.
we therefore now describe a procedure for reducing the effect of extreme observations, while the effects of this method will be considered in the results and discussion section.
winsorization is a simple transformation reducing the influence of outliers by moving observations outside a certain fractile in the distribution to that fractile .
for identically distributed observations y <dig> …,yp, the corresponding winsorized observations are defined as yjw=Ψ where
Ψ=Ψ=−θ,y<−θθ,y>θy,otherwise.
here, θ> <dig> determines how extreme an observation must be to be relocated, as well as the replacement value.
a common choice is θ=τs, where typically τ∈ and s is a robust estimate of the standard deviation .
a robust scale estimator is the median absolute deviation , defined as the median of the values |yj−m^|, where m^ is the median of y <dig> …,yp.
for normally distributed observations, sm= <dig> ·mad corresponds to sd.
winsorization of copy number data may be achieved by first estimating the trend in the data and then winsorizing the residuals.
let the observations representing copy numbers in p genomic loci be y=, ordered according to genomic position.
a simple estimator of the trend is the median filter.
the trend estimate m^j in the jth locus is then given by the median of yj−k,…,yj + k for some k> <dig>  e.g.
k= <dig>  the sd of the residuals yj−m^j may then be estimated with the mad estimator sm, and winsorized observations y1w,…,ypw obtained by yjw=m^j+Ψ.
often, such simple and fast winsorization is sufficient.
however, copynumber also includes an iterative procedure with improved trend estimation based on the segmentation procedures described below .
single sample segmentation
consider first the basic problem of obtaining individual segmentations for each of a number of samples.
suppose attention is restricted to one chromosome arm on one sample.
for each of the p loci, the obtained measurement can be conceived of as a sum of two contributions:
yj=zj+εj
where zjis an unknown parameter reflecting the actual amount of sample dna at the j’th locus and εjrepresents measurement noise.
a breakpoint is said to occur between probe j and j +  <dig> if zj≠zj +  <dig>  the sequence z <dig> …,zp thus implies a segmentation s={i <dig> …,im} of the chromosome arm, where i <dig> consists of the probes before the first breakpoint, i <dig> consists of the subsequent probes until the second breakpoint, and so on.
to fit model , we minimize the penalized least squares criterion
∑j=1p2+γ·|s|
with respect to the sequence z <dig> …,zp.
here, |s| denotes the number of segments in s, and γ> <dig> is a constant that controls the trade-off between seeking a good fit to the data  and restraining the number of level shifts .
the minimizer z^ <dig> …,z^p of  is fully determined by the segmentation s, since the best fit z^j on a given segment i is the average y¯i of the observations on that segment.
substituting the latter into  we obtain the equivalent criterion:
l=∑i∈s∑j∈i2+γ·|s|
=∑i∈s∑j∈iyj2−∑i∈s2/ni+γ·|s|
where ni denotes the number of probes in segment i. note that the first term in  does not depend on the segmentation s, hence minimization of  is equivalent to minimizing
l′=−∑i∈s2/ni+γ·|s|.
naive optimization of the cost function  with respect to the segmentation s requires examination of every possible division of the probes on a chromosome arm into segments.
for large p, this is not practically feasible.
however, a much more efficient implementation based on dynamic programming and requiring only o operations is available.
dynamic programming is a method for solving complex problems by breaking them down into simpler subproblems, and specifically for problems where global decisions can be decomposed into a series of nested smaller decision problems.
the crucial observation that allows the use of dynamic programming to solve the present segmentation problem is that the optimal segmentations on each side of a breakpoint are mutually independent.
this can be used to iteratively build up a solution to the global segmentation problem.
suppose we know the optimal segmentations from the first probe up until the st probe.
assume furthermore that the optimal segmentation for the k first probes contains breakpoints.
then the optimal segmentations from the last of these breakpoints and downwards has already been computed.
thus, by solving the above subproblems iteratively for increasing k, each step can utilize the results from the previous steps .
more formally, assume that the optimal segmentation of 1…r and the corresponding total error er are known for all probes r<k.
to extend the solution to r=k, first note that there must be a last segment starting at some index j≤k.
from  we find that the cost term associated with that segment is:
djk=1j−k−1∑r=jkyr <dig>
then the total error for the optimal solution up until index k is found by minimizing the cost over the possible start positions j of the last segment.
this cost consists of three terms: the cost of the last segment , the optimal cost of the segmentation up until that point  and the penalty for the break point :
ek=minj∈{ <dig> …,k}
where e0= <dig>  the main work load of the above computation is to determine djk for all 1≤j≤k≤p.
in interpreted languages  where loop execution is often quite inefficient, a considerable improvement of performance may be obtained by utilizing native-language vector operations.
let ajk=∑r=jkyr, ak= and dk=.
then we may calculate all required coefficients through a simple recursion:
ak=+ykdk=−ak∗ak/
where = and operators are vector-based.
hence, addition of a vector and a scalar adds the latter to each component of the former, and multiplications and divisions are performed component-wise on the operands, e.g., ak∗ak=.
algorithm  <dig> summarizes the computations.
algorithm 1: single sample pcf
input: log-transformed copy numbers y <dig> …,yp; penalty γ> <dig> output: segment start indices s <dig> …,sm and segment averages y¯ <dig> …,y¯m.
<dig>  calculate scores by letting a0= and e0= <dig>  and iterate for k=1…p:
•
ak= + yk
•
dk=−ak∗ak/
•ek=
storing also the index tk∈{ <dig> ,…,k} at which the minimum in the last step is achieved.
<dig>  find segment start indices  s1=tp,s2=ts1−1…,sm= <dig>  where m≥ <dig>
<dig>  find segment averages y¯m=ave−1) for m= <dig> …,m, where s0=p +  <dig>
throughout the paper we will tacitly assume that the penalty for the ith sample is γi=γσ^i <dig>  where σ^i <dig> is the estimated sample specific residual variance.
in this way, we avoid scale dependency, and obtain consistent results for samples with equal signal-to-noise ratios.
such rescaling is also done by default in copynumber.
note that replacing the data yji for the ith sample with yji/σ^i for j= <dig> …,p, and rescaling after estimation, has the same effect.
in copynumber, the algorithm has also been extended to allow a constraint on the least number of probes in a segment.
multi-sample segmentation
detection of very short or very low amplitude segments requires a small penalty γ, with low specificity as a potential result.
however, when such segments are common to several samples, joint segmentation of multiple samples is an additional mechanism to increase sensitivity.
this is a main motivation for introducing multi-sample segmentation methods that impose common breakpoints across all samples.
such methods are potentially useful for discovery of copy number variations  and in those instances where the origin of the samples implies that segment boundaries are partly shared.
multi-sample segmentation with high penalty on breakpoints may also be used to obtain low-dimensional descriptions of the data, which may form the basis for defining variables to be used in statistical procedures relating aberration patterns to clinical outcome.
in the following, we describe a direct generalization of single sample pcf to handle multiple samples simultaneously, obtaining common breakpoints for all the samples with minimal residual sum of squares for a given number of breakpoints.
suppose copy number measurements yi= for samples i= <dig> ,…,nare obtained at the same loci in each sample.
by direct generalization of the criterion , we seek in multi-sample pcf the minimizer of
l=∑i=1nl
where l is defined as in  and s is a given segmentation common to all samples.
algorithm 2: multi-sample pcf
input: log-transformed copy numbers for n samples y <dig> …,yp∈rn; penalty γ> <dig> output: common segment start indices s <dig> …,sm and segment averages y¯ <dig> …,y¯m∈rn.
<dig>  calculate scores by letting a0= and e0= <dig>  and iterate for k=1…p:
•
ak= + yk
•
dk=−1t/
•ek=
storing also the index tk∈{ <dig> ,…,k} at which the minimum in the last step is achieved.
<dig>  find segment start indices  s1=tp,s2=ts1−1…,sm= <dig>  where m≥ <dig>
<dig>  find segment averages y¯m=ave−1) for m= <dig> …,m, where s0=p +  <dig>
the multi-sample pcf algorithm  is in principle quite similar to single sample pcf.
however, when updating the solution from k− <dig> to k, the sums and sums of squares for the segments must be accumulated and stored separately for each sample.
this can still be done iteratively, implying that the computational effort will be approximately equal to carrying out single sample pcf on the same set of samples.
since the noise level may vary between samples, normalisation of the samples prior to segmentation and corresponding rescaling after estimation is advisable.
it may also be desirable to scale the samples, e.g.
to adjust for different tumor percentages.
thus, prior to running multi-sample pcf, we may replace yi by wiyi/σ^i for i= <dig> …,n, where wi are weights and σ^i is an estimate of the sd.
in copynumber normalization is performed by default for multi-sample pcf while further weighting is left as an option for the user.
allele-specific segmentation
the pcf algorithm is easily adapted to variants of the basic segmentation problem discussed above.
here, we consider an adaptation to handle snp genotype data.
we then have for each snp locus a measurement of  copy number  as well as the b allele frequency .
we may also have measurements of copy number only for a number of additional loci.
the b allele frequency is a number between  <dig> and  <dig> indicating the allelic imbalance of a snp.
for a homozygous locus we have baf close to  <dig> or  <dig>  while for a heterozygous locus with an equal number of the two alleles a and b, baf will be close to  <dig> .
an imbalance between the number of a’s and b’s results in a baf value deviating from  <dig> .
a change in the total number of copies of a segment will alter the logr value, hence result in a level shift in the logr track.
unless the copy number change is balanced with respect to the two alleles, the baf value will also change.
in cases involving multiple copy number events at the same locus, the change may manifest itself only in one of the two tracks.
for example, a loss of one copy of a followed by a gain of one copy of b would lead to unchanged logr and changed baf.
the purpose of the allele-specific pcf algorithm is to detect breakpoints for all such events.
it fits piecewise constant curves simultaneously to the logr and the baf data, forcing breakpoints to occur at the same positions in both.
we emphasize that the purpose of the allele-specific pcf algorithm is segmentation only and not to make allele-specific copy number calls.
however, such calls can be made on the basis of the segmentation described below, and this is done e.g.
in the ascat algorithm  which estimates allele-specific copy numbers as well as the percentage of cells with aberrant dna and the tumor ploidy  <cit> .
suppose the data are given by  for j= <dig> …,p, where rj denotes the logr value and bj the baf value at the jth locus.
for copy number probes, only rj is given and bj will be missing .
for germline homozygous probes, the baf values are noninformative and should be omitted from the analysis.
if the germline genotype is known , the user should replace the corresponding baf values by na.
if the genotype is not known, the algorithm will apply a proxy to handle this issue .
prior to segmentation, the allele-specific pcf algorithm performs the following steps:
the baf data are mirrored around  <dig>  by replacing bjwith 1−bjif bj> <dig> .
baf values bj<θare replaced by na.
by default θ= <dig> .
if germline homozygous probes have previously been replaced by na’s, let θ= <dig>
let b~ <dig> …,b~m denote the nonmissing b allele frequencies.
corresponding copy number values r~ <dig> …,r~m are found by pairing each logr probe with the nearest b-allele probe  and then averaging logr values paired to the same b-allele probe.
finally, let y1= and y2=.
the remaining part of the allele-specific pcf algorithm is then essentially an adaptation of the multi-sample pcf algorithm applied to two samples.
it finds a common segmentation s for the two tracks by minimizing the penalized criterion
l=l+l
where l is defined as in .
fast implementations of pcf
the pcf algorithms may be generalized to allow breakpoints only at certain prespecified positions.
combined with simple heuristics, this may be used to further enhance the computational speed of pcf.
for brevity we describe only the single sample segmentation case here; however the copynumber package contains fast implementations of both single- and multi-sample pcf.
computationally inexpensive methods can be used to identify a set of potential breakpoints among which the breakpoints of the solution to  are highly likely to be found.
suppose we restrict our attention to such a set of potential breakpoints.
all relevant information for solving the optimization problem in  may then be condensed into three arrays containing the number of observations between two potential breakpoints, the corresponding sum of the observations and the sum of squares.
based on these quantities, pcf may be used with straightforward modifications.
since the algorithm is of order o, where q is the number of potential breakpoints, the potential increase in speed is substantial.
algorithm  <dig> outlines the procedure, while possible heuristics for finding potential breakpoints are discussed below.
one way to identify potential breakpoints is to use high-pass filters, i.e.
a filter obtaining high absolute values when passing over a breakpoint.
the simplest such filter uses for each position i the difference ∑j=i+1i+kyj−∑j=i−k+1iyj for some k. to reduce artifacts due to the abrupt edges of such a filter, the copynumber implementation assigns half weight to the outer 1/ <dig> of the observations on each side.
fast implementations of such filters in r may be obtained using the cumsum function.
we currently use two filters with k= <dig> and  <dig>  respectively; additionally the single sample pcf implementation includes a filter searching for aberrations of length equal to the lowest accepted one.
these filters together identify about 15% of the probe positions as potential breakpoints.
an additional way to speed up the computations on long sequences is to initially divide the sequence into overlapping subsequences, and iteratively find the solution.
having found the solution for the m first subsequences, we use high-pass filters to detect potential breakpoints for subsequence m +  <dig>  and then use the fast pcf algorithm with the latter potential breakpoints as well as those found by pcf on earlier subsequences.
the intention behind this iterative approach is to reduce potential boundary effects.
due to the quadratic order of the algorithm, this division into subsequences implies a substantial efficiency gain.
in copynumber, subsequences are used when the chromosomal arm length exceeds  <dig> probes, with subsequences of length  <dig> and overlap  <dig>
algorithm 3: fast pcf
input: log-transformed copy numbers y <dig> …,yp; penalty γ> <dig> output: segment start indices s <dig> …,sm and segment averages y¯ <dig> …,y¯m.
<dig>  apply heuristics to find potential breakpoints r <dig> r <dig> …,rq, where r0= <dig> and rq=p +  <dig>
<dig>  form aggregates by letting uk=∑j=rk−1rk−1yj, where k= <dig> …,q.
<dig>  calculate scores by letting a0=, c0=, e0= <dig>  and iterate for k= <dig> …,q:
•
ak= + uk
•
ck= + rk−rk−1
•
dk=−ak∗ak/ck
•
ek=
storing also the index tk∈{ <dig> ,…,k} at which the minimum in the last step is achieved.
<dig>  find segment start indices  s1=rtq,s2=rts1−1…,sm= <dig>  where m≥ <dig>
<dig>  find segment averages y¯m=ave−1) for m= <dig> …,m, where s0=p +  <dig>
RESULTS
selection of penalty
the selection of parameters determining the trade-off between high sensitivity  and high specificity  is important in all segmentation procedures.
in pcf, this is controlled by the single penalty parameter γ. a number of general model selection criteria exist, such as cross-validation, the akaike information criterion  and the related schwarz’s bayesian information criterion .
however, model selection for copy number segmentation is complicated by several factors.
first, the distribution of the data at hand may vary substantially.
an important example is the presence of local trends mimicking smaller aberrations; such low-amplitude “waves” in the data may e.g.
be due to variations in gc-content .
second, the purpose of the analysis may favor either higher sensitivity or higher specificity.
for example, in clinical studies aimed at finding prognostic markers, the main focus may be on the most pronounced and commonly occurring deviations, while detecting more sporadic aberrations may simply increase the noise level.
in our experience, the above model selection criteria tend to give too small penalty estimates and thus undersmooth the data.
this is consistent with previous investigations showing that aic and bic are not appropriate for the breakpoint problem .
simulation studies of specificity may suggest a lower bound on the penalty γ. for this purpose, sequences of independent and normally distributed observations without underlying aberrations were generated, and pcf was applied with different choices of γ. at γ= <dig> the number of falsely called aberrations is about  <dig>  per  <dig>  probes, at γ= <dig> roughly  <dig> per  <dig>  probes, at γ= <dig> roughly  <dig> per  <dig>  probes, and for γ≤ <dig> the number of falsely called aberrations is substantial.
this suggests γ≈8− <dig> as a lower bound.
since the number of false aberrations per chromosome increases with increasing probe density, low values are most relevant for arrays with low probe density.
in the presence of local trends, the number of false calls tends to inflate and the penalty should thus be increased above the lower bound.
a fairly conservative penalty of γ= <dig> is the default in the copynumber package.
this provides a starting point for exploration of the best penalty value for the specific problem at hand, however a systematic inspection of results obtained for different penalties is advisable.
figure  <dig> illustrates the effect of changing γ. notice that the main features in the data are captured across the whole range of γ-values, while finer details are only evident for smaller values.
aberration calling
aberration calling is used for detection of recurring alterations and in many other analyses.
introducing a parameter θ> <dig> that determines the sensitivity of the aberration calling , we call probes for which z^j<−θ as losses and probes for which z^j>θ as gains.
optionally, different thresholds θ + and θ− may be used for gains and losses.
to examine how well pcf aberration calling manages to distinguish between normal and aberrant regions, performance was compared with a very accurate measurement method.
specifically, aberration calls obtained with pcf on the basis of  <dig> m snp array data on  <dig> samples were compared with calls obtained with mlpa .
since mlpa is limited to a small set of genomic positions, only  <dig> loci were used for the comparison.
in all samples combined, mlpa identified  <dig> aberrant and  <dig> normal loci .
using the mlpa-classification as the gold standard, the sensitivity and specificity of pcf aberration calling were calculated for a range of threshold values θ. figure  <dig> shows the resulting roc curves, and panel  illustrates how the results for pcf depend on the choice of γ. importantly, aberration calling appears to be only moderately dependent on the choice of parameter values over a fairly wide range of γ-values.
single- versus multi-sample segmentation
whether the initial segmentation of a dataset is most appropriately done using single- or multi-sample methods depends both on the purpose and the data.
using methods with common breakpoints for samples will increase the power for detecting concordant but quantitatively weak segments, while it will reduce the ability of detecting  discrepant breakpoints.
a well known example of aberrations with common boundaries is germline copy number variants , thus some proposed algorithms for cnv detection utilize segmentation with joint segment borders .
another important example of samples with  common segment boundaries arises when the samples originate from different clones of the same  tumor.
this is illustrated below in two examples, one on disseminated tumor cells from breast carcinomas, the other on tumor clones found at successive biopsies from lymphoma patients.
recent reports  <cit>  on marked variations in aberration patterns within the same tumor is likely to increase the number of studies using several samples taken from each tumor.
what is common as well as what differs in the aberration patterns will then be of interest, motivating the combined use of single- and multi-sample methods.
in applications searching for genomic copy number hot spots with relevance to cancer development, it may be important to utilize the precise delineation of the aberrations found in each sample, and thus the use of single-sample methods is most appropriate.
the identification of the relevant recurrent aberrations may then utilize post processing tools like gistic  <cit> , kcsmart  <cit>  or cghmcr  <cit>  .
if focus is on clustering samples or on constructing regression variables for relating more broad aberrations to clinical outcome, one may consider using multi-sample methods.
however, to be useful, the estimates from the multi-sample methods should in a proper way reflect the main information content in each sample.
this implies that a multi-sample analysis should result in estimates approximating those obtained from single-sample analyses.
figure  <dig> shows heatmaps of results from single- and multi-sample pcf for  <dig> breast cancers from the so-called micma data set  analyzed on 244k agilent arrays.
the main features appear to be well reflected in the multi-sample analysis.
on a more detailed level, differences can be observed: the multi-sample solution misses some short aberrations occurring in only a few samples, aberration borders are sometimes slightly shifted, and longer segments obtained with single sample pcf are often divided into subsegments with slightly different copy number estimates.
the moderate difference between the results of single- and multi-sample pcf was also confirmed by a comparison of the ability to detect specific aberrations as revealed by comparison to mlpa analyses, see figure 3b.
this indicates that at least for cancer types where aberrations are focused in certain areas of the genome, methods using joint boundaries might be considered for constructing variables to be used in further statistical analysis.
comparing tracks: analysis of disseminated tumor cells
disseminated tumor cells  are detected in the bone marrow of some patients with breast carcinomas.
the presence of dtcs in the bone marrow identifies patients with less favorable outcome , and genomic characterization of such cells is of substantial interest.
it is still an open question to what extent the aberration patterns in dtcs correspond to those found in the primary tumor; the dtcs may potentially have obtained new aberrations or, alternatively, the cells may have originated from  subclones of the tumor with less aberrations.
it is possible to analyze single cells using acgh; however, currently the noise level is high, making it difficult to draw definitive conclusions from a single cell.
however, since segment boundaries are assumed to be partly common, we tested the multi-sample pcf algorithm on breast cancer cases from which dtcs were available .
figure  <dig> shows the results on a set of dtcs and the corresponding primary tumor from one such patient.
since multi-sample pcf is used, segment boundaries are common, while the estimated level in each segment is determined by the individual dtc/primary tumor.
in figure  <dig>  two of the single cells seem to have a pattern similar to the primary tumor.
the last one has an essentially flat  profile and is likely to be a hematopoietic cell misclassified as a tumor cell .
these data thus indicate that the aberration pattern of the dtcs quite closely reflect that of the primary tumor.
with only two single cells present, figure  <dig> primarily shows that dtcs inherit the aberrations of the primary tumor; with higher numbers of cells, multi-sample pcf may also be used to search for aberrations found in dtcs but not in the tumor.
defining variables: genetic evolution in follicular lymphoma
follicular lymphoma is normally a slowly progressing malignancy, but relapses are common and the disease is usually fatal.
in a recent study,  <dig> biopsies from  <dig> patients diagnosed with follicular lymphoma were evaluated using a custom-made acgh platform consisting of 3k bac/pac probes  <cit> .
a whole-genome view of aberration frequencies  and highly correlated aberrations  are shown in figure  <dig>
although the delineation of segments varied between biopsies, several areas with a high frequency of aberrations could be detected.
to try to identify aberrations with prognostic potential, we therefore found a common segmentation for the initial biopsies taken from each of the  <dig> patients using the multi-sample pcf algorithm.
removing very low variance segments,  <dig> segments remained.
the corresponding copy number estimates were used as covariates in a multivariate cox proportional hazards regression.
this revealed  <dig> segments for which gains were significantly associated with a survival disadvantage.
a particularly strong association was detected for gains on chromosome x in male patients.
to study the relation between successive biopsies taken from the same patient, multi-sample pcf was applied to each patient individually .
as expected, many aberrations are common, but interestingly, some aberrations are present in early biopsies and not in later ones.
this contradicts the hypothesis of linear development which states that late tumor clones arise from earlier ones, and supports the alternative hypothesis of parallel evolution in different lymph nodes.
allele-specific copy number analysis in breast cancer
copy number alterations have been extensively studied in breast cancer.
to what degree gains and losses are associated only with certain alleles has been less studied.
in a recent study, genotyping of  <dig> breast carcinoma samples was performed using illumina 109k snp arrays, and the ascat method was used to infer the allele-specific copy numbers at each locus  <cit> .
however, to do this we first had to segment the data; for this purpose we applied allele-specific pcf segmentation to all samples.
in figure  <dig>  the result of this segmentation is shown for one particular sample and two different chromosomes.
in figure 7a, the segmentation of chromosome  <dig> is shown, and we clearly identify three segments on the p-arm with copy numbers less than two, larger than two, and identical to two .
suppose we consider only germline heterozygous loci, in which case the allelic ratio is 1/ <dig> when no aberrations are present .
the baf track reveals allelic imbalance in the first two segments, and more pronounced in the first segment than in the second.
this is consistent with a loss of one copy in the first segment , and a gain of one copy in the second segment .
the third segment has an allelic ratio of 1/ <dig>  notice that in case of allelic imbalance, the observed allele ratio is substantially closer to  <dig>  than expected by the above theoretical ratios.
this attenuation of the signal  is due to technical issues like cross-hybridization, as well as the fact that in reality the tumor is a mixture of cells with normal dna  and tumor cells with aberrant dna.
in figure 7b we notice a sharp trough in the logr track on 17p, accompanied by an allelic ratio close to  <dig> .
if for a certain snp locus one allele is substantially more frequently gained than the other allele, one may hypothesize that the former allele is subject to a larger selective pressure to change copy number.
this, in turn, may be an indication of different roles being played by the two alleles with respect to cancer progression and evolution, suggesting that loci subject to allelic skewness can be potential unique markers for breast cancer development.
even from a relatively small number of samples, probes with highly significant allelic skewness have been identified in a genome-wide statistical evaluation  <cit> .
outliers and winsorization
while least squares methods are often favored due to their optimality properties, they are also known to be sensitive to extreme observations.
thus, except if the purpose is to search for short aberrations of biological origins , we advise the use of an outlier handling procedure.
to evaluate the proposed winsorization scheme, we first established a suitable way of simulating extreme observations.
a classical way is to use “contaminated normals”, where the error distribution is a mixture of two normal distributions  <cit> .
with probability 1−α the error is drawn from a distribution n, and with probability αfrom n, typically with d= <dig> and α= <dig> .
we compared the fraction of outliers in observed copy number data to the corresponding fractions in normals and contaminated normals, using mad to estimate sds.
for the normal distribution, the fraction of observations outside  <dig> sd is  <dig> % and outside  <dig> sd <  <dig> %, while these fractions for the 5% contaminated normal are  <dig> % and  <dig> %.
for the agilent 244k used on the micma dataset, the fractions were  <dig> % and  <dig> %, that is, slightly above the values for the contaminated normal.
for the 44k agilent and illumina 109k applied to the same data, the percentages were slightly lower , however still indicating that the 5% contaminated normal is an appropriate distribution when evaluating robustness of copy number assessment procedures.
inspection of data obtained by the 318k illumina,  <dig> x 180k agilent and nimblegen  <dig> m arrays also confirmed the existence of substantial amounts of outliers.
the pcf algorithm was tested with and without winsorization on simulated data with outliers .
outliers in the contaminated distributions may cause the detection of short false aberrations; such spikes occurred roughly ten times per  <dig> probes, as compared to less than two times per  <dig> for uncontaminated data.
table  <dig> further shows that winsorization efficiently reduces the number of falsely detected aberrations and make results for the contaminated distribution roughly equal to the ones for the normal distribution.
in line with these observations, outliers tend to change the form of aberrations , while winsorization brings the distribution fairly close to the one found for normal data .
shown is the effect of winsorization on simulated data with outliers and artificial  aberrations.
two types of aberrations are considered:  aberrations of height  <dig>  and length  <dig> probes and  aberrations of height  <dig>  and length  <dig>  the contamination consists of normals with sd= <dig> and the mad estimate of sd equals  <dig> .
sensitivity is the percentage of amplified probes that are detected as amplified, while specificity is the percentage of non-amplified probes classified as such.
the false aberration column gives the percentage of aberrations not covering the central part of the real amplifications.
another way to avoid that a few extreme observations result in a segment is to impose a lower limit on the length  of a segment.
with a lower length limit of five probes, we found about twice as many false spikes as with winsorization when adjusting γto give equal sensitivity for true aberrations.
still, simulations indicate that a lower limit on segment length is valuable in combination with winsorization.
note that outliers of biological origin will be more extreme if the technology has an inherent low noise level, as is the case, e.g., for bac arrays and for high throughput sequencing.
thus, outliers are not a sign of inappropriate functioning of a technique, but a characteristic of the data requiring consideration in the analysis.
in summary, copy number data tend to contain a high fraction of outliers.
these outliers often induce false aberrations, but simple procedures like winsorization will efficiently reduce these undesired effects.
computational performance
in r, using the vector based pcf implementation described in algorithm  <dig> implies a substantial efficiency gain over loop based implementation, roughly a 10- <dig> times reduction in time requirements.
the fast implementation of pcf  gives a further marked reduction in computing time.
on the micma 244k dataset , the implemented fast version is about  <dig> times faster than the exact one, and uses around  <dig>  minutes to process the  <dig> samples .
the multi-sample method was slightly faster than the single sample version.
the average computation time  per sample is shown for copynumber , dnacopy v <dig> .
<dig>  and cghflasso v <dig> - <dig>  on the micma  <dig> k data set  and on the logr values from an illumina  <dig>  m snp array data set .
the iqr over samples is given in parenthesis.
the methods were applied to both raw data and data with outliers removed using the winsorization method.
all tests were performed on a pc with a  <dig> ghz intel i <dig> cpu with  <dig> gb of memory running windows  <dig> and r  <dig> .
<dig> .
the deviations between the solutions found by the exact pcf and fast pcf on the micma set were small; in terms of reduction in variance  below  <dig> %.
the differences observed for the curves were typically small shifts in the border of aberrations.
thus, we conclude that the results from the fast procedure for practical purposes may be regarded as global solutions to , and the fast version is therefore used by default in copynumber.
we also compared the performance of pcf with two other segmentation methods: circular binary segmentation   <cit>  and fused lasso regression   <cit> .
in comparison studies  <cit> , cbs has shown good performance in terms of sensitivity and false discovery rate.
it is probably the most commonly used freely available algorithm and is also implemented in several commercial analysis tools.
cbs is available in the r package dnacopy, which is used for this comparison.
fl is a more recent proposal implemented in the r package cghflasso, and is one of three preferred methods in the web-based segmentation tool cghweb  <cit> .
using default parameter settings, we compared the computing times of pcf, cbs and fl on the  <dig> samples in the  <dig> k micma data set, and on  <dig> samples from a  <dig>  m illumina snp array .
table  <dig> gives the average computation time  per sample.
with no preprocessing of the data, pcf is on average 3- <dig> times faster than cbs on both data sets, and about  <dig> times faster than fl on the largest data set.
note that copynumber detects and operates on chromosome arms, while dnacopy operates on whole chromosomes.
this partly explains the difference in performance between pcf and cbs for the micma data set; for the illumina data this has little impact due to the iterative approach used in pcf for the longest sequences.
pcf was also markedly faster in evaluations based on simulated data; however, comparisons are complicated by the fact that the speed of cbs depends on the data in a nontrivial manner.
as seen from the iqrs listed in parentheses in table  <dig>  the speed of cbs is quite variable from sample to sample while pcf and fl is nearly constant.
moreover, the table shows that cbs runs 2- <dig> times slower when outliers have been removed using winsorization, underlining that the performance of cbs is highly data dependent.
we underline that the above-mentioned results only relate to the current r implementations.
as mentioned in the introduction, pcf is conceptually similar to the cghseg method described by picard et al.
<cit> , and we also examined the computational performance of this method using the implementation in the r package cghseg.
using the version of cghseg that requires a prespecified number of segments for each chromosome, the algorithm is fast, although the speed depends on the number of segments.
using the full cghseg algorithm that automatically determines the number of segments, the algorithm is very slow for high-resolution data.
hence, making a fair comparison between pcf and cghseg is difficult.
segmentation accuracy
we further compare the accuracy of the segmentation solutions found by pcf and cbs.
figure 3c shows roc curves using mlpa classifications as the truth, and then applying a range of aberration calling thresholds to pcf estimates, cbs estimates, a running median with window size  <dig> and raw copy number data .
results for pcf and cbs are similar, both achieving high sensitivity and specificity.
the running median also gives good results, illustrating that many probes are fairly easy to classify and that the gain obtained by using methods like cbs and pcf is mainly an improved classification close to borders between segments.
we also repeated the simulation study in  <cit>  where cbs was found to be the most sensitive method while also having the lowest false discovery rate.
again, we found that pcf and cbs had very similar performance .
a more detailed comparison of segmentation results shows that overall results are quite similar for single sample pcf and cbs, however for both methods the results depend on the choice of parameter values and the handling of extreme observations, see additional file  <dig>  in conclusion, pcf and cbs typically provide similar results and have equivalent accuracy when parameters are tuned appropriately.
CONCLUSIONS
copy number segmentation based on least squares principles and combined with a suitable penalization scheme is appealing, since the solution will be optimal in a least squares sense for a given number of breakpoints.
we have proposed a suite of platform independent algorithms based on this principle for independent as well as joint segmentation of copy number data.
the algorithms perform similarly as other leading segmentation methods in terms of sensitivity and specificity.
furthermore, the proposed algorithms are easy to generalize and are computationally very efficient also on high-resolution data.
the bioconductor package copynumber offers a user-friendly interface to the proposed algorithms.
several extensions and modifications of the proposed least-squares framework are possible.
in principle, the l2-based distance measure used in the current implementation of pcf is easily extended to general lp-distances.
however the current implementation is highly optimized for l <dig>  and other distance measures would require substantial heuristics to obtain comparable computational performance.
another extension is to introduce locus specific penalties for breakpoints, thus essentially introducing a prior on the location of breakpoints.
work in progress includes specialized routines to handle high throughput sequencing data more efficiently and joint analysis of multiple samples in allele-specific pcf.
availability and requirements
project name: copynumber
project home page:http://heim.ifi.uio.no/bioinf/projects/copynumber/
operating system: all systems supporting the r environment
programming language: r
other requirements: no
license: gnu artistic license  <dig> .
abbreviations
acgh: array comparative genomic hybridization; aic: akaike’s information criterion.
ascat: allele-specific copy number analysis of tumors; bac: bacterial artificial chromosome; baf: b-allele frequency; bic: schwarz’s bayesian information criterion; cbs: circular binary segmentation; cnv: copy number variation; dtc: disseminated tumor cells; fl: fused lasso; hts: high-throughput sequencing; iqr: interquartile range; mad: median absolute deviation; mlpa: multiplex ligation-dependent probe amplification; pcf: piecewise constant fitting ; roc: receiver operating characteristic curve; snp: single-nucleotide polymorphism.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
the study was initiated by kl, albd and ocl.
gn, kl and ocl drafted the manuscript.
the software was written by gn with contributions from kl based on algorithms developed by gn, kl and ocl.
pvl, hkmv, mbe and lob contributed with examples and in discussions of the manuscript and software.
omr, sfc, rr and cc provided and analysed the mlpa data.
all authors have read, commented on and accepted the final manuscript.
supplementary material
additional file 1
this pdf-file contains a formal description of the iterative pcf-based winsorization algorithm.
click here for file
additional file 2
this pdf-file contains a description of the three data sets used in this paper.
click here for file
additional file 3
this pdf-file describes a comparison of segmentations performed by cbs and pcf on a micma sample.
click here for file
acknowledgements
gn, kl and ocl received funding from the centre of cancer biomedicine  at the university of oslo for equipment and travelling.
pvl is a postdoctoral researcher of the research foundation - flanders .

 in cancer, the path from normal to malignant cell involves multiple genomic alterations including losses and gains of genomic dna. a long series of studies have demonstrated the biological and clinical relevance of studying such genomic alterations. genome-wide scans of copy number alterations may be obtained with array-based comparative genomic hybridization, snp arrays and high-throughput sequencing. after proper normalization and transformation of the raw signal intensities obtained from such technologies, the next step is usually to perform segmentation to identify regions of constant copy number. many segmentation algorithms are designed to analyse samples individually, while most studies involve multiple samples, multiple tracks, or both. joint handling of multiple samples is computationally and conceptually challenging, see e.g.. most systematic approaches for this problem are based on individual segmentation of each sample followed by post-processing to combine results across samples, while some recent publications propose strategies for joint segmentation of all samples. recently, the emergence of new technologies have pushed the limit of genomic resolution, opening new vistas for studying very short aberrations, including aberrations affecting only part of a gene or gene regulatory sites in the dna. a major challenge raised by these novel technologies is the steadily growing length of the data tracks, which drastically increases the demand for computationally efficient algorithms. the occurrence of extreme observations of biological or technical origin pose an additional challenge, as most segmentation methods are substantially affected by such observations. picard et al. propose a least squares based segmentation method that results in a piecewise constant fit to the copy number data. their approach assumes that the user either supplies the desired number of segments or leaves to the method to automatically determine this number. in this paper, we describe a related approach. in particular, the proposed method utilizes penalized least squares regression to determine a piecewise constant fit to the data. introducing a fixed penalty γ> for any difference in the fitted values of two neighboring observations induces an optimal solution of particular relevance to copy number data: a piecewise constant curve fully determined by the breakpoints and the average copy number values on each segment. the user defined penalty γessentially controls the level of empirical evidence required to introduce a breakpoint. given the number of breakpoints, the solution will be optimal in terms of least squares error. to achieve high processing efficiency, dynamic programming is used. to further increase computational efficiency, a novel vector based algorithm is proposed, and even further speed optimization is obtained through heuristics. a central aim of the present work has been to provide methodology and high-performance algorithms for solving single- and multiple-track problems within a statistically and computationally unified framework. all proposed algorithms are embedded in a comprehensive software suite for copy number segmentation and visualization, available as the bioconductor package copynumber. main features of the package include: • independent as well as joint segmentation of multiple samples • segmentation of allele-specific snp array data • preprocessing tools for outlier detection and handling, and missing value imputation. • visualization tools implementation systems overview the copynumber package provides functionality for many of the tasks typically encountered in copy number analysis: data preprocessing tools, segmentation methods for various analysis scenarios, and visualization tools. figure shows an overview of the typical work flow. input is normalized and log2-transformed copy number measurements from one or more acgh, snp-array or hts experiments. allele-frequencies may also be specified for the segmentation of snp-array data. it is strongly recommended to detect and appropriately modify extreme observations prior to segmentation, as these can have a substantial negative effect on the analysis. for this purpose, a specially designed winsorization method is included in the software package. a missing-value imputation method appropriate for copy number data is also available. segmentation methods for three different scenarios are implemented in the package. all these methods are referred to as piecewise constant fitting algorithms and seek to minimize a penalized least squares criterion. in single sample pcf, individual segmentation curves are fitted to each sample. in multi-sample pcf, segmentation curves with common segment borders are simultaneously fitted to all samples. in allele-specific pcf, the segmentation curves are fitted to bivariate snp-array data, providing identical segment borders for both data tracks. a set of graphical tools are also available in the package to visualize data and segmentation results, and to plot aberration frequencies and heatmaps. also included are diagnostics to explore different trade-offs between goodness-of-fit and parsimony in terms of the number of segments. in the remaining part of this section, a formal description of the algorithms is given. however, note that these details are not a prerequisite for reading later sections or for using the copynumber package. preprocessing: outlier handling a challenging factor in copy number analysis is the frequent occurrence of outliers - single probe values that differ markedly from their neighbors. such extreme observations can be due to the presence of very short segments of dna with deviant copy numbers, to technical aberrations, or a combination. when identification of cnvs is a purpose of the study, the multi-sample method described below may be applied for such detection. however, when the focus is on detection of broader aberrations, the potentially harmful effect of extreme observations on aberration detection methods induces a need for outlier handling procedures. since the copynumber package is based on least squares, an extreme observation will tend to cause the detection of a short segment. when searching for broader segments, such short segments will represent noise and may also affect the identification of other segments. we therefore now describe a procedure for reducing the effect of extreme observations, while the effects of this method will be considered in the results and discussion section. winsorization is a simple transformation reducing the influence of outliers by moving observations outside a certain fractile in the distribution to that fractile. for identically distributed observations y …,yp, the corresponding winsorized observations are defined as yjw=Ψ where Ψ=Ψ=−θ,y<−θθ,y>θy,otherwise. here, θ> determines how extreme an observation must be to be relocated, as well as the replacement value. a common choice is θ=τs, where typically τ∈ and s is a robust estimate of the standard deviation. a robust scale estimator is the median absolute deviation, defined as the median of the values |yj−m^|, where m^ is the median of y …,yp. for normally distributed observations, sm= ·mad corresponds to sd. winsorization of copy number data may be achieved by first estimating the trend in the data and then winsorizing the residuals. let the observations representing copy numbers in p genomic loci be y=, ordered according to genomic position. a simple estimator of the trend is the median filter. the trend estimate m^j in the jth locus is then given by the median of yj−k,…,yj + k for some k> e.g. k= the sd of the residuals yj−m^j may then be estimated with the mad estimator sm, and winsorized observations y1w,…,ypw obtained by yjw=m^j+Ψ. often, such simple and fast winsorization is sufficient. however, copynumber also includes an iterative procedure with improved trend estimation based on the segmentation procedures described below. single sample segmentation consider first the basic problem of obtaining individual segmentations for each of a number of samples. suppose attention is restricted to one chromosome arm on one sample. for each of the p loci, the obtained measurement can be conceived of as a sum of two contributions: yj=zj+εj where zjis an unknown parameter reflecting the actual amount of sample dna at the j’th locus and εjrepresents measurement noise. a breakpoint is said to occur between probe j and j + if zj≠zj + the sequence z …,zp thus implies a segmentation s={i …,im} of the chromosome arm, where i consists of the probes before the first breakpoint, i consists of the subsequent probes until the second breakpoint, and so on. to fit model, we minimize the penalized least squares criterion ∑j=1p2+γ·|s| with respect to the sequence z …,zp. here, |s| denotes the number of segments in s, and γ> is a constant that controls the trade-off between seeking a good fit to the data and restraining the number of level shifts. the minimizer z^ …,z^p of is fully determined by the segmentation s, since the best fit z^j on a given segment i is the average y¯i of the observations on that segment. substituting the latter into we obtain the equivalent criterion: l=∑i∈s∑j∈i2+γ·|s| =∑i∈s∑j∈iyj2−∑i∈s2/ni+γ·|s| where ni denotes the number of probes in segment i. note that the first term in does not depend on the segmentation s, hence minimization of is equivalent to minimizing l′=−∑i∈s2/ni+γ·|s|. naive optimization of the cost function with respect to the segmentation s requires examination of every possible division of the probes on a chromosome arm into segments. for large p, this is not practically feasible. however, a much more efficient implementation based on dynamic programming and requiring only o operations is available. dynamic programming is a method for solving complex problems by breaking them down into simpler subproblems, and specifically for problems where global decisions can be decomposed into a series of nested smaller decision problems. the crucial observation that allows the use of dynamic programming to solve the present segmentation problem is that the optimal segmentations on each side of a breakpoint are mutually independent. this can be used to iteratively build up a solution to the global segmentation problem. suppose we know the optimal segmentations from the first probe up until the st probe. assume furthermore that the optimal segmentation for the k first probes contains breakpoints. then the optimal segmentations from the last of these breakpoints and downwards has already been computed. thus, by solving the above subproblems iteratively for increasing k, each step can utilize the results from the previous steps. more formally, assume that the optimal segmentation of 1…r and the corresponding total error er are known for all probes r<k. to extend the solution to r=k, first note that there must be a last segment starting at some index j≤k. from we find that the cost term associated with that segment is: djk=1j−k−1∑r=jkyr then the total error for the optimal solution up until index k is found by minimizing the cost over the possible start positions j of the last segment. this cost consists of three terms: the cost of the last segment, the optimal cost of the segmentation up until that point and the penalty for the break point: ek=minj∈{ …,k} where e0= the main work load of the above computation is to determine djk for all 1≤j≤k≤p. in interpreted languages where loop execution is often quite inefficient, a considerable improvement of performance may be obtained by utilizing native-language vector operations. let ajk=∑r=jkyr, ak= and dk=. then we may calculate all required coefficients through a simple recursion: ak=+ykdk=−ak∗ak/ where = and operators are vector-based. hence, addition of a vector and a scalar adds the latter to each component of the former, and multiplications and divisions are performed component-wise on the operands, e.g., ak∗ak=. algorithm summarizes the computations. algorithm 1: single sample pcf input: log-transformed copy numbers y …,yp; penalty γ> output: segment start indices s …,sm and segment averages y¯ …,y¯m. calculate scores by letting a0= and e0= and iterate for k=1…p: • ak= + yk • dk=−ak∗ak/ •ek= storing also the index tk∈{,…,k} at which the minimum in the last step is achieved. find segment start indices s1=tp,s2=ts1−1…,sm= where m≥ find segment averages y¯m=ave−1) for m= …,m, where s0=p + throughout the paper we will tacitly assume that the penalty for the ith sample is γi=γσ^i where σ^i is the estimated sample specific residual variance. in this way, we avoid scale dependency, and obtain consistent results for samples with equal signal-to-noise ratios. such rescaling is also done by default in copynumber. note that replacing the data yji for the ith sample with yji/σ^i for j= …,p, and rescaling after estimation, has the same effect. in copynumber, the algorithm has also been extended to allow a constraint on the least number of probes in a segment. multi-sample segmentation detection of very short or very low amplitude segments requires a small penalty γ, with low specificity as a potential result. however, when such segments are common to several samples, joint segmentation of multiple samples is an additional mechanism to increase sensitivity. this is a main motivation for introducing multi-sample segmentation methods that impose common breakpoints across all samples. such methods are potentially useful for discovery of copy number variations and in those instances where the origin of the samples implies that segment boundaries are partly shared. multi-sample segmentation with high penalty on breakpoints may also be used to obtain low-dimensional descriptions of the data, which may form the basis for defining variables to be used in statistical procedures relating aberration patterns to clinical outcome. in the following, we describe a direct generalization of single sample pcf to handle multiple samples simultaneously, obtaining common breakpoints for all the samples with minimal residual sum of squares for a given number of breakpoints. suppose copy number measurements yi= for samples i=,…,nare obtained at the same loci in each sample. by direct generalization of the criterion, we seek in multi-sample pcf the minimizer of l=∑i=1nl where l is defined as in and s is a given segmentation common to all samples. algorithm 2: multi-sample pcf input: log-transformed copy numbers for n samples y …,yp∈rn; penalty γ> output: common segment start indices s …,sm and segment averages y¯ …,y¯m∈rn. calculate scores by letting a0= and e0= and iterate for k=1…p: • ak= + yk • dk=−1t/ •ek= storing also the index tk∈{,…,k} at which the minimum in the last step is achieved. find segment start indices s1=tp,s2=ts1−1…,sm= where m≥ find segment averages y¯m=ave−1) for m= …,m, where s0=p + the multi-sample pcf algorithm is in principle quite similar to single sample pcf. however, when updating the solution from k− to k, the sums and sums of squares for the segments must be accumulated and stored separately for each sample. this can still be done iteratively, implying that the computational effort will be approximately equal to carrying out single sample pcf on the same set of samples. since the noise level may vary between samples, normalisation of the samples prior to segmentation and corresponding rescaling after estimation is advisable. it may also be desirable to scale the samples, e.g. to adjust for different tumor percentages. thus, prior to running multi-sample pcf, we may replace yi by wiyi/σ^i for i= …,n, where wi are weights and σ^i is an estimate of the sd. in copynumber normalization is performed by default for multi-sample pcf while further weighting is left as an option for the user. allele-specific segmentation the pcf algorithm is easily adapted to variants of the basic segmentation problem discussed above. here, we consider an adaptation to handle snp genotype data. we then have for each snp locus a measurement of copy number as well as the b allele frequency. we may also have measurements of copy number only for a number of additional loci. the b allele frequency is a number between and indicating the allelic imbalance of a snp. for a homozygous locus we have baf close to or while for a heterozygous locus with an equal number of the two alleles a and b, baf will be close to. an imbalance between the number of a’s and b’s results in a baf value deviating from. a change in the total number of copies of a segment will alter the logr value, hence result in a level shift in the logr track. unless the copy number change is balanced with respect to the two alleles, the baf value will also change. in cases involving multiple copy number events at the same locus, the change may manifest itself only in one of the two tracks. for example, a loss of one copy of a followed by a gain of one copy of b would lead to unchanged logr and changed baf. the purpose of the allele-specific pcf algorithm is to detect breakpoints for all such events. it fits piecewise constant curves simultaneously to the logr and the baf data, forcing breakpoints to occur at the same positions in both. we emphasize that the purpose of the allele-specific pcf algorithm is segmentation only and not to make allele-specific copy number calls. however, such calls can be made on the basis of the segmentation described below, and this is done e.g. in the ascat algorithm which estimates allele-specific copy numbers as well as the percentage of cells with aberrant dna and the tumor ploidy. suppose the data are given by for j= …,p, where rj denotes the logr value and bj the baf value at the jth locus. for copy number probes, only rj is given and bj will be missing. for germline homozygous probes, the baf values are noninformative and should be omitted from the analysis. if the germline genotype is known, the user should replace the corresponding baf values by na. if the genotype is not known, the algorithm will apply a proxy to handle this issue. prior to segmentation, the allele-specific pcf algorithm performs the following steps: the baf data are mirrored around by replacing bjwith 1−bjif bj>. baf values bj<θare replaced by na. by default θ=. if germline homozygous probes have previously been replaced by na’s, let θ= let b~ …,b~m denote the nonmissing b allele frequencies. corresponding copy number values r~ …,r~m are found by pairing each logr probe with the nearest b-allele probe and then averaging logr values paired to the same b-allele probe. finally, let y1= and y2=. the remaining part of the allele-specific pcf algorithm is then essentially an adaptation of the multi-sample pcf algorithm applied to two samples. it finds a common segmentation s for the two tracks by minimizing the penalized criterion l=l+l where l is defined as in. fast implementations of pcf the pcf algorithms may be generalized to allow breakpoints only at certain prespecified positions. combined with simple heuristics, this may be used to further enhance the computational speed of pcf. for brevity we describe only the single sample segmentation case here; however the copynumber package contains fast implementations of both single- and multi-sample pcf. computationally inexpensive methods can be used to identify a set of potential breakpoints among which the breakpoints of the solution to are highly likely to be found. suppose we restrict our attention to such a set of potential breakpoints. all relevant information for solving the optimization problem in may then be condensed into three arrays containing the number of observations between two potential breakpoints, the corresponding sum of the observations and the sum of squares. based on these quantities, pcf may be used with straightforward modifications. since the algorithm is of order o, where q is the number of potential breakpoints, the potential increase in speed is substantial. algorithm outlines the procedure, while possible heuristics for finding potential breakpoints are discussed below. one way to identify potential breakpoints is to use high-pass filters, i.e. a filter obtaining high absolute values when passing over a breakpoint. the simplest such filter uses for each position i the difference ∑j=i+1i+kyj−∑j=i−k+1iyj for some k. to reduce artifacts due to the abrupt edges of such a filter, the copynumber implementation assigns half weight to the outer 1/ of the observations on each side. fast implementations of such filters in r may be obtained using the cumsum function. we currently use two filters with k= and respectively; additionally the single sample pcf implementation includes a filter searching for aberrations of length equal to the lowest accepted one. these filters together identify about 15% of the probe positions as potential breakpoints. an additional way to speed up the computations on long sequences is to initially divide the sequence into overlapping subsequences, and iteratively find the solution. having found the solution for the m first subsequences, we use high-pass filters to detect potential breakpoints for subsequence m + and then use the fast pcf algorithm with the latter potential breakpoints as well as those found by pcf on earlier subsequences. the intention behind this iterative approach is to reduce potential boundary effects. due to the quadratic order of the algorithm, this division into subsequences implies a substantial efficiency gain. in copynumber, subsequences are used when the chromosomal arm length exceeds probes, with subsequences of length and overlap algorithm 3: fast pcf input: log-transformed copy numbers y …,yp; penalty γ> output: segment start indices s …,sm and segment averages y¯ …,y¯m. apply heuristics to find potential breakpoints r r …,rq, where r0= and rq=p + form aggregates by letting uk=∑j=rk−1rk−1yj, where k= …,q. calculate scores by letting a0=, c0=, e0= and iterate for k= …,q: • ak= + uk • ck= + rk−rk−1 • dk=−ak∗ak/ck • ek= storing also the index tk∈{,…,k} at which the minimum in the last step is achieved. find segment start indices s1=rtq,s2=rts1−1…,sm= where m≥ find segment averages y¯m=ave−1) for m= …,m, where s0=p +  selection of penalty the selection of parameters determining the trade-off between high sensitivity and high specificity is important in all segmentation procedures. in pcf, this is controlled by the single penalty parameter γ. a number of general model selection criteria exist, such as cross-validation, the akaike information criterion and the related schwarz’s bayesian information criterion. however, model selection for copy number segmentation is complicated by several factors. first, the distribution of the data at hand may vary substantially. an important example is the presence of local trends mimicking smaller aberrations; such low-amplitude “waves” in the data may e.g. be due to variations in gc-content. second, the purpose of the analysis may favor either higher sensitivity or higher specificity. for example, in clinical studies aimed at finding prognostic markers, the main focus may be on the most pronounced and commonly occurring deviations, while detecting more sporadic aberrations may simply increase the noise level. in our experience, the above model selection criteria tend to give too small penalty estimates and thus undersmooth the data. this is consistent with previous investigations showing that aic and bic are not appropriate for the breakpoint problem. simulation studies of specificity may suggest a lower bound on the penalty γ. for this purpose, sequences of independent and normally distributed observations without underlying aberrations were generated, and pcf was applied with different choices of γ. at γ= the number of falsely called aberrations is about per probes, at γ= roughly per probes, at γ= roughly per probes, and for γ≤ the number of falsely called aberrations is substantial. this suggests γ≈8− as a lower bound. since the number of false aberrations per chromosome increases with increasing probe density, low values are most relevant for arrays with low probe density. in the presence of local trends, the number of false calls tends to inflate and the penalty should thus be increased above the lower bound. a fairly conservative penalty of γ= is the default in the copynumber package. this provides a starting point for exploration of the best penalty value for the specific problem at hand, however a systematic inspection of results obtained for different penalties is advisable. figure illustrates the effect of changing γ. notice that the main features in the data are captured across the whole range of γ-values, while finer details are only evident for smaller values. aberration calling aberration calling is used for detection of recurring alterations and in many other analyses. introducing a parameter θ> that determines the sensitivity of the aberration calling, we call probes for which z^j<−θ as losses and probes for which z^j>θ as gains. optionally, different thresholds θ + and θ− may be used for gains and losses. to examine how well pcf aberration calling manages to distinguish between normal and aberrant regions, performance was compared with a very accurate measurement method. specifically, aberration calls obtained with pcf on the basis of m snp array data on samples were compared with calls obtained with mlpa. since mlpa is limited to a small set of genomic positions, only loci were used for the comparison. in all samples combined, mlpa identified aberrant and normal loci. using the mlpa-classification as the gold standard, the sensitivity and specificity of pcf aberration calling were calculated for a range of threshold values θ. figure shows the resulting roc curves, and panel illustrates how the results for pcf depend on the choice of γ. importantly, aberration calling appears to be only moderately dependent on the choice of parameter values over a fairly wide range of γ-values. single- versus multi-sample segmentation whether the initial segmentation of a dataset is most appropriately done using single- or multi-sample methods depends both on the purpose and the data. using methods with common breakpoints for samples will increase the power for detecting concordant but quantitatively weak segments, while it will reduce the ability of detecting discrepant breakpoints. a well known example of aberrations with common boundaries is germline copy number variants, thus some proposed algorithms for cnv detection utilize segmentation with joint segment borders. another important example of samples with common segment boundaries arises when the samples originate from different clones of the same tumor. this is illustrated below in two examples, one on disseminated tumor cells from breast carcinomas, the other on tumor clones found at successive biopsies from lymphoma patients. recent reports on marked variations in aberration patterns within the same tumor is likely to increase the number of studies using several samples taken from each tumor. what is common as well as what differs in the aberration patterns will then be of interest, motivating the combined use of single- and multi-sample methods. in applications searching for genomic copy number hot spots with relevance to cancer development, it may be important to utilize the precise delineation of the aberrations found in each sample, and thus the use of single-sample methods is most appropriate. the identification of the relevant recurrent aberrations may then utilize post processing tools like gistic, kcsmart or cghmcr. if focus is on clustering samples or on constructing regression variables for relating more broad aberrations to clinical outcome, one may consider using multi-sample methods. however, to be useful, the estimates from the multi-sample methods should in a proper way reflect the main information content in each sample. this implies that a multi-sample analysis should result in estimates approximating those obtained from single-sample analyses. figure shows heatmaps of results from single- and multi-sample pcf for breast cancers from the so-called micma data set analyzed on 244k agilent arrays. the main features appear to be well reflected in the multi-sample analysis. on a more detailed level, differences can be observed: the multi-sample solution misses some short aberrations occurring in only a few samples, aberration borders are sometimes slightly shifted, and longer segments obtained with single sample pcf are often divided into subsegments with slightly different copy number estimates. the moderate difference between the results of single- and multi-sample pcf was also confirmed by a comparison of the ability to detect specific aberrations as revealed by comparison to mlpa analyses, see figure 3b. this indicates that at least for cancer types where aberrations are focused in certain areas of the genome, methods using joint boundaries might be considered for constructing variables to be used in further statistical analysis. comparing tracks: analysis of disseminated tumor cells disseminated tumor cells are detected in the bone marrow of some patients with breast carcinomas. the presence of dtcs in the bone marrow identifies patients with less favorable outcome, and genomic characterization of such cells is of substantial interest. it is still an open question to what extent the aberration patterns in dtcs correspond to those found in the primary tumor; the dtcs may potentially have obtained new aberrations or, alternatively, the cells may have originated from subclones of the tumor with less aberrations. it is possible to analyze single cells using acgh; however, currently the noise level is high, making it difficult to draw definitive conclusions from a single cell. however, since segment boundaries are assumed to be partly common, we tested the multi-sample pcf algorithm on breast cancer cases from which dtcs were available. figure shows the results on a set of dtcs and the corresponding primary tumor from one such patient. since multi-sample pcf is used, segment boundaries are common, while the estimated level in each segment is determined by the individual dtc/primary tumor. in figure two of the single cells seem to have a pattern similar to the primary tumor. the last one has an essentially flat profile and is likely to be a hematopoietic cell misclassified as a tumor cell. these data thus indicate that the aberration pattern of the dtcs quite closely reflect that of the primary tumor. with only two single cells present, figure primarily shows that dtcs inherit the aberrations of the primary tumor; with higher numbers of cells, multi-sample pcf may also be used to search for aberrations found in dtcs but not in the tumor. defining variables: genetic evolution in follicular lymphoma follicular lymphoma is normally a slowly progressing malignancy, but relapses are common and the disease is usually fatal. in a recent study, biopsies from patients diagnosed with follicular lymphoma were evaluated using a custom-made acgh platform consisting of 3k bac/pac probes. a whole-genome view of aberration frequencies and highly correlated aberrations are shown in figure although the delineation of segments varied between biopsies, several areas with a high frequency of aberrations could be detected. to try to identify aberrations with prognostic potential, we therefore found a common segmentation for the initial biopsies taken from each of the patients using the multi-sample pcf algorithm. removing very low variance segments, segments remained. the corresponding copy number estimates were used as covariates in a multivariate cox proportional hazards regression. this revealed segments for which gains were significantly associated with a survival disadvantage. a particularly strong association was detected for gains on chromosome x in male patients. to study the relation between successive biopsies taken from the same patient, multi-sample pcf was applied to each patient individually. as expected, many aberrations are common, but interestingly, some aberrations are present in early biopsies and not in later ones. this contradicts the hypothesis of linear development which states that late tumor clones arise from earlier ones, and supports the alternative hypothesis of parallel evolution in different lymph nodes. allele-specific copy number analysis in breast cancer copy number alterations have been extensively studied in breast cancer. to what degree gains and losses are associated only with certain alleles has been less studied. in a recent study, genotyping of breast carcinoma samples was performed using illumina 109k snp arrays, and the ascat method was used to infer the allele-specific copy numbers at each locus. however, to do this we first had to segment the data; for this purpose we applied allele-specific pcf segmentation to all samples. in figure the result of this segmentation is shown for one particular sample and two different chromosomes. in figure 7a, the segmentation of chromosome is shown, and we clearly identify three segments on the p-arm with copy numbers less than two, larger than two, and identical to two. suppose we consider only germline heterozygous loci, in which case the allelic ratio is 1/ when no aberrations are present. the baf track reveals allelic imbalance in the first two segments, and more pronounced in the first segment than in the second. this is consistent with a loss of one copy in the first segment, and a gain of one copy in the second segment. the third segment has an allelic ratio of 1/ notice that in case of allelic imbalance, the observed allele ratio is substantially closer to than expected by the above theoretical ratios. this attenuation of the signal is due to technical issues like cross-hybridization, as well as the fact that in reality the tumor is a mixture of cells with normal dna and tumor cells with aberrant dna. in figure 7b we notice a sharp trough in the logr track on 17p, accompanied by an allelic ratio close to. if for a certain snp locus one allele is substantially more frequently gained than the other allele, one may hypothesize that the former allele is subject to a larger selective pressure to change copy number. this, in turn, may be an indication of different roles being played by the two alleles with respect to cancer progression and evolution, suggesting that loci subject to allelic skewness can be potential unique markers for breast cancer development. even from a relatively small number of samples, probes with highly significant allelic skewness have been identified in a genome-wide statistical evaluation. outliers and winsorization while least squares methods are often favored due to their optimality properties, they are also known to be sensitive to extreme observations. thus, except if the purpose is to search for short aberrations of biological origins, we advise the use of an outlier handling procedure. to evaluate the proposed winsorization scheme, we first established a suitable way of simulating extreme observations. a classical way is to use “contaminated normals”, where the error distribution is a mixture of two normal distributions. with probability 1−α the error is drawn from a distribution n, and with probability αfrom n, typically with d= and α=. we compared the fraction of outliers in observed copy number data to the corresponding fractions in normals and contaminated normals, using mad to estimate sds. for the normal distribution, the fraction of observations outside sd is % and outside sd < %, while these fractions for the 5% contaminated normal are % and %. for the agilent 244k used on the micma dataset, the fractions were % and %, that is, slightly above the values for the contaminated normal. for the 44k agilent and illumina 109k applied to the same data, the percentages were slightly lower, however still indicating that the 5% contaminated normal is an appropriate distribution when evaluating robustness of copy number assessment procedures. inspection of data obtained by the 318k illumina, x 180k agilent and nimblegen m arrays also confirmed the existence of substantial amounts of outliers. the pcf algorithm was tested with and without winsorization on simulated data with outliers. outliers in the contaminated distributions may cause the detection of short false aberrations; such spikes occurred roughly ten times per probes, as compared to less than two times per for uncontaminated data. table further shows that winsorization efficiently reduces the number of falsely detected aberrations and make results for the contaminated distribution roughly equal to the ones for the normal distribution. in line with these observations, outliers tend to change the form of aberrations, while winsorization brings the distribution fairly close to the one found for normal data. shown is the effect of winsorization on simulated data with outliers and artificial aberrations. two types of aberrations are considered: aberrations of height and length probes and aberrations of height and length the contamination consists of normals with sd= and the mad estimate of sd equals. sensitivity is the percentage of amplified probes that are detected as amplified, while specificity is the percentage of non-amplified probes classified as such. the false aberration column gives the percentage of aberrations not covering the central part of the real amplifications. another way to avoid that a few extreme observations result in a segment is to impose a lower limit on the length of a segment. with a lower length limit of five probes, we found about twice as many false spikes as with winsorization when adjusting γto give equal sensitivity for true aberrations. still, simulations indicate that a lower limit on segment length is valuable in combination with winsorization. note that outliers of biological origin will be more extreme if the technology has an inherent low noise level, as is the case, e.g., for bac arrays and for high throughput sequencing. thus, outliers are not a sign of inappropriate functioning of a technique, but a characteristic of the data requiring consideration in the analysis. in summary, copy number data tend to contain a high fraction of outliers. these outliers often induce false aberrations, but simple procedures like winsorization will efficiently reduce these undesired effects. computational performance in r, using the vector based pcf implementation described in algorithm implies a substantial efficiency gain over loop based implementation, roughly a 10- times reduction in time requirements. the fast implementation of pcf gives a further marked reduction in computing time. on the micma 244k dataset, the implemented fast version is about times faster than the exact one, and uses around minutes to process the samples. the multi-sample method was slightly faster than the single sample version. the average computation time per sample is shown for copynumber, dnacopy v. and cghflasso v - on the micma k data set and on the logr values from an illumina m snp array data set. the iqr over samples is given in parenthesis. the methods were applied to both raw data and data with outliers removed using the winsorization method. all tests were performed on a pc with a ghz intel i cpu with gb of memory running windows and r.. the deviations between the solutions found by the exact pcf and fast pcf on the micma set were small; in terms of reduction in variance below %. the differences observed for the curves were typically small shifts in the border of aberrations. thus, we conclude that the results from the fast procedure for practical purposes may be regarded as global solutions to, and the fast version is therefore used by default in copynumber. we also compared the performance of pcf with two other segmentation methods: circular binary segmentation and fused lasso regression. in comparison studies, cbs has shown good performance in terms of sensitivity and false discovery rate. it is probably the most commonly used freely available algorithm and is also implemented in several commercial analysis tools. cbs is available in the r package dnacopy, which is used for this comparison. fl is a more recent proposal implemented in the r package cghflasso, and is one of three preferred methods in the web-based segmentation tool cghweb. using default parameter settings, we compared the computing times of pcf, cbs and fl on the samples in the k micma data set, and on samples from a m illumina snp array. table gives the average computation time per sample. with no preprocessing of the data, pcf is on average 3- times faster than cbs on both data sets, and about times faster than fl on the largest data set. note that copynumber detects and operates on chromosome arms, while dnacopy operates on whole chromosomes. this partly explains the difference in performance between pcf and cbs for the micma data set; for the illumina data this has little impact due to the iterative approach used in pcf for the longest sequences. pcf was also markedly faster in evaluations based on simulated data; however, comparisons are complicated by the fact that the speed of cbs depends on the data in a nontrivial manner. as seen from the iqrs listed in parentheses in table the speed of cbs is quite variable from sample to sample while pcf and fl is nearly constant. moreover, the table shows that cbs runs 2- times slower when outliers have been removed using winsorization, underlining that the performance of cbs is highly data dependent. we underline that the above-mentioned results only relate to the current r implementations. as mentioned in the introduction, pcf is conceptually similar to the cghseg method described by picard et al., and we also examined the computational performance of this method using the implementation in the r package cghseg. using the version of cghseg that requires a prespecified number of segments for each chromosome, the algorithm is fast, although the speed depends on the number of segments. using the full cghseg algorithm that automatically determines the number of segments, the algorithm is very slow for high-resolution data. hence, making a fair comparison between pcf and cghseg is difficult. segmentation accuracy we further compare the accuracy of the segmentation solutions found by pcf and cbs. figure 3c shows roc curves using mlpa classifications as the truth, and then applying a range of aberration calling thresholds to pcf estimates, cbs estimates, a running median with window size and raw copy number data. results for pcf and cbs are similar, both achieving high sensitivity and specificity. the running median also gives good results, illustrating that many probes are fairly easy to classify and that the gain obtained by using methods like cbs and pcf is mainly an improved classification close to borders between segments. we also repeated the simulation study in where cbs was found to be the most sensitive method while also having the lowest false discovery rate. again, we found that pcf and cbs had very similar performance. a more detailed comparison of segmentation results shows that overall results are quite similar for single sample pcf and cbs, however for both methods the results depend on the choice of parameter values and the handling of extreme observations, see additional file in conclusion, pcf and cbs typically provide similar results and have equivalent accuracy when parameters are tuned appropriately.  copy number segmentation based on least squares principles and combined with a suitable penalization scheme is appealing, since the solution will be optimal in a least squares sense for a given number of breakpoints. we have proposed a suite of platform independent algorithms based on this principle for independent as well as joint segmentation of copy number data. the algorithms perform similarly as other leading segmentation methods in terms of sensitivity and specificity. furthermore, the proposed algorithms are easy to generalize and are computationally very efficient also on high-resolution data. the bioconductor package copynumber offers a user-friendly interface to the proposed algorithms. several extensions and modifications of the proposed least-squares framework are possible. in principle, the l2-based distance measure used in the current implementation of pcf is easily extended to general lp-distances. however the current implementation is highly optimized for l and other distance measures would require substantial heuristics to obtain comparable computational performance. another extension is to introduce locus specific penalties for breakpoints, thus essentially introducing a prior on the location of breakpoints. work in progress includes specialized routines to handle high throughput sequencing data more efficiently and joint analysis of multiple samples in allele-specific pcf. availability and requirements project name: copynumber project home page:http://heim.ifi.uio.no/bioinf/projects/copynumber/ operating system: all systems supporting the r environment programming language: r other requirements: no license: gnu artistic license. abbreviations acgh: array comparative genomic hybridization; aic: akaike’s information criterion. ascat: allele-specific copy number analysis of tumors; bac: bacterial artificial chromosome; baf: b-allele frequency; bic: schwarz’s bayesian information criterion; cbs: circular binary segmentation; cnv: copy number variation; dtc: disseminated tumor cells; fl: fused lasso; hts: high-throughput sequencing; iqr: interquartile range; mad: median absolute deviation; mlpa: multiplex ligation-dependent probe amplification; pcf: piecewise constant fitting; roc: receiver operating characteristic curve; snp: single-nucleotide polymorphism. competing interests the authors declare that they have no competing interests. authors’ contributions the study was initiated by kl, albd and ocl. gn, kl and ocl drafted the manuscript. the software was written by gn with contributions from kl based on algorithms developed by gn, kl and ocl. pvl, hkmv, mbe and lob contributed with examples and in discussions of the manuscript and software. omr, sfc, rr and cc provided and analysed the mlpa data. all authors have read, commented on and accepted the final manuscript. supplementary material additional file 1 this pdf-file contains a formal description of the iterative pcf-based winsorization algorithm. click here for file additional file 2 this pdf-file contains a description of the three data sets used in this paper. click here for file additional file 3 this pdf-file describes a comparison of segmentations performed by cbs and pcf on a micma sample. click here for file acknowledgements gn, kl and ocl received funding from the centre of cancer biomedicine at the university of oslo for equipment and travelling. pvl is a postdoctoral researcher of the research foundation - flanders.
a comprehensive bioconductor package for copy number analysis is presented. the growing variety of tools for measuring genomic copy numbers, including various types of array-cgh, snp arrays and high-throughput sequencing, calls for a coherent framework offering unified and consistent handling of single- and multi-track segmentation problems. in addition, there is a demand for highly computationally efficient segmentation algorithms, due to the emergence of very high density scans of copy number. the r package copynumber is a software suite for segmentation of single- and multi-track copy number data using algorithms based on coherent least squares principles. the package offers a unified framework for single sample, multi-sample and multi-track segmentation and is based on statistically sound penalized least squares principles. conditional on the number of breakpoints, the estimates are optimal in the least squares sense. cancer progression is associated with genomic instability and an accumulation of gains and losses of dna.
 cancer progression is associated with genomic instability and an accumulation of gains and losses of dna. the growing variety of tools for measuring genomic copy numbers, including various types of array-cgh, snp arrays and high-throughput sequencing, calls for a coherent framework offering unified and consistent handling of single- and multi-track segmentation problems. in addition, there is a demand for highly computationally efficient segmentation algorithms, due to the emergence of very high density scans of copy number.  a comprehensive bioconductor package for copy number analysis is presented. the package offers a unified framework for single sample, multi-sample and multi-track segmentation and is based on statistically sound penalized least squares principles. conditional on the number of breakpoints, the estimates are optimal in the least squares sense. a novel and computationally highly efficient algorithm is proposed that utilizes vector-based operations in r. three case studies are presented.  the r package copynumber is a software suite for segmentation of single- and multi-track copy number data using algorithms based on coherent least squares principles. copy numberacghsegmentationallele-specific segmentationpenalized regressionleast squaresbioconductor 