BACKGROUND
lesions at dna level represent the cause of cancer and of many congenital or hereditary disorders.
the change of the number of copies of dna in a genomic region is one of the most common aberrations.
in normal cells each genomic segment is present in two copies, but, for example, in tumor cells the genome can present regions of deletions , gains  or amplifications .
thus, in general, the dna copy number along the genome can be represented as a piecewise constant function.
with microarray technology it is possible to simultaneously measure the copy number along the genome at hundred thousands of positions .
however, raw copy number data are generally very noisy.
hence, it is important to define a method which allows one to estimate the number of regions with different copy number, the endpoints of these regions  and their copy number.
several methods have been developed to solve this issue.
many methods consider the log <dig> ratio of the data  and model it as a normal random variable, since they assume that the errors are normally distributed.
we can roughly subdivide all of these methods into two classes: those ones that estimate the copy numbers as a piecewise constant function and the others that estimate the copy numbers as a continuous curve.
the methods belonging to the latter group are called smoothing methods.
among the methods belonging to the first class, we can find the following.
the circular binary segmentation  approach is a recursive method in which the breakpoints are determined on the basis of a test of hypothesis, with null hypothesis that in the interval considered there is no change in copy number  <cit> .
picard et al.
<cit>  used a piecewise constant regression model, where the parameters are estimated maximizing a penalized likelihood .
this method is usually denoted with the abbreviation cghseg.
the glad method is another piecewise constant regression method, but in this case the parameters are estimated maximizing a weighted likelihood  <cit> .
fridlyand et al.
<cit>  applied hidden markov models , while marioni et al.
<cit>  defined an hmm method which takes into account the distance among the data points .
recently, nilsson et al  <cit>  derived a segmentation method based on total variation minimization, called rendersome.
it is optimized for gene expression data, but the authors affirm that it can be used also on copy number data.
among the smoothing methods, hsu et al.
<cit>  used a wavelet regression method with haar wavelet.
eilers and de menez  <cit>  applied a quantile smoothing regression , with the solution found by minimizing a loss function based on the l <dig> norm, to obtain a flatter curve.
huang et al.
<cit>  proposed smoothseg, i.e.
a smooth segmentation method based on a doubly heavy-tailed random-effect model.
we propose a piecewise constant regression method, using bayesian statistics, which appears appropriate when regions contain only few data points.
the original version of the method  was presented by hutter  <cit> .
in this paper we propose improved bayesian estimators of the parameters involved and we apply the model to dna copy number estimation.
finally, we compare our algorithm with some among the most cited or more recent methods, on artificial and real data.
our method was implemented in r and is freely available at  or in additional file  <dig>  furthermore, an r package will be soon available.
methods
in the first two subsections, we briefly describe the original bayesian piecewise constant regression method, explaining the hypothesis of the model and the estimation of its parameters with bayesian regression.
we emphasize the definitions of the original parameter estimators in order to show how we changed some of these estimators in the other subsections.
a brief explanation of the dynamic program for the computation of the estimators can be found in the additional file  <dig> .
regarding notations, we do not indicate explicitly the random variable to which a distribution is referred, if it is clear from the context.
for example, pk ≡ p or fy, m ≡ f.
hypotheses of the model
let y ∊ ℝn be a random vector, such that each component  is conditionally normally distributed:
yi|μ˜i <dig> σ2~n.
suppose also that y represents a noisy observation of a piecewise constant function, which consists of k <dig> horizontal segments.
then, the segment level at a generic position i  does not assume different values for each i, but the data are divided into k <dig> intervals  where μ˜tq−1+10=...=μ˜tq0=:μq <dig> for each q =  <dig> ...,k <dig>  hence, μq <dig> represents the level of the qth segment.
our goal is to estimate the levels μ0= of all the segments.
in order to do that, we first need to estimate the number of the segments k <dig> and the partition of the data t <dig>  from a bayesian point of view, μ <dig>  t <dig> and k <dig> are treated as random variables, hence we denote them with the corresponding upper case letters .
moreover, because of their randomness, we need to define a prior distribution for each of them to complete the model.
for the number of segments and the boundaries, we assume noninformative prior distributions:
p=1kmax⁡k∈k
p=1t∈tk,n,
where k = { <dig> ...,kmax} and tk,n is the subspace of ℕ0k+ <dig> such that t <dig> =  <dig>  tk = n and tq ∊ { <dig> ...,n - 1} for all q =  <dig> ...,k -  <dig>  in an ordered way and without repetitions.
about m, we assume that all its components are mutually independent and identically normally distributed,
m|ν,ρ <dig> k=k~n,
where ν ∊ ℝk, such that νq = ν for each q =  <dig> ...,k, and i ∊ ℝk × k, such that ip, q = δp,q for each p, q =  <dig> ...,k.
instead of these assumptions, we could assume a cauchy distribution for each yi or mq in order to model an observation whose noise has heavier tails, as previously done by hutter  <cit> .
original estimation: the bpcr method
the statistical procedure consists in a sequence of estimations due to the relationship among the parameters.
bpcr estimates the number of segments with the map  estimate given the sample point y,
k^:=arg⁡max⁡k∈kp,
and, given k^, also each boundary is estimated separately with its corresponding map estimate,
t^p:=arg⁡max⁡h∈{p,...,n−}p
for all p =  <dig> ...,k^ -  <dig>  finally, the rth moment of the level of the mth segment is estimated with its posterior mean.
since its computation needs the knowledge of the number of segments and the partition of the data, we replace them with the estimated ones,
μ^mr:=e
for all m =  <dig> ...,k^.
when r =  <dig> and we assume that y and m are normally distributed, the estimate turns out to be
μ^m=ρ2∑i=tm−1+1tmyi+σ2νρ2+σ <dig>
for all m =  <dig> ...,k^.
when the sample contains only one segment, the bayesian estimation of the posterior distribution of the levels should theoretically lead to a normal distribution, similar to a dirac delta function centered at ν^, since the levels can assume only one value from the data.
in fact, in this case, if we estimate ρ <dig> only using the data , then this value will be close to zero  and so the level will be estimated with ν^, the mean of the data ).
we can estimate the segment level m˜s at a generic position s, using the fact that it belongs to some segment m and in this segment it is equal to the corresponding mm.
then, summing over all the possible segments, we can compute its posterior distribution in the following way:
f=∑m=1k0∑i=0s−1∑j=snf
and the corresponding estimate of m˜s given k^ is
μ˜^s:=e,
for all s =  <dig> ...,n. the vector μ˜^ is called bayesian regression curve .
the probability distributions defined previously depend on the hyper-parameters ν, ρ <dig> and σ <dig> .
hutter  <cit>  suggested the following estimators:
ν^:=1n∑i=1nyi=y¯
ρ^2:=1n∑i=1n <dig>
σ^2:=12∑i=1n− <dig>
improved estimators of the number of segments
to understand the real meaning of the map estimator k^, we need to introduce the theory of the construction of a generic bayesian estimator.
in general, a bayesian estimator is defined in the following way.
let us suppose that z is a random variable whose distribution depends on an unknown parameter θ, which we want to estimate.
since we cannot exactly know the true value of the parameter, we consider it as a random variable Θ with a given prior probability distribution.
in order to measure the goodness of the estimation, we define an error  and we choose the estimator that minimizes the expected error given the sample z,
Θ^:=arg⁡min⁡θ′e.
the 0– <dig> error  is commonly used for a parameter which can assume only a discrete number of values.
the estimator corresponding to this error is the map estimator,
arg⁡min⁡θ′e=arg⁡max⁡θ′∑θδθ,θ′p=arg⁡max⁡θ′p.
obviously, if we use different types of errors, we can obtain different estimators.
in the following, we will use k^ to denote any estimator of k, while k^ <dig> to denote the parameter estimator k^ based on the 0– <dig> error.
using the 0– <dig> error, we give the same penalty to each value different from the true value, whether it is close to or far away from the true one.
to take into account the distance of the estimated value from the true one, we need to use other types of errors, which are based on different definitions of distance, such as,
absolute error := |θ - θ'|
squared error :=  <dig>
if the parameter θ ∊ ℝ, then the estimators corresponding to these errors are the median and the mean of its posterior distribution, respectively.
in our case, we denote these estimators of k <dig> with k^ <dig> and k^ <dig>
improved estimators of the boundaries
similarly to the previous subsection, we derive alternative boundary estimators, by considering different types of errors.
we denote the map boundary estimator defined in equation  with t^ <dig>
meaning of the estimator t^01
the estimator t^ <dig> is defined in such a way that each component minimizes the 0– <dig> error of the corresponding boundary, separately.
explicitly, given the sample point y and the segment number k <dig>  its estimate is
t^01=,...,arg⁡max⁡tk0−1∈tp,n),
where t = { <dig> ...,n - 1}.
t^ <dig> may be regarded as an approximation of the bayesian estimator that minimizes the error which counts the number of wrongly estimated boundaries:
sum 0-1 error=∑p=1k0−1=k0−1−∑p=1k0−1δtp <dig> tp,
that is
t^sum=arg⁡max⁡t∈tk <dig> n∑p=1k0−1p.
definition of the estimator Τ^joint
a problem of the latter estimator ) is its computational complexity, because it needs the computation of all the ordered combinations of the boundaries.
on the other hand, there are two reasons for which t^ <dig> is not a suitable estimator of the boundaries.
first, it does not take into account that the boundaries are dependent, because they have to be ordered, and second, in principle, it can have more than one component with the same value.
as a consequence, a theoretically more correct way to estimate the boundaries is minimizing the 0– <dig> error with respect to the joint boundary probability distribution .
then, given k <dig> and y, the boundary estimator becomes
t^joint=arg⁡max⁡t∈tk <dig> np.
definition of the estimators Τ^binerr and Τ^binerrak
we must notice that the estimators considered until now have the same length of the true vector of the boundaries.
in practice, the number of segments k <dig> is unknown, so that we should use k^.
as a consequence, if k^ is different from k <dig>  then, strictly speaking, we cannot minimize the previous types of error because the vectors have different length.
a way to solve this issue is to map each boundary vector into a vector τ∈ℝ0n+ <dig> in the following way:
t↦τ such that τi={1if ∃p such that tp=i0otherwise.
we denote with ttk,n the set of all the possible τ with τ <dig> =  <dig>  τn =  <dig> and k -  <dig> of the other components equal to  <dig>
now, for the two new vectors τ <dig> and τ, we define the following binary error,
binary error=k0−1−〈τ <dig> τ〉=k0−1−∑i=1n−1τi0τi.
since the two-norm of the vectors involved is fixed, minimizing  is the same as minimizing the euclidean distance between the two vectors or the sum 0– <dig> error.
furthermore, error  is consistent with the russell-rao dissimilarity measure defined on the space of the binary vectors.
its corresponding estimator is
Τ^binerr:=arg⁡min⁡τ′∈ttk <dig> ne=arg⁡max⁡τ′∈ttk <dig> ne.
since we do not know the real value of k <dig>  we should replace it with k^ to compute equation .
doing this, we could amplify the error of the boundary estimation because of the addition of the error of the segment number estimation.
a way to attenuate this issue is to integrate out the number of segments in the conditional expected value.
then the estimator becomes
Τ^binerrak:=arg⁡max⁡τ′∈ttk^,ne.
improved regression curve
as we saw in the previous subsections, there are cases in which the estimation of a parameter of our interest can be made independently of other parameters by integration.
the computation of the brc  and ) suggests to average also over the number of segments by considering the posterior probability of m˜s, given only the sample point y,
μ˜^s:=e.
unfortunately, the computation of this quantity requires time o , hence it could be a problem with samples of big size.
this new type of m˜s estimation is referred to as bayesian regression curve averaging over k .
the same procedure cannot be applied for the level estimation, because in that case we need to know the partition of the whole interval.
properties of the hyper-parameter estimators and definition of new estimators
in order to study the properties of the hyper-parameter estimators defined in equations ,  and , first we need to compute the moments of the data y|ν,ρ <dig> σ <dig>  in the following, we will denote with nq the number of data points in the qth segment.
at first, let us consider only the data which belong to the qth segment.
from the hypothesis of the model, we know that
yj|mq,σ2~nj=tq−1+ <dig> ...,tqmq|ν,ρ2~n,
hence the marginal distribution of any two data points yi and yj belonging to the qth segment is n, where
Σ=.
it follows that the covariance between two data points, which belong to the same segment, is
cov = ρ <dig> i ≠ j,
and
e = ν
var = σ <dig> + ρ <dig>
for each j =  <dig> ...,n, independently of the segment to which it belongs.
furthermore, from the hypotheses of the model, given the segmentation t <dig>  data points belonging to different segments are independent.
expected value and variance of the estimator ν^
the estimator of ν is defined as ν^=y¯ ).
from equation , we can see that this estimator is unbiased and its variance turns out to be
var=σ2n+ρ2∑p=1k0np2n <dig>
hence, the variance is always greater than o, even if we use a denser sampling, i.e.
we augment the number of data points in the interval in which we are estimating the piecewise constant function.
new definition of the estimator σ2^ and its expected value
a circular version of the σ <dig> estimator defined in equation  is
σ2^:=12n∑i=1n <dig>
where yn+ <dig> := y <dig>  using the values of the moments of the data points, its expected value is
e^=σ2+ρ2k0ni{k0≥2},
where we considered two cases in the computation:  when k <dig> =  <dig>  y <dig> and yn belong to the same segment ,  when k <dig> ≥  <dig>  we supposed that the first and the last segments have different levels and so y <dig> and yn are independent.
if the first and the last segments had the same level, then the two segments would be joined together and thus y <dig> and yn would be dependent.
in this case, the expected value would be the same but with k <dig> -  <dig> instead of k <dig>  since the number of segments would be k <dig> -  <dig>  in any case, for k <dig> =  <dig>  the estimator σ2^ is unbiased, while for k <dig> ≪ n but k <dig> ≠  <dig>  σ2^ is almost unbiased.
expected value of the estimator ρ2^
the expected value of the estimator ρ2^ ) is
e=σ2+ρ <dig>
note that when k <dig> =  <dig> , e=σ <dig>  in this degenerate case, the variance of the segment levels ρ <dig> should be estimated with zero but ρ2^ estimates it with the variance of the data points.
moreover, since ∑p=1k0np2≥n , we obtain that
σ2≤e≤.
hence, if n is large the expected value is between σ <dig> and σ <dig> + ρ <dig>  so that, if ρ <dig> ≪ σ <dig>  the estimator is almost unbiased for σ <dig> .
definition of alternative estimator of ρ2: ρ12^
since the covariance between data points belonging to the same segment is ρ <dig>  we could try to use a circular version of the estimator of the autocovariance of a stationary time series
ρ12^:=1n∑i=1n,
where yn+ <dig> := y <dig>  the expected value of the estimator turns out to be
e={−σ2nif k0=1ρ2n−σ2nif k0≥ <dig>
in the computation we considered two cases: k <dig> =  <dig> and k <dig> ≥  <dig>  when k <dig> =  <dig>  y <dig> and yn belong to the same segment and so they are dependent; when k <dig> ≥  <dig>  we suppose that the first and the last segment have not the same level value and so y <dig> and yn are independent.
if k <dig> ≥  <dig> and the first and the last segment had the same level value , then the first and the last segments would be joined together and so y <dig> and yn would be dependent.
in this case, the expected value of the estimator would have the same formula, but with k <dig> -  <dig> instead of k <dig>
we can observe that, when k <dig> =  <dig>  the expected value is negative, while, when k <dig> ≥  <dig>  it can be negative or positive.
moreover, the coefficient of σ <dig> is −1n and so this addendum does not contribute so much to the unbiasedness of the estimator.
the negativity of the expected value happens because the estimator is a generic estimator of the covariance and, in general, this quantity can be negative.
to prevent the negativity of the estimator, we can use its absolute value.
in this way, when the quantity in  is negative, we use the same estimator but with the sign changed in one of the factors of each product, ρ12^=1n∑i=1n.
hence, the meaning of the estimator is the same.
we are interested only in the absolute value of the estimate and not in its sign.
in fact, we already know that the correlation is positive and the negativity of the estimate is due only to the property of the estimator.
our final definition of the estimator is then
ρ12^:=1n|∑i=1n|.
RESULTS
in this section we show and discuss results obtained on both the simulated and the real data.
we used the simulated data with a twofold aim.
the first was to choose empirically the best estimators among those proposed in the previous section.
the second was to compare the original version of bpcr/brc and their modified versions with each other and with other existing methods estimating dna copy number value  <cit> .
on the basis of the results, we selected the best modified version of bpcr, called mbpcr, and of brc.
finally, we compared the performance of mbpcr, cbs, cghseg, glad, hmm, biohmm and rendersome on the real data.
simulation description
in the comparisons, we used several types of artificial data.
we call sample a sequence of data which represents the copy number data of a genomic region, we call dataset a set of samples, while collection a set of datasets.
in order to experimentally evaluate the behavior of all the estimators proposed, we used the artificial datasets sampled from the priors, defined in the hypotheses of the model.
we always chose ν =  <dig> , while we changed the values of σ <dig> and ρ <dig> for each dataset, in order to study different situations of noise .
the most problematic cases were the ones with ρ <dig> <σ <dig> , because in these cases it was hard to identify the true profile of the levels.
we always used n =  <dig>  similar to the mean number of probes of a small chromosome in the affymetrix genechip mapping 10k array , and kmax =  <dig>  in order to have at least  <dig> probes per segment on average.
sometimes we needed datasets where all samples had the same true profile of the segment levels .
this type of dataset is called dataset with replicates.
otherwise, the dataset is called without replicates .
the number of samples per dataset was  <dig>  for datasets with replicates, and  <dig> otherwise.
we considered datasets with replicates in order to be able to compare the goodness of different types of estimations for a given profile.
we also compared the behavior of our boundary estimators using the artificial dataset already employed in  <cit> , where three methods for copy number estimation were examined.
this dataset contained  <dig> samples consisting of  <dig> chromosomes, each of  <dig> probes, which emulated the real copy number data.
this dataset is referred to as simulated chromosomes.
to assess the performance of the several methods, we used three types of artificial datasets.
the first type consisted of four datasets with replicates used in the comparison among the estimators.
this collection of datasets is called cases.
the second type consisted of datasets adapted from the datasets used in  <cit>  to compare several methods for copy number estimation.
in these datasets, each sample was an artificial chromosome of  <dig> probes, where the copy number value was zero apart from the central part where there was an aberration.
the authors considered several widths of aberration:  <dig>   <dig>   <dig> and  <dig> probes.
the noise was always distributed as n, while the signal to noise ratio  was  <dig>   <dig>   <dig> or  <dig>  the snr was defined as the ratio between the height of the aberration and the standard deviation of the noise.
the data of the paper consisted of datasets of  <dig> samples for each combination of width and snr.
we defined our datasets in the following way.
for a fixed snr value, we constructed a chromosome with four aberrations of width of  <dig>   <dig>   <dig> and  <dig> probes, respectively, by joining the corresponding four types of chromosome of the previous datasets.
this collection of datasets is called four aberrations.
in the following, we will consider only the datasets with snr =  <dig>  and snr =  <dig> .
the third type of dataset used was the simulated chromosomes dataset.
comparison among the estimators on simulated data
in this subsection, we present how we selected the best estimators among those proposed in the section methods, on the basis of their results obtained on the artificial datasets.
the comparisons were accomplished using both the true and the estimated values of the other parameters involved in the estimation.
comparison among the hyper-parameter estimators
we applied the hyper-parameter estimators on  <dig> datasets without replicates, considering different values for σ <dig> and ρ <dig>  to evaluate the behavior of the hyper-parameter estimators in all these cases, for each dataset we computed the  mean square error, mse, with respect to the true value of the parameter .
to measure the accuracy of the estimators, we used the estimated mean relative error over all datasets .
the table shows the estimated mean square error of the estimators ν^, σ^ <dig>  ρ^ <dig> and ρ^ <dig> applied to datasets without replicates, for different values of σ <dig> and ρ <dig>  the table shows that mseν^ increased with ρ <dig> and mseσ^ <dig> increased with σ <dig>  the estimator ρ^ <dig> was generally more accurate than ρ^ <dig> with respect to the mse error measure.
from the results, we can deduce that σ^ <dig> is a good estimator because it was quite precise in all situations, while ν^ was sometimes poor but in general acceptable.
about the ρ <dig> estimation, it is better to use ρ^ <dig> than ρ^ <dig>  when the variance of the noise is higher than the variance of the levels.
otherwise, it seems better to use ρ^ <dig> because it does not underestimate ρ <dig> .
comparison among the segment number estimators
we evaluated the quality of the estimators of the number of segments, using datasets with and without replicates for different values of the hyper-parameters σ <dig> and ρ <dig>  the estimations were made using either the true values of the hyper-parameters or the estimated ones.
in this way, we could also observe the behavior of the boundary estimators without the influence of the hyper-parameter estimation.
comparing the absolute, squared and 0– <dig> errors, we found that k^ <dig> generally had the lowest upper bound  of the confidence interval at level 95%, for any type of error and any type of value of the parameter ρ <dig> .
moreover, all estimators always had a similar confidence interval of the 0– <dig> error, while using ρ^ <dig>  in most cases the upper bound of the confidence interval of the absolute and the squared error was lower than using ρ^ <dig>  all these results support the suggestion to use k^ <dig> with ρ^ <dig>
we should also observe that in general k^ <dig> underestimates k <dig>  while k^ <dig> and k^ <dig> overestimate it.
in addition, the percentage of the underestimations increases using ρ^ <dig>
comparison among the boundary estimators
we compared the boundary estimators on the same datasets, previously used for the estimators of the number of segments, with both the estimated and the true value of the parameters involved.
the following errors were taken into account: the sum 0– <dig> error, the joint 0– <dig> error and the binary error, defined in section methods, and the average square error,
aver.
square error:=1k−1∑i=1k−1min⁡j= <dig> ...,k− <dig>
which corresponds to the mean square error over the whole vector of estimated boundaries.
as observed in section methods, when the estimated segment number is used in the estimation of the boundaries, we are only able to compute the binary error, because it does not require that the vector of estimated boundaries has the same length as the vector of the true boundaries.
we found  that the best estimators were t^binerr, t^binerrak and t^joint, but it seemed that t^binerr, t^binerrak were slightly better than t^joint, when we estimated the parameters.
moreover, between these ones, we preferred the latter, because its computation depended less on the estimation of k <dig> and thus it was more stable.
in addition, in case of σ <dig> > ρ <dig>  we observed that there was a great difference between the errors obtained using ρ^ <dig> and ρ^12: using the latter, the errors were significantly lower .
the table shows the estimated mean binary errors ± standard deviation of the estimators of t <dig> applied to datasets with replicates for different values of σ <dig> and ρ <dig>  we always used k^ <dig>  ν^, σ^ <dig> and both ρ <dig> estimators as estimators of the other parameters involved.
the estimators t^binerr and t^binerrak always had the lowest errors.
when σ <dig> > ρ <dig> the error was lower by using ρ^ <dig>  otherwise the two ρ <dig> estimators gave similar errors.
to choose between the estimators t^binerrak and t^joint, we performed the segment level estimation, considering both ρ <dig> estimators.
then, we computed the mse of the estimated segment level per probe and its upper bound of the confidence interval at level 95% .
the results showed that, for a given estimator of ρ <dig>  in general t^binerrak was better than t^joint and the upper bound of the confidence interval of the mse of the former estimator was lower than that one of the latter.
moreover, using ρ^ <dig> the error was generally lower.
in addition, we compared the behavior of our boundary estimators on dataset simulated chromosomes.
to assess the goodness of their estimation, we measured the sensitivity  and the false discovery rate , while to assess the influence of the boundary estimation on the profile estimation, we calculated the sum of squared distance , the median absolute deviation  and the accuracy .
the sensitivity and the fdr were computed not only looking at the exact position of the breakpoints , but also accounting for a neighborhood of  <dig> or  <dig> probes around the true positions .
we also computed the accuracy inside and outside the aberrations separately, since the samples of dataset simulated chromosomes contained only few small copy number changes and thus the accuracy depended more on the correct estimation/classification of the probes in the "normal" regions.
a more detailed explanation of these measures can be found in  <cit> .
since we estimated σ^ <dig> =  <dig>  and ρ^ <dig> =  <dig>  , we expected that using ρ^ <dig> we should obtain better results because σ <dig> was lower than ρ <dig> and indeed the results obtained with ρ^ <dig> were slightly better than the others .
moreover, we found in general that t^binerrak had the highest sensitivity but also a higher fdr than t^ <dig> .
this was due to the fact that, although the estimated number of segments was the same, t^ <dig> could estimate some breakpoints with the same position, reducing the total number of breakpoints and so reducing the fdr.
we can see in table  <dig> that the false estimated breakpoints did not negatively influence the profile estimation.
in fact, the false breakpoints are often used by the algorithm in two ways: either to divide a long segment into two or more segments with close levels or, if it is difficult to determine the position of a breakpoint, to add, before or after the aberration, a segment of one point with the value of the level between zero and the aberration level.
overall, t^binerrak with ρ^ <dig> performed best on this dataset.
this table shows the comparison among the error measures for profile estimation, obtained on dataset simulated chromosomes, for all estimators of t <dig> .
we always used k^ <dig>  ν^, σ^ <dig> and both ρ <dig> estimators as estimators of the other parameters involved.
the estimator t^binerrak always had the lowest ssq and mad errors and the highest accuracy both inside and outside the regions of aberration.
using ρ^ <dig> the performance was slightly better, because σ <dig> <ρ <dig>
in conclusion, we suggest to use t^binerrak, even if t^joint is also a quite good estimator in some cases.
regarding the estimation of ρ <dig>  it seems that it is better to use ρ^ <dig> in presence of high noise.
comparison among the regression curves
we compared the estimation of the levels of brc with the one of brcak, also taking into account the influence of the different estimators of the parameters on the final results.
to valuate the performance of the methods, we used the root mean square error  per probe and per sample, computed with respect to the true profile of the levels.
for this purpose, we necessarily needed datasets with replicates.
using brcak we generally obtained a better or equal result with respect to the brc .
moreover, we observed that, using brc, it was better to estimate the segment number with k^ <dig> or k^ <dig>
note that we still have to solve the problem to determine which is the best estimator of ρ <dig>  in most cases, the profile obtained by using ρ^ <dig> was better than using ρ^ <dig> .
this is due to the fact that sometimes ρ^ <dig> slightly underestimated ρ <dig>  leading to overfitting.
still we recommend to use ρ^ <dig>  even if it could lead to a slight overfitting especially in the case of few segments.
comparison with other methods on simulated data
in this subsection we compare the original and modified versions of bpcr and brc, with other existing methods for genomic copy number estimation  <cit> .
error measures used in the comparison
we used two different measures to examine the behavior of the different methods on the collections cases and four aberrations: the root mean square error  and the roc curve .
each point of the roc curve has as coordinates the false positive rate  and the true positive rate  for a certain threshold.
the tpr is defined as the fraction of probes in the true aberrations whose estimated value is above the threshold considered, while the fpr consists in the fraction of probes outside the true aberrations whose estimated value is above the threshold.
hence, the roc curve measures the accuracy of the method in the detection of the true aberrations.
instead, the evaluation of the several methods on dataset simulated chromosomes was accomplished using the error measures described in  <cit> , already used in the study of the boundary estimators.
before showing the results, we need to remember that some methods estimate the copy numbers as a piecewise constant function, while other algorithms estimate them as a continuous curve.
hence, bpcr was compared to the former group of methods, while the bayesian regression curves to the latter.
since some error measures of  <cit>  suppose that the estimated profile is piecewise constant, we applied only the former group of methods on dataset simulated chromosomes.
the piecewise constant estimation
we compared the original and the modified versions of bpcr with cbs  <cit> , cghseg  <cit> , glad  <cit> , hmm  <cit> , biohmm  <cit>  and rendersome  <cit> .
for thoroughness, in the modified versions of bpcr, we used both ρ^ <dig> and ρ^ <dig> as estimators of ρ <dig>  k^ <dig> as estimator of k <dig> and both t^joint and t^binerrak as estimators of the boundaries.
we used ρ^ <dig> when the noise was low  and otherwise ρ^ <dig>  examples of estimated profiles can be found in figure  <dig>
in general, in presence of medium noise, the glad method performed worst, since it had a high error in the level estimation of the small peaks, while, for high noise, often both glad and rendersome failed to detect the aberrations .
the cghseg method did not usually exhibit an appropriate level estimation except sometimes for segments of large width .
this is due to the fact that cghseg estimates the level of a segment with the arithmetic mean of the data points in the segment and this estimator performs poorly if the segment contains few data points and the breakpoint estimation is not accurate.
the cbs method, in general, performed quite well, but it was unable to detect aberrations of small width, especially when the noise was high .
on the collection cases and the dataset of four aberrations with snr =  <dig>  the rmse plots and the roc curves of the hmm method showed that it generally estimated the profile well, but sometimes it exhibited high errors near breakpoint positions , likely because it was unable to determine the true position of the breakpoints precisely.
moreover, on the dataset with snr =  <dig>  we recognized the true issues of the estimation with hmm.
the rmse plot showed that it had a high error outside the regions of the aberrations, while inside these regions the error was always more or less the same.
hence, it often failed also in the estimation of the largest aberration, the easiest one to detect .
the reason of this behavior of the rmse is the following.
the method estimated the true profile either with only one segment, or more often with a profile consisting of a lot of small segments, but all with the same level.
since in the latter case the estimated levels were close to that one of the true aberrations, the rmse was low in the regions of the aberrations but high outside them.
however, the estimated profiles were not similar to the true one.
in presence of medium noise , the method biohmm was more precise than hmm in the determination of the breakpoint positions and in the level estimation , while for high noise it behaved similarly to hmm .
in general we found that, when σ <dig> <ρ <dig> , the version of bpcr with t^binerrak and ρ^ <dig>  generally gave the best estimation compared to the other versions of bpcr and to the other methods .
in the following, we will call this modified version of bpcr, mbpcr.
only on the dataset with snr =  <dig>  we could not choose the "best" method, because this case showed the limits of all the methods considered.
the problem regarding the modified versions of bpcr was essentially the estimation of the number of segments.
the roc curves  of the modified versions of bpcr with ρ^ <dig> were the closest to the left and the top sides of the box, while the rmse plot  showed that these methods were the best methods in the estimation of the levels inside the aberrations, but not outside them.
in general, in case of very high noise, all the modified versions of bpcr well detected the aberrations, but had problems in the estimation of the profile outside them because of the poor estimation of the number of segments.
in fact, k^ <dig> tends to overestimate the number of segments and this problem worsens using ρ^ <dig>  in conclusion, in a situation with very high noise, using ρ^ <dig> the bpcr methods detect better the small segments, but, at the same time, the large ones are divided in small segments.
on the other hand, using ρ^ <dig>  smaller segments are not detected and are joined to the closest large segment.
finally, the comparison performed on dataset simulated chromosomes showed that cbs and mbpcr better estimated the profiles .
regarding the breakpoint error measures , we found that mbpcr had the highest sensitivity , but also a higher fdr than cbs.
we have already explained in the previous subsection the possible reason of the high fdr of mbpcr and we can observe again that this fact did not influence negatively the profile estimation .
the glad method showed a low sensitivity and low fdr, apart from the case regarding the exact position of the breakpoints , which means that it underestimated the segment number and the estimated breakpoints were not located exactly at their true positions.
also cghseg underestimated the number of segments because of low sensitivity and fdr, while hmm had low sensitivity and high fdr when w =  <dig> and vice versa in the other cases, which means that it often detected the true segment number, but it was unable to put the breakpoints at their exact position.
instead, biohmm solved the issue of hmm with w =  <dig>  but overall had a lower sensitivity than hmm.
rendersome missed several true aberrations  and detected some false aberrations .
the table shows the comparison of the level estimations obtained using several piecewise constant methods on dataset simulated chromosomes.
in this comparison, the methods cbs and mbpcr exhibited the lowest ssq error in the profile estimation and the highest accuracy inside the aberrated regions.
on the other hand, hmm, biohmm and rendersome had the highest accuracy outside the aberrations, but a high ssq error.
therefore, the former group of algorithms globally estimated a better profile than the latter.
because of its definition, the mad error is less informative: it does not take into account if a small number of probes are wrongly estimated, but these probes could correspond to breakpoints or small aberrations.
estimation with a continuous curve
we compared the several versions of the bayesian regression curves with methods which estimate the copy number as a continuous curve .
lowess is the acronym of "locally weighted smoothing"  and it is one of the methods considered in the comparison performed in  <cit> .
as we saw previously, both the brc, which uses k^ <dig>  and the brcak perform well, so we tested both versions with both estimators of ρ <dig>  figure  <dig> shows examples of estimated profiles with these smoothing methods.
in general, we found that all methods detected the regions of aberration quite well .
the wavelet method showed a higher error in the level estimation of the aberrations in the datasets snr =  <dig> and snr =  <dig> .
the methods lowess and quantreg had the highest rmse in the collection cases, while their error was not significantly different outside and inside the aberrations on datasets with snr =  <dig>   <dig>  therefore, in the last cases the error was low inside the aberrations and high outside them in comparison with the other methods.
the method smoothseg showed a similar behavior, but with a lower error.
moreover, we found that the roc measure was affected by oscillations in the estimated curve, which lead to roc curves intersected and difficult to be interpreted .
this complex behavior is a consequence of the way in which lowess, wavelet, quantreg and smoothseg yielded oscillating curves with positive and negative values outside the aberrations; while brcs estimated the true profile with a line almost flat and close to zero .
in conclusion, the version of brc with k^ <dig> and brcak gave in general a better estimation than the other brcs and the other smoothing methods considered.
regarding the ρ <dig> estimation, we found that it is better to use ρ^ <dig>  if σ <dig> <ρ <dig>  and ρ^ <dig>  if σ <dig> > ρ <dig>
application to real data
in this subsection, we show how mbpcr performed compared to other piecewise constant estimation methods on the real data.
we used samples from three mantle cell lymphoma cell lines  previously analyzed by us with the affymetrix genechip mapping 10k array , an oligonucleotides-based microarray  <cit> .
we also used the data obtained on jeko- <dig>  by using the higher density affymetrix genechip mapping 250k nsp array.
we considered eight recurrent gene regions of aberration in lymphoma plus other two gene regions  and we compared the corresponding copy numbers obtained by the several piecewise constant methods with those obtained by the fish technique in  <cit> .
in the end, we also show a comparison among the estimated profiles of chromosome  <dig> of jeko- <dig>
the 10k array data used are freely available at the public repository gene expression omnibus  <cit>  with geo accession: gsm <dig>  gsm <dig> and gsm <dig>  the 250k array data of cell line jeko- <dig> will be soon available in the same repository.
with the current implementation, on a computer with dual cpu  and  <dig> gb ram, the algorithm needed about  <dig> minutes to analyze a 10k array sample, while about  <dig> day to estimate the profile of a 250k array sample.
the computations were performed by chromosome  and using kmax =  <dig>
gene copy number estimation
to properly evaluate the methods, the knowledge of the true underling profile is required.
in general, large aberrations on chromosomes can be detected with conventional karyotype analysis or with sky-fish and one could use this information for the evaluation procedure, but the width of these aberrations is so large that all the methods can detect them well, leading to a useless comparison.
for this reason, we decided to take into account only genes to compare the piecewise constant methods.
in the comparison, as previously published  <cit> , when two fish copy numbers had been assigned to one gene, the first number should correspond to the copy number detected in the majority of the cells.
we assigned two estimated copy numbers to one gene, when the position of the gene is between two snps and the method assigned two different values to these snps.
the results on rec- <dig>  did not show any significant difference among the methods, instead those on granta- <dig>  showed that glad was unable to detect the true copy number in five cases, while hmm, biohmm and rendersome detected an amplification on malt <dig> greater than what detected by fish analysis.
all methods did not detect the true copy number of atm, probably because the snps around atm are far away from the corresponding fish region  and the deletion affects only this region.
only mbpcr with ρ^ <dig> and hmm detected a breakpoint between the two snps around atm region, indicating a copy number change.
regarding the jeko- <dig> data, since the cell line is triploid, to obtain more realistic copy number value, we centered the estimated log2ratio around log <dig>  <dig>  with the denser 250k array data, all methods behaved equally good.
only hmm had a problem in the detection of the breakpoint corresponding to the c-myc amplification .
on both arrays, all methods identified a gain  at the ccnd <dig> position, while the copy number detected by fish is  <dig>  this fact cannot be explained as previously for atm, because this region is well covered by snps.
instead, on the jeko- <dig> 10k array data , the noisiest among all samples, we can see several cases in which cbs, hmm and glad did not detect correctly the gene copy number .
this occurred more frequently to biohmm and rendersome, while only once to cghseg .
the method mbpcr with ρ^ <dig> always estimated gene copy numbers correctly, apart from ccnd <dig>
on this noisy data, biohmm and rendersome often estimated the gene copy number wrongly, while this occurred only sometimes to cbs, hmm and glad.
the method mbpcr with ρ^ <dig> correctly estimated the gene copy numbers, apart from ccnd <dig> whose copy number was estimated by all methods differently from the fish technique.
profile estimation
to compare the profile estimations, we chose the sample jeko- <dig> because, using the results obtained on both types of array, we could at least understand which regions were more realistically estimated.
for now, whole validated chromosomic profiles are not available.
among all chromosomes, we chose chromosome  <dig> since three of the previous genes belong to that: ccnd <dig> , birc <dig>  and atm .
from the graphs in figure  <dig> we can observe that, among all the piecewise constant methods, only mbpcr with ρ^ <dig> was able to detect the high amplification after position  <dig> mb on the 10k array data, while it was recognized by all methods  on the 250k array data.
moreover, on the 10k array data, almost all methods detected a false deletion around position  <dig> mb, due to the presence of a sequence of outliers, and biohmm did not find any copy number change in the chromosome.
on the 250k array data, hmm and rendersome had problems in recognizing the last part of the chromosome as a flat region.
moreover, on the 10k array data, rendersome estimated several outliers as true aberrations and, on the 250k array data, it was unable  to identify the whole region from about  <dig> mb to  <dig> mb as gained.
CONCLUSIONS
we introduced new estimators for the parameters involved in bpcr and we selected the best ones on the basis of theoretic and empirical results.
in particular, we found that the best way is to estimate the segment number with k^ <dig> and the boundaries with t^binerrak .
we call mbpcr the bpcr version which uses k^ <dig> and t^binerrak.
concerning the estimation of the variance of the segment levels, we found that the original estimator ρ^ <dig> overestimates ρ <dig>  by an addendum proportional to σ <dig> , see equation .
hence, the estimation fails when σ <dig> > ρ <dig>  the new estimator ρ^ <dig> is more precise but slightly underestimates ρ <dig>  leading to an overestimation of the segment number.
applying both estimators on artificial datasets, we found that, in general, the best way is to use ρ^ <dig> when σ <dig> <ρ <dig> , but to use ρ^ <dig> when σ <dig> > ρ <dig> , even if it could lead to a slight overfitting.
on real dna copy number data, commonly σ <dig> > ρ <dig>
we also compared mbpcr with other methods which also estimate the copy number as a piecewise constant function: cbs, hmm, cghseg, glad, biohmm and rendersome.
as a whole, the results showed that mbpcr gave the best estimation on the dataset used.
however, when σ <dig> ≫ ρ <dig> it is hard to understand which method is the most appropriate.
most of the other methods were not able to detect aberrations with a small width  and the same was true for mbpcr using ρ^ <dig>  on the other hand, the use of ρ^ <dig> led to the detection of the smaller segments, but the larger ones were often divided in small segments and sometimes the segments consisted of only one point.
the optimal choice of the ρ <dig> estimator is still not fully determined.
the new estimators improved also brc, which is a bayesian regression with a smooth curve.
moreover, we derived a formula to estimate brc without employing the estimated number of segments.
we referred to it as brcak.
applying these methods to artificial data, the best estimators were found to be the brc version with k^ <dig> and brcak.
about the choice of the two estimators of ρ <dig>  we found a similar conclusion as before, with the advantage of being less problematic when σ <dig> ≫ ρ <dig>  we compared these two regression methods with other methods which estimate the copy number data as a continuous curve: wavelet, lowess, quantreg and smoothseg.
the results showed that our modified regression methods were the most appropriate for the estimation of the segment levels on the datasets considered.
even though these smoothing methods seem to have less problems in the estimation and the error measures  suggest that their estimation is even better than the piecewise constant estimation, these methods do not detect the position of the breakpoints explicitly and hence the changes in the value of the segment level.
thus, they seem to be less adequate in practice.
for this reason, we feel that the roc curve cannot be used as the only measure to compare methods, as previously done in  <cit> .
the rmse is generally an acceptable measure, but we observed that in some cases even this is not sufficient, due to the overfitting.
willenbrock and fridlyand  <cit>  proposed other measures to compare methods for copy number estimation, regarding both breakpoint and level estimation.
in particular, the sensitivity measure of breakpoint estimation is useful to select which methods should be used, because it quantifies the precision of the methods in determining the position of the breakpoints.
finally, we have applied mbpcr and all other piecewise constant regression methods to real data.
the comparisons showed that mbpcr estimated well the copy number of the genes.
on these data, in most cases the choice of the ρ <dig> estimator did not affect the analysis.
in comparison with the other methods, the current implementation of our algorithm is computationally intensive.
the real computational time can be reduced linearly diminishing kmax and quadratically diminishing the number of probes.
moreover, the computation can be easily parallelized by arm and by chromosome, reducing further the calculation time.
in cancer research, the accuracy in the dna copy number estimation is crucial for the correct determination of the mutations that characterize the disease.
in particular, the estimation of the breakpoints must be precise to detect correctly which genes are affected by these genomic aberrations.
as recently shown  <cit> , snp microarrays can also potentially detect the breakpoints involved in unbalanced translocations, allowing the identification of fusion genes .
in this context, the use of our method can highly improve the disease investigation, because it accurately determines breakpoints, is less sensitive to high noise and generally outperforms all the methods considered in our comparisons.
moreover, smoothing algorithms are clearly not suitable for such analysis.
availability and requirements
project name: mbpcr.
project home page: .
operating systems: the software should run in linux, mac-os or windows.
tests were performed on windows and linux systems.
programming language: r.
other requirements: none.
licence: gnu gpl.
any restrictions to use by non-academics: none.
authors' contributions
pmvr carried out the study and wrote the manuscript.
mh and ik supervised the statistical analysis.
fb supervised the validation study and provided the real data.
all authors read and approved the final manuscript.
supplementary material
additional file 1
mbpcr source code.
this zipped file contains the source code of the mbpcr algorithm in r, including help files, sample data and examples.
click here for file
additional file 2
supplementary material.
this file contains a brief description of the dynamic program used to implement the method, a formal definition of some error measure used, some supplementary tables and some supplementary figures.
click here for file
acknowledgements
this work was supported by the grant 205321- <dig> of the swiss national science foundation and by oncosuisse grant ocs 1939-08- <dig>

 lesions at dna level represent the cause of cancer and of many congenital or hereditary disorders. the change of the number of copies of dna in a genomic region is one of the most common aberrations. in normal cells each genomic segment is present in two copies, but, for example, in tumor cells the genome can present regions of deletions, gains or amplifications. thus, in general, the dna copy number along the genome can be represented as a piecewise constant function. with microarray technology it is possible to simultaneously measure the copy number along the genome at hundred thousands of positions. however, raw copy number data are generally very noisy. hence, it is important to define a method which allows one to estimate the number of regions with different copy number, the endpoints of these regions and their copy number. several methods have been developed to solve this issue. many methods consider the log ratio of the data and model it as a normal random variable, since they assume that the errors are normally distributed. we can roughly subdivide all of these methods into two classes: those ones that estimate the copy numbers as a piecewise constant function and the others that estimate the copy numbers as a continuous curve. the methods belonging to the latter group are called smoothing methods. among the methods belonging to the first class, we can find the following. the circular binary segmentation approach is a recursive method in which the breakpoints are determined on the basis of a test of hypothesis, with null hypothesis that in the interval considered there is no change in copy number. picard et al. used a piecewise constant regression model, where the parameters are estimated maximizing a penalized likelihood. this method is usually denoted with the abbreviation cghseg. the glad method is another piecewise constant regression method, but in this case the parameters are estimated maximizing a weighted likelihood. fridlyand et al. applied hidden markov models, while marioni et al. defined an hmm method which takes into account the distance among the data points. recently, nilsson et al derived a segmentation method based on total variation minimization, called rendersome. it is optimized for gene expression data, but the authors affirm that it can be used also on copy number data. among the smoothing methods, hsu et al. used a wavelet regression method with haar wavelet. eilers and de menez applied a quantile smoothing regression, with the solution found by minimizing a loss function based on the l norm, to obtain a flatter curve. huang et al. proposed smoothseg, i.e. a smooth segmentation method based on a doubly heavy-tailed random-effect model. we propose a piecewise constant regression method, using bayesian statistics, which appears appropriate when regions contain only few data points. the original version of the method was presented by hutter. in this paper we propose improved bayesian estimators of the parameters involved and we apply the model to dna copy number estimation. finally, we compare our algorithm with some among the most cited or more recent methods, on artificial and real data. our method was implemented in r and is freely available at or in additional file furthermore, an r package will be soon available. methods in the first two subsections, we briefly describe the original bayesian piecewise constant regression method, explaining the hypothesis of the model and the estimation of its parameters with bayesian regression. we emphasize the definitions of the original parameter estimators in order to show how we changed some of these estimators in the other subsections. a brief explanation of the dynamic program for the computation of the estimators can be found in the additional file. regarding notations, we do not indicate explicitly the random variable to which a distribution is referred, if it is clear from the context. for example, pk ≡ p or fy, m ≡ f. hypotheses of the model let y ∊ ℝn be a random vector, such that each component is conditionally normally distributed: yi|μ˜i σ2~n. suppose also that y represents a noisy observation of a piecewise constant function, which consists of k horizontal segments. then, the segment level at a generic position i does not assume different values for each i, but the data are divided into k intervals where μ˜tq−1+10=...=μ˜tq0=:μq for each q =...,k hence, μq represents the level of the qth segment. our goal is to estimate the levels μ0= of all the segments. in order to do that, we first need to estimate the number of the segments k and the partition of the data t from a bayesian point of view, μ t and k are treated as random variables, hence we denote them with the corresponding upper case letters. moreover, because of their randomness, we need to define a prior distribution for each of them to complete the model. for the number of segments and the boundaries, we assume noninformative prior distributions: p=1kmax⁡k∈k p=1t∈tk,n, where k = {...,kmax} and tk,n is the subspace of ℕ0k+ such that t = tk = n and tq ∊ {...,n - 1} for all q =...,k - in an ordered way and without repetitions. about m, we assume that all its components are mutually independent and identically normally distributed, m|ν,ρ k=k~n, where ν ∊ ℝk, such that νq = ν for each q =...,k, and i ∊ ℝk × k, such that ip, q = δp,q for each p, q =...,k. instead of these assumptions, we could assume a cauchy distribution for each yi or mq in order to model an observation whose noise has heavier tails, as previously done by hutter. original estimation: the bpcr method the statistical procedure consists in a sequence of estimations due to the relationship among the parameters. bpcr estimates the number of segments with the map estimate given the sample point y, k^:=arg⁡max⁡k∈kp, and, given k^, also each boundary is estimated separately with its corresponding map estimate, t^p:=arg⁡max⁡h∈{p,...,n−}p for all p =...,k^ - finally, the rth moment of the level of the mth segment is estimated with its posterior mean. since its computation needs the knowledge of the number of segments and the partition of the data, we replace them with the estimated ones, μ^mr:=e for all m =...,k^. when r = and we assume that y and m are normally distributed, the estimate turns out to be μ^m=ρ2∑i=tm−1+1tmyi+σ2νρ2+σ for all m =...,k^. when the sample contains only one segment, the bayesian estimation of the posterior distribution of the levels should theoretically lead to a normal distribution, similar to a dirac delta function centered at ν^, since the levels can assume only one value from the data. in fact, in this case, if we estimate ρ only using the data, then this value will be close to zero and so the level will be estimated with ν^, the mean of the data ). we can estimate the segment level m˜s at a generic position s, using the fact that it belongs to some segment m and in this segment it is equal to the corresponding mm. then, summing over all the possible segments, we can compute its posterior distribution in the following way: f=∑m=1k0∑i=0s−1∑j=snf and the corresponding estimate of m˜s given k^ is μ˜^s:=e, for all s =...,n. the vector μ˜^ is called bayesian regression curve. the probability distributions defined previously depend on the hyper-parameters ν, ρ and σ. hutter suggested the following estimators: ν^:=1n∑i=1nyi=y¯ ρ^2:=1n∑i=1n σ^2:=12∑i=1n− improved estimators of the number of segments to understand the real meaning of the map estimator k^, we need to introduce the theory of the construction of a generic bayesian estimator. in general, a bayesian estimator is defined in the following way. let us suppose that z is a random variable whose distribution depends on an unknown parameter θ, which we want to estimate. since we cannot exactly know the true value of the parameter, we consider it as a random variable Θ with a given prior probability distribution. in order to measure the goodness of the estimation, we define an error and we choose the estimator that minimizes the expected error given the sample z, Θ^:=arg⁡min⁡θ′e. the 0– error is commonly used for a parameter which can assume only a discrete number of values. the estimator corresponding to this error is the map estimator, arg⁡min⁡θ′e=arg⁡max⁡θ′∑θδθ,θ′p=arg⁡max⁡θ′p. obviously, if we use different types of errors, we can obtain different estimators. in the following, we will use k^ to denote any estimator of k, while k^ to denote the parameter estimator k^ based on the 0– error. using the 0– error, we give the same penalty to each value different from the true value, whether it is close to or far away from the true one. to take into account the distance of the estimated value from the true one, we need to use other types of errors, which are based on different definitions of distance, such as, absolute error:= |θ - θ'| squared error:= if the parameter θ ∊ ℝ, then the estimators corresponding to these errors are the median and the mean of its posterior distribution, respectively. in our case, we denote these estimators of k with k^ and k^ improved estimators of the boundaries similarly to the previous subsection, we derive alternative boundary estimators, by considering different types of errors. we denote the map boundary estimator defined in equation with t^ meaning of the estimator t^01 the estimator t^ is defined in such a way that each component minimizes the 0– error of the corresponding boundary, separately. explicitly, given the sample point y and the segment number k its estimate is t^01=,...,arg⁡max⁡tk0−1∈tp,n), where t = {...,n - 1}. t^ may be regarded as an approximation of the bayesian estimator that minimizes the error which counts the number of wrongly estimated boundaries: sum 0-1 error=∑p=1k0−1=k0−1−∑p=1k0−1δtp tp, that is t^sum=arg⁡max⁡t∈tk n∑p=1k0−1p. definition of the estimator Τ^joint a problem of the latter estimator ) is its computational complexity, because it needs the computation of all the ordered combinations of the boundaries. on the other hand, there are two reasons for which t^ is not a suitable estimator of the boundaries. first, it does not take into account that the boundaries are dependent, because they have to be ordered, and second, in principle, it can have more than one component with the same value. as a consequence, a theoretically more correct way to estimate the boundaries is minimizing the 0– error with respect to the joint boundary probability distribution. then, given k and y, the boundary estimator becomes t^joint=arg⁡max⁡t∈tk np. definition of the estimators Τ^binerr and Τ^binerrak we must notice that the estimators considered until now have the same length of the true vector of the boundaries. in practice, the number of segments k is unknown, so that we should use k^. as a consequence, if k^ is different from k then, strictly speaking, we cannot minimize the previous types of error because the vectors have different length. a way to solve this issue is to map each boundary vector into a vector τ∈ℝ0n+ in the following way: t↦τ such that τi={1if ∃p such that tp=i0otherwise. we denote with ttk,n the set of all the possible τ with τ = τn = and k - of the other components equal to now, for the two new vectors τ and τ, we define the following binary error, binary error=k0−1−〈τ τ〉=k0−1−∑i=1n−1τi0τi. since the two-norm of the vectors involved is fixed, minimizing is the same as minimizing the euclidean distance between the two vectors or the sum 0– error. furthermore, error is consistent with the russell-rao dissimilarity measure defined on the space of the binary vectors. its corresponding estimator is Τ^binerr:=arg⁡min⁡τ′∈ttk ne=arg⁡max⁡τ′∈ttk ne. since we do not know the real value of k we should replace it with k^ to compute equation. doing this, we could amplify the error of the boundary estimation because of the addition of the error of the segment number estimation. a way to attenuate this issue is to integrate out the number of segments in the conditional expected value. then the estimator becomes Τ^binerrak:=arg⁡max⁡τ′∈ttk^,ne. improved regression curve as we saw in the previous subsections, there are cases in which the estimation of a parameter of our interest can be made independently of other parameters by integration. the computation of the brc and ) suggests to average also over the number of segments by considering the posterior probability of m˜s, given only the sample point y, μ˜^s:=e. unfortunately, the computation of this quantity requires time o, hence it could be a problem with samples of big size. this new type of m˜s estimation is referred to as bayesian regression curve averaging over k. the same procedure cannot be applied for the level estimation, because in that case we need to know the partition of the whole interval. properties of the hyper-parameter estimators and definition of new estimators in order to study the properties of the hyper-parameter estimators defined in equations, and, first we need to compute the moments of the data y|ν,ρ σ in the following, we will denote with nq the number of data points in the qth segment. at first, let us consider only the data which belong to the qth segment. from the hypothesis of the model, we know that yj|mq,σ2~nj=tq−1+...,tqmq|ν,ρ2~n, hence the marginal distribution of any two data points yi and yj belonging to the qth segment is n, where Σ=. it follows that the covariance between two data points, which belong to the same segment, is cov = ρ i ≠ j, and e = ν var = σ + ρ for each j =...,n, independently of the segment to which it belongs. furthermore, from the hypotheses of the model, given the segmentation t data points belonging to different segments are independent. expected value and variance of the estimator ν^ the estimator of ν is defined as ν^=y¯ ). from equation, we can see that this estimator is unbiased and its variance turns out to be var=σ2n+ρ2∑p=1k0np2n hence, the variance is always greater than o, even if we use a denser sampling, i.e. we augment the number of data points in the interval in which we are estimating the piecewise constant function. new definition of the estimator σ2^ and its expected value a circular version of the σ estimator defined in equation is σ2^:=12n∑i=1n where yn+:= y using the values of the moments of the data points, its expected value is e^=σ2+ρ2k0ni{k0≥2}, where we considered two cases in the computation: when k = y and yn belong to the same segment, when k ≥ we supposed that the first and the last segments have different levels and so y and yn are independent. if the first and the last segments had the same level, then the two segments would be joined together and thus y and yn would be dependent. in this case, the expected value would be the same but with k - instead of k since the number of segments would be k - in any case, for k = the estimator σ2^ is unbiased, while for k ≪ n but k ≠ σ2^ is almost unbiased. expected value of the estimator ρ2^ the expected value of the estimator ρ2^ ) is e=σ2+ρ note that when k =, e=σ in this degenerate case, the variance of the segment levels ρ should be estimated with zero but ρ2^ estimates it with the variance of the data points. moreover, since ∑p=1k0np2≥n, we obtain that σ2≤e≤. hence, if n is large the expected value is between σ and σ + ρ so that, if ρ ≪ σ the estimator is almost unbiased for σ. definition of alternative estimator of ρ2: ρ12^ since the covariance between data points belonging to the same segment is ρ we could try to use a circular version of the estimator of the autocovariance of a stationary time series ρ12^:=1n∑i=1n, where yn+:= y the expected value of the estimator turns out to be e={−σ2nif k0=1ρ2n−σ2nif k0≥ in the computation we considered two cases: k = and k ≥ when k = y and yn belong to the same segment and so they are dependent; when k ≥ we suppose that the first and the last segment have not the same level value and so y and yn are independent. if k ≥ and the first and the last segment had the same level value, then the first and the last segments would be joined together and so y and yn would be dependent. in this case, the expected value of the estimator would have the same formula, but with k - instead of k we can observe that, when k = the expected value is negative, while, when k ≥ it can be negative or positive. moreover, the coefficient of σ is −1n and so this addendum does not contribute so much to the unbiasedness of the estimator. the negativity of the expected value happens because the estimator is a generic estimator of the covariance and, in general, this quantity can be negative. to prevent the negativity of the estimator, we can use its absolute value. in this way, when the quantity in is negative, we use the same estimator but with the sign changed in one of the factors of each product, ρ12^=1n∑i=1n. hence, the meaning of the estimator is the same. we are interested only in the absolute value of the estimate and not in its sign. in fact, we already know that the correlation is positive and the negativity of the estimate is due only to the property of the estimator. our final definition of the estimator is then ρ12^:=1n|∑i=1n|.  in this section we show and discuss results obtained on both the simulated and the real data. we used the simulated data with a twofold aim. the first was to choose empirically the best estimators among those proposed in the previous section. the second was to compare the original version of bpcr/brc and their modified versions with each other and with other existing methods estimating dna copy number value. on the basis of the results, we selected the best modified version of bpcr, called mbpcr, and of brc. finally, we compared the performance of mbpcr, cbs, cghseg, glad, hmm, biohmm and rendersome on the real data. simulation description in the comparisons, we used several types of artificial data. we call sample a sequence of data which represents the copy number data of a genomic region, we call dataset a set of samples, while collection a set of datasets. in order to experimentally evaluate the behavior of all the estimators proposed, we used the artificial datasets sampled from the priors, defined in the hypotheses of the model. we always chose ν =, while we changed the values of σ and ρ for each dataset, in order to study different situations of noise. the most problematic cases were the ones with ρ <σ, because in these cases it was hard to identify the true profile of the levels. we always used n = similar to the mean number of probes of a small chromosome in the affymetrix genechip mapping 10k array, and kmax = in order to have at least probes per segment on average. sometimes we needed datasets where all samples had the same true profile of the segment levels. this type of dataset is called dataset with replicates. otherwise, the dataset is called without replicates. the number of samples per dataset was for datasets with replicates, and otherwise. we considered datasets with replicates in order to be able to compare the goodness of different types of estimations for a given profile. we also compared the behavior of our boundary estimators using the artificial dataset already employed in, where three methods for copy number estimation were examined. this dataset contained samples consisting of chromosomes, each of probes, which emulated the real copy number data. this dataset is referred to as simulated chromosomes. to assess the performance of the several methods, we used three types of artificial datasets. the first type consisted of four datasets with replicates used in the comparison among the estimators. this collection of datasets is called cases. the second type consisted of datasets adapted from the datasets used in to compare several methods for copy number estimation. in these datasets, each sample was an artificial chromosome of probes, where the copy number value was zero apart from the central part where there was an aberration. the authors considered several widths of aberration: and probes. the noise was always distributed as n, while the signal to noise ratio was or the snr was defined as the ratio between the height of the aberration and the standard deviation of the noise. the data of the paper consisted of datasets of samples for each combination of width and snr. we defined our datasets in the following way. for a fixed snr value, we constructed a chromosome with four aberrations of width of and probes, respectively, by joining the corresponding four types of chromosome of the previous datasets. this collection of datasets is called four aberrations. in the following, we will consider only the datasets with snr = and snr =. the third type of dataset used was the simulated chromosomes dataset. comparison among the estimators on simulated data in this subsection, we present how we selected the best estimators among those proposed in the section methods, on the basis of their results obtained on the artificial datasets. the comparisons were accomplished using both the true and the estimated values of the other parameters involved in the estimation. comparison among the hyper-parameter estimators we applied the hyper-parameter estimators on datasets without replicates, considering different values for σ and ρ to evaluate the behavior of the hyper-parameter estimators in all these cases, for each dataset we computed the mean square error, mse, with respect to the true value of the parameter. to measure the accuracy of the estimators, we used the estimated mean relative error over all datasets. the table shows the estimated mean square error of the estimators ν^, σ^ ρ^ and ρ^ applied to datasets without replicates, for different values of σ and ρ the table shows that mseν^ increased with ρ and mseσ^ increased with σ the estimator ρ^ was generally more accurate than ρ^ with respect to the mse error measure. from the results, we can deduce that σ^ is a good estimator because it was quite precise in all situations, while ν^ was sometimes poor but in general acceptable. about the ρ estimation, it is better to use ρ^ than ρ^ when the variance of the noise is higher than the variance of the levels. otherwise, it seems better to use ρ^ because it does not underestimate ρ. comparison among the segment number estimators we evaluated the quality of the estimators of the number of segments, using datasets with and without replicates for different values of the hyper-parameters σ and ρ the estimations were made using either the true values of the hyper-parameters or the estimated ones. in this way, we could also observe the behavior of the boundary estimators without the influence of the hyper-parameter estimation. comparing the absolute, squared and 0– errors, we found that k^ generally had the lowest upper bound of the confidence interval at level 95%, for any type of error and any type of value of the parameter ρ. moreover, all estimators always had a similar confidence interval of the 0– error, while using ρ^ in most cases the upper bound of the confidence interval of the absolute and the squared error was lower than using ρ^ all these results support the suggestion to use k^ with ρ^ we should also observe that in general k^ underestimates k while k^ and k^ overestimate it. in addition, the percentage of the underestimations increases using ρ^ comparison among the boundary estimators we compared the boundary estimators on the same datasets, previously used for the estimators of the number of segments, with both the estimated and the true value of the parameters involved. the following errors were taken into account: the sum 0– error, the joint 0– error and the binary error, defined in section methods, and the average square error, aver. square error:=1k−1∑i=1k−1min⁡j=...,k− which corresponds to the mean square error over the whole vector of estimated boundaries. as observed in section methods, when the estimated segment number is used in the estimation of the boundaries, we are only able to compute the binary error, because it does not require that the vector of estimated boundaries has the same length as the vector of the true boundaries. we found that the best estimators were t^binerr, t^binerrak and t^joint, but it seemed that t^binerr, t^binerrak were slightly better than t^joint, when we estimated the parameters. moreover, between these ones, we preferred the latter, because its computation depended less on the estimation of k and thus it was more stable. in addition, in case of σ > ρ we observed that there was a great difference between the errors obtained using ρ^ and ρ^12: using the latter, the errors were significantly lower. the table shows the estimated mean binary errors ± standard deviation of the estimators of t applied to datasets with replicates for different values of σ and ρ we always used k^ ν^, σ^ and both ρ estimators as estimators of the other parameters involved. the estimators t^binerr and t^binerrak always had the lowest errors. when σ > ρ the error was lower by using ρ^ otherwise the two ρ estimators gave similar errors. to choose between the estimators t^binerrak and t^joint, we performed the segment level estimation, considering both ρ estimators. then, we computed the mse of the estimated segment level per probe and its upper bound of the confidence interval at level 95%. the results showed that, for a given estimator of ρ in general t^binerrak was better than t^joint and the upper bound of the confidence interval of the mse of the former estimator was lower than that one of the latter. moreover, using ρ^ the error was generally lower. in addition, we compared the behavior of our boundary estimators on dataset simulated chromosomes. to assess the goodness of their estimation, we measured the sensitivity and the false discovery rate, while to assess the influence of the boundary estimation on the profile estimation, we calculated the sum of squared distance, the median absolute deviation and the accuracy. the sensitivity and the fdr were computed not only looking at the exact position of the breakpoints, but also accounting for a neighborhood of or probes around the true positions. we also computed the accuracy inside and outside the aberrations separately, since the samples of dataset simulated chromosomes contained only few small copy number changes and thus the accuracy depended more on the correct estimation/classification of the probes in the "normal" regions. a more detailed explanation of these measures can be found in. since we estimated σ^ = and ρ^ =, we expected that using ρ^ we should obtain better results because σ was lower than ρ and indeed the results obtained with ρ^ were slightly better than the others. moreover, we found in general that t^binerrak had the highest sensitivity but also a higher fdr than t^. this was due to the fact that, although the estimated number of segments was the same, t^ could estimate some breakpoints with the same position, reducing the total number of breakpoints and so reducing the fdr. we can see in table that the false estimated breakpoints did not negatively influence the profile estimation. in fact, the false breakpoints are often used by the algorithm in two ways: either to divide a long segment into two or more segments with close levels or, if it is difficult to determine the position of a breakpoint, to add, before or after the aberration, a segment of one point with the value of the level between zero and the aberration level. overall, t^binerrak with ρ^ performed best on this dataset. this table shows the comparison among the error measures for profile estimation, obtained on dataset simulated chromosomes, for all estimators of t. we always used k^ ν^, σ^ and both ρ estimators as estimators of the other parameters involved. the estimator t^binerrak always had the lowest ssq and mad errors and the highest accuracy both inside and outside the regions of aberration. using ρ^ the performance was slightly better, because σ <ρ in conclusion, we suggest to use t^binerrak, even if t^joint is also a quite good estimator in some cases. regarding the estimation of ρ it seems that it is better to use ρ^ in presence of high noise. comparison among the regression curves we compared the estimation of the levels of brc with the one of brcak, also taking into account the influence of the different estimators of the parameters on the final results. to valuate the performance of the methods, we used the root mean square error per probe and per sample, computed with respect to the true profile of the levels. for this purpose, we necessarily needed datasets with replicates. using brcak we generally obtained a better or equal result with respect to the brc. moreover, we observed that, using brc, it was better to estimate the segment number with k^ or k^ note that we still have to solve the problem to determine which is the best estimator of ρ in most cases, the profile obtained by using ρ^ was better than using ρ^. this is due to the fact that sometimes ρ^ slightly underestimated ρ leading to overfitting. still we recommend to use ρ^ even if it could lead to a slight overfitting especially in the case of few segments. comparison with other methods on simulated data in this subsection we compare the original and modified versions of bpcr and brc, with other existing methods for genomic copy number estimation. error measures used in the comparison we used two different measures to examine the behavior of the different methods on the collections cases and four aberrations: the root mean square error and the roc curve. each point of the roc curve has as coordinates the false positive rate and the true positive rate for a certain threshold. the tpr is defined as the fraction of probes in the true aberrations whose estimated value is above the threshold considered, while the fpr consists in the fraction of probes outside the true aberrations whose estimated value is above the threshold. hence, the roc curve measures the accuracy of the method in the detection of the true aberrations. instead, the evaluation of the several methods on dataset simulated chromosomes was accomplished using the error measures described in, already used in the study of the boundary estimators. before showing the results, we need to remember that some methods estimate the copy numbers as a piecewise constant function, while other algorithms estimate them as a continuous curve. hence, bpcr was compared to the former group of methods, while the bayesian regression curves to the latter. since some error measures of suppose that the estimated profile is piecewise constant, we applied only the former group of methods on dataset simulated chromosomes. the piecewise constant estimation we compared the original and the modified versions of bpcr with cbs, cghseg, glad, hmm, biohmm and rendersome. for thoroughness, in the modified versions of bpcr, we used both ρ^ and ρ^ as estimators of ρ k^ as estimator of k and both t^joint and t^binerrak as estimators of the boundaries. we used ρ^ when the noise was low and otherwise ρ^ examples of estimated profiles can be found in figure in general, in presence of medium noise, the glad method performed worst, since it had a high error in the level estimation of the small peaks, while, for high noise, often both glad and rendersome failed to detect the aberrations. the cghseg method did not usually exhibit an appropriate level estimation except sometimes for segments of large width. this is due to the fact that cghseg estimates the level of a segment with the arithmetic mean of the data points in the segment and this estimator performs poorly if the segment contains few data points and the breakpoint estimation is not accurate. the cbs method, in general, performed quite well, but it was unable to detect aberrations of small width, especially when the noise was high. on the collection cases and the dataset of four aberrations with snr = the rmse plots and the roc curves of the hmm method showed that it generally estimated the profile well, but sometimes it exhibited high errors near breakpoint positions, likely because it was unable to determine the true position of the breakpoints precisely. moreover, on the dataset with snr = we recognized the true issues of the estimation with hmm. the rmse plot showed that it had a high error outside the regions of the aberrations, while inside these regions the error was always more or less the same. hence, it often failed also in the estimation of the largest aberration, the easiest one to detect. the reason of this behavior of the rmse is the following. the method estimated the true profile either with only one segment, or more often with a profile consisting of a lot of small segments, but all with the same level. since in the latter case the estimated levels were close to that one of the true aberrations, the rmse was low in the regions of the aberrations but high outside them. however, the estimated profiles were not similar to the true one. in presence of medium noise, the method biohmm was more precise than hmm in the determination of the breakpoint positions and in the level estimation, while for high noise it behaved similarly to hmm. in general we found that, when σ <ρ, the version of bpcr with t^binerrak and ρ^ generally gave the best estimation compared to the other versions of bpcr and to the other methods. in the following, we will call this modified version of bpcr, mbpcr. only on the dataset with snr = we could not choose the "best" method, because this case showed the limits of all the methods considered. the problem regarding the modified versions of bpcr was essentially the estimation of the number of segments. the roc curves of the modified versions of bpcr with ρ^ were the closest to the left and the top sides of the box, while the rmse plot showed that these methods were the best methods in the estimation of the levels inside the aberrations, but not outside them. in general, in case of very high noise, all the modified versions of bpcr well detected the aberrations, but had problems in the estimation of the profile outside them because of the poor estimation of the number of segments. in fact, k^ tends to overestimate the number of segments and this problem worsens using ρ^ in conclusion, in a situation with very high noise, using ρ^ the bpcr methods detect better the small segments, but, at the same time, the large ones are divided in small segments. on the other hand, using ρ^ smaller segments are not detected and are joined to the closest large segment. finally, the comparison performed on dataset simulated chromosomes showed that cbs and mbpcr better estimated the profiles. regarding the breakpoint error measures, we found that mbpcr had the highest sensitivity, but also a higher fdr than cbs. we have already explained in the previous subsection the possible reason of the high fdr of mbpcr and we can observe again that this fact did not influence negatively the profile estimation. the glad method showed a low sensitivity and low fdr, apart from the case regarding the exact position of the breakpoints, which means that it underestimated the segment number and the estimated breakpoints were not located exactly at their true positions. also cghseg underestimated the number of segments because of low sensitivity and fdr, while hmm had low sensitivity and high fdr when w = and vice versa in the other cases, which means that it often detected the true segment number, but it was unable to put the breakpoints at their exact position. instead, biohmm solved the issue of hmm with w = but overall had a lower sensitivity than hmm. rendersome missed several true aberrations and detected some false aberrations. the table shows the comparison of the level estimations obtained using several piecewise constant methods on dataset simulated chromosomes. in this comparison, the methods cbs and mbpcr exhibited the lowest ssq error in the profile estimation and the highest accuracy inside the aberrated regions. on the other hand, hmm, biohmm and rendersome had the highest accuracy outside the aberrations, but a high ssq error. therefore, the former group of algorithms globally estimated a better profile than the latter. because of its definition, the mad error is less informative: it does not take into account if a small number of probes are wrongly estimated, but these probes could correspond to breakpoints or small aberrations. estimation with a continuous curve we compared the several versions of the bayesian regression curves with methods which estimate the copy number as a continuous curve. lowess is the acronym of "locally weighted smoothing" and it is one of the methods considered in the comparison performed in. as we saw previously, both the brc, which uses k^ and the brcak perform well, so we tested both versions with both estimators of ρ figure shows examples of estimated profiles with these smoothing methods. in general, we found that all methods detected the regions of aberration quite well. the wavelet method showed a higher error in the level estimation of the aberrations in the datasets snr = and snr =. the methods lowess and quantreg had the highest rmse in the collection cases, while their error was not significantly different outside and inside the aberrations on datasets with snr = therefore, in the last cases the error was low inside the aberrations and high outside them in comparison with the other methods. the method smoothseg showed a similar behavior, but with a lower error. moreover, we found that the roc measure was affected by oscillations in the estimated curve, which lead to roc curves intersected and difficult to be interpreted. this complex behavior is a consequence of the way in which lowess, wavelet, quantreg and smoothseg yielded oscillating curves with positive and negative values outside the aberrations; while brcs estimated the true profile with a line almost flat and close to zero. in conclusion, the version of brc with k^ and brcak gave in general a better estimation than the other brcs and the other smoothing methods considered. regarding the ρ estimation, we found that it is better to use ρ^ if σ <ρ and ρ^ if σ > ρ application to real data in this subsection, we show how mbpcr performed compared to other piecewise constant estimation methods on the real data. we used samples from three mantle cell lymphoma cell lines previously analyzed by us with the affymetrix genechip mapping 10k array, an oligonucleotides-based microarray. we also used the data obtained on jeko- by using the higher density affymetrix genechip mapping 250k nsp array. we considered eight recurrent gene regions of aberration in lymphoma plus other two gene regions and we compared the corresponding copy numbers obtained by the several piecewise constant methods with those obtained by the fish technique in. in the end, we also show a comparison among the estimated profiles of chromosome of jeko- the 10k array data used are freely available at the public repository gene expression omnibus with geo accession: gsm gsm and gsm the 250k array data of cell line jeko- will be soon available in the same repository. with the current implementation, on a computer with dual cpu and gb ram, the algorithm needed about minutes to analyze a 10k array sample, while about day to estimate the profile of a 250k array sample. the computations were performed by chromosome and using kmax = gene copy number estimation to properly evaluate the methods, the knowledge of the true underling profile is required. in general, large aberrations on chromosomes can be detected with conventional karyotype analysis or with sky-fish and one could use this information for the evaluation procedure, but the width of these aberrations is so large that all the methods can detect them well, leading to a useless comparison. for this reason, we decided to take into account only genes to compare the piecewise constant methods. in the comparison, as previously published, when two fish copy numbers had been assigned to one gene, the first number should correspond to the copy number detected in the majority of the cells. we assigned two estimated copy numbers to one gene, when the position of the gene is between two snps and the method assigned two different values to these snps. the results on rec- did not show any significant difference among the methods, instead those on granta- showed that glad was unable to detect the true copy number in five cases, while hmm, biohmm and rendersome detected an amplification on malt greater than what detected by fish analysis. all methods did not detect the true copy number of atm, probably because the snps around atm are far away from the corresponding fish region and the deletion affects only this region. only mbpcr with ρ^ and hmm detected a breakpoint between the two snps around atm region, indicating a copy number change. regarding the jeko- data, since the cell line is triploid, to obtain more realistic copy number value, we centered the estimated log2ratio around log with the denser 250k array data, all methods behaved equally good. only hmm had a problem in the detection of the breakpoint corresponding to the c-myc amplification. on both arrays, all methods identified a gain at the ccnd position, while the copy number detected by fish is this fact cannot be explained as previously for atm, because this region is well covered by snps. instead, on the jeko- 10k array data, the noisiest among all samples, we can see several cases in which cbs, hmm and glad did not detect correctly the gene copy number. this occurred more frequently to biohmm and rendersome, while only once to cghseg. the method mbpcr with ρ^ always estimated gene copy numbers correctly, apart from ccnd on this noisy data, biohmm and rendersome often estimated the gene copy number wrongly, while this occurred only sometimes to cbs, hmm and glad. the method mbpcr with ρ^ correctly estimated the gene copy numbers, apart from ccnd whose copy number was estimated by all methods differently from the fish technique. profile estimation to compare the profile estimations, we chose the sample jeko- because, using the results obtained on both types of array, we could at least understand which regions were more realistically estimated. for now, whole validated chromosomic profiles are not available. among all chromosomes, we chose chromosome since three of the previous genes belong to that: ccnd, birc and atm. from the graphs in figure we can observe that, among all the piecewise constant methods, only mbpcr with ρ^ was able to detect the high amplification after position mb on the 10k array data, while it was recognized by all methods on the 250k array data. moreover, on the 10k array data, almost all methods detected a false deletion around position mb, due to the presence of a sequence of outliers, and biohmm did not find any copy number change in the chromosome. on the 250k array data, hmm and rendersome had problems in recognizing the last part of the chromosome as a flat region. moreover, on the 10k array data, rendersome estimated several outliers as true aberrations and, on the 250k array data, it was unable to identify the whole region from about mb to mb as gained.  we introduced new estimators for the parameters involved in bpcr and we selected the best ones on the basis of theoretic and empirical results. in particular, we found that the best way is to estimate the segment number with k^ and the boundaries with t^binerrak. we call mbpcr the bpcr version which uses k^ and t^binerrak. concerning the estimation of the variance of the segment levels, we found that the original estimator ρ^ overestimates ρ by an addendum proportional to σ, see equation. hence, the estimation fails when σ > ρ the new estimator ρ^ is more precise but slightly underestimates ρ leading to an overestimation of the segment number. applying both estimators on artificial datasets, we found that, in general, the best way is to use ρ^ when σ <ρ, but to use ρ^ when σ > ρ, even if it could lead to a slight overfitting. on real dna copy number data, commonly σ > ρ we also compared mbpcr with other methods which also estimate the copy number as a piecewise constant function: cbs, hmm, cghseg, glad, biohmm and rendersome. as a whole, the results showed that mbpcr gave the best estimation on the dataset used. however, when σ ≫ ρ it is hard to understand which method is the most appropriate. most of the other methods were not able to detect aberrations with a small width and the same was true for mbpcr using ρ^ on the other hand, the use of ρ^ led to the detection of the smaller segments, but the larger ones were often divided in small segments and sometimes the segments consisted of only one point. the optimal choice of the ρ estimator is still not fully determined. the new estimators improved also brc, which is a bayesian regression with a smooth curve. moreover, we derived a formula to estimate brc without employing the estimated number of segments. we referred to it as brcak. applying these methods to artificial data, the best estimators were found to be the brc version with k^ and brcak. about the choice of the two estimators of ρ we found a similar conclusion as before, with the advantage of being less problematic when σ ≫ ρ we compared these two regression methods with other methods which estimate the copy number data as a continuous curve: wavelet, lowess, quantreg and smoothseg. the results showed that our modified regression methods were the most appropriate for the estimation of the segment levels on the datasets considered. even though these smoothing methods seem to have less problems in the estimation and the error measures suggest that their estimation is even better than the piecewise constant estimation, these methods do not detect the position of the breakpoints explicitly and hence the changes in the value of the segment level. thus, they seem to be less adequate in practice. for this reason, we feel that the roc curve cannot be used as the only measure to compare methods, as previously done in. the rmse is generally an acceptable measure, but we observed that in some cases even this is not sufficient, due to the overfitting. willenbrock and fridlyand proposed other measures to compare methods for copy number estimation, regarding both breakpoint and level estimation. in particular, the sensitivity measure of breakpoint estimation is useful to select which methods should be used, because it quantifies the precision of the methods in determining the position of the breakpoints. finally, we have applied mbpcr and all other piecewise constant regression methods to real data. the comparisons showed that mbpcr estimated well the copy number of the genes. on these data, in most cases the choice of the ρ estimator did not affect the analysis. in comparison with the other methods, the current implementation of our algorithm is computationally intensive. the real computational time can be reduced linearly diminishing kmax and quadratically diminishing the number of probes. moreover, the computation can be easily parallelized by arm and by chromosome, reducing further the calculation time. in cancer research, the accuracy in the dna copy number estimation is crucial for the correct determination of the mutations that characterize the disease. in particular, the estimation of the breakpoints must be precise to detect correctly which genes are affected by these genomic aberrations. as recently shown, snp microarrays can also potentially detect the breakpoints involved in unbalanced translocations, allowing the identification of fusion genes. in this context, the use of our method can highly improve the disease investigation, because it accurately determines breakpoints, is less sensitive to high noise and generally outperforms all the methods considered in our comparisons. moreover, smoothing algorithms are clearly not suitable for such analysis. availability and requirements project name: mbpcr. project home page:. operating systems: the software should run in linux, mac-os or windows. tests were performed on windows and linux systems. programming language: r. other requirements: none. licence: gnu gpl. any restrictions to use by non-academics: none. authors' contributions pmvr carried out the study and wrote the manuscript. mh and ik supervised the statistical analysis. fb supervised the validation study and provided the real data. all authors read and approved the final manuscript. supplementary material additional file 1 mbpcr source code. this zipped file contains the source code of the mbpcr algorithm in r, including help files, sample data and examples. click here for file additional file 2 supplementary material. this file contains a brief description of the dynamic program used to implement the method, a formal definition of some error measure used, some supplementary tables and some supplementary figures. click here for file acknowledgements this work was supported by the grant 205321- of the swiss national science foundation and by oncosuisse grant ocs 1939-08- 
some diseases, like tumors, can be related to chromosomal aberrations, leading to changes of dna copy number. we propose an improved method for dna copy number estimation, mbpcr, which performed very well compared to previously published algorithms. in particular, mbpcr was more powerful in the detection of the true position of the breakpoints and of small aberrations in very noisy data. we also proposed an alternative estimator of the variance of the segment levels, which is useful in case of data with high noise. the copy number of an aberrant genome can be represented as a piecewise constant function, since it can exhibit regions of deletions or gains. instead, in a healthy cell the copy number is two because we inherit one copy of each chromosome from each our parents. for example, the boundary estimator did not take into account the dependency among the boundaries and succeeded in estimating more than one breakpoint at the same position, losing segments. we derived an improved version of the bpcr and brc, changing the segment number estimator and the boundary estimator to enhance the fitting procedure. using artificial data, we compared the original and the modified version of bpcr and brc with other regression methods, showing that our improved version of bpcr generally outperformed all the others. similar results were also observed on real data. the method estimates the unknown segment number, the endpoints of the segments and the value of the segment levels of the underlying piecewise constant function. hence, from a biological point of view, our method can be very useful, for example, to find targets of genomic aberrations in clinical cancer samples. bayesian piecewise constant regression is a bayesian regression method for data that are noisy observations of a piecewise constant function.
 some diseases, like tumors, can be related to chromosomal aberrations, leading to changes of dna copy number. the copy number of an aberrant genome can be represented as a piecewise constant function, since it can exhibit regions of deletions or gains. instead, in a healthy cell the copy number is two because we inherit one copy of each chromosome from each our parents. bayesian piecewise constant regression is a bayesian regression method for data that are noisy observations of a piecewise constant function. the method estimates the unknown segment number, the endpoints of the segments and the value of the segment levels of the underlying piecewise constant function. the bayesian regression curve estimates the same data with a smoothing curve. however, in the original formulation, some estimators failed to properly determine the corresponding parameters. for example, the boundary estimator did not take into account the dependency among the boundaries and succeeded in estimating more than one breakpoint at the same position, losing segments.  we derived an improved version of the bpcr and brc, changing the segment number estimator and the boundary estimator to enhance the fitting procedure. we also proposed an alternative estimator of the variance of the segment levels, which is useful in case of data with high noise. using artificial data, we compared the original and the modified version of bpcr and brc with other regression methods, showing that our improved version of bpcr generally outperformed all the others. similar results were also observed on real data.  we propose an improved method for dna copy number estimation, mbpcr, which performed very well compared to previously published algorithms. in particular, mbpcr was more powerful in the detection of the true position of the breakpoints and of small aberrations in very noisy data. hence, from a biological point of view, our method can be very useful, for example, to find targets of genomic aberrations in clinical cancer samples. 