Enhancing Estimates of Breakpoints in Genome Copy Number Alteration using Confidence Masks

Abstract

Chromosomal structural changes in human body known as copy number alteration (CNA) are often associated with diseases, such as various forms of cancer. Therefore, accurate estimation of breakpoints of the CNAs is important to understand the genetic basis of many diseases. The high‐resolution comparative genomic hybridization (HR‐CGH) and single‐nucleotide polymorphism (SNP) technologies enable cost‐efficient and high‐throughput CNA detection. However, probing provided using these profiles gives data highly contaminated by intensive Gaussian noise having white properties. We observe the probabilistic properties of CNA in HR‐CGH and SNP measurements and show that jitter in the breakpoints can statistically be described with either the discrete skew Laplace distribution when the segmental signal‐to‐noise ratio (SNR) exceeds unity or modified Bessel function‐based approximation when SNR is <1. Based upon these approaches, the confidence masks can be developed and used to enhance the estimates of the CNAs for the given confidence probability by removing some unlikely existing breakpoints.

1. Introduction

It is well known that the deoxyribonucleic acid (DNA) of a genome essential for human life often demonstrates structural changes [1–3] called genome copy number alterations (CNAs) [4–6], which are associated with disease such as cancer [7]. Analysis of the breakpoint locations in the CAN structure is still an important issue because it helps detecting structural alterations, load of alterations in the tumor genome, and absolute segment copy numbers. Thus, efficient estimators are required to extract information about the breakpoints with accuracy acceptable for medical needs. To produce CNA profile, several technologies have been developed such as comparative genomic hybridization (CGH) [8], high‐resolution CGH (HR‐CGH) [9], whole genome sequencing [10], and most recently single‐nucleotide polymorphism (SNP) [11]. The HR‐CGH technology is still used widely in spite of its low resolution [12]. It has been reported in [13] that the HR‐CGH arrays are accurate to detect structural variations (SVs) at the resolution of 200 bp (base pairs). Most recently, the single‐nucleotide polymorphism technology was developed in the study of Wang et al. [11] to provide high‐resolution measurements of the CNAs. In spite of their high resolution, the modern methods still demonstrate the inability in obtaining good estimates of the breakpoint locations because of the following factors: (1) the nature of biological material (tumor is contaminated by normal tissue, relative values, and unknown baseline for copy number estimation), (2) technological biases (quality of material and hybridization/sequencing), and (3) intensive random noise. The HR‐CGH and SNP profiles have demonstrated deficiency in detecting the CNAs, but noise in the detected changes still remains at a high level [14] and accurate estimators are required to extract information about structural changes.

In the HR‐CGH microarray technique, the CNAs are often normalized and plotted as ${\log }_{2}R/G={\log }_2 \text{ratio}$, where R and G are the fluorescent Red and Green intensities, respectively [12]. The CNA measurements using SNP technologies are represented by the Log‐R ratios (LRRs), which are the log‐transformed ratios of experimental and normal reference SNP intensities centered at zero for each sample [14]. From the standpoint of signal processing, the following properties of the CNA function are of importance [15]:

• It is piecewise constant (PWC) and sparse with a small number of alterations on a long base‐pair length.

• Constant values are integer, although this property is not survived in the log‐R ratio.

• The measurement noise in the log‐R ratio is highly intensive and can be modeled as additive white Gaussian.

The CNA estimation problem is thus to predict the breakpoint locations and the segmental levels with a maximum possible accuracy and precision acceptable for medical applications. In this work, we developed our methods to two types of cancer: B‐cell chronic lymphocytic leukemia (B‐CLL) and BLC primary breast carcinoma. Nevertheless, the methods were designed to any samples of cancer with the characteristics described above.

2. Methods

2.1. CNA model and problem statement

Consider a chromosome section observed with some resolution $\Delta$, bp at $M$ discrete breakpoints, $n\in \left[1,M\right]$. An example of the CNA probes with a single breakpoint and two segments is shown in Figure 1. Suppose that the copy numbers change at $K$ breakpoints, $1<i_1<\dots <i_K<M$, united in a vector $I={[i_1{i}_2\dots i_K]}^{\text{T}}$. The measurement can thus be represented with a vector $y\in {ℛ}^M$ as

$y={\left[y_1{y}_2\dots y_{i_1}{y}_{i_1}+1\dots y_{i_2}\dots y_{i_K}\dots y_M\right]}^T.$E1

Figure 1.

Typical CNA measurements with white Gaussian noise with a single breakpoint, between two segments $l$ and $l+1$ having different segmental variances. The pdf for neighboring segments are depicted as $p_l(x)$ and $p_{l+1}(x)$.

Introduce a vector $a\in {ℛ}^{K+1}$ of segmental levels,$a={[a_1{a}_2\dots a_{K+1}]}^{\text{T}}$ , where $a_1$ corresponds to the interval $[1,i_1],a_{K+1}$ to $\left[i_K,M\right]$, and $a_K$, $k\ge 2$, to $[i_{k-1},i_k]$. In such a formulation, $y$ obeys the linear regression model

$y=A\left(I\right)a+\nu$E2

where the regression matrix $A\in {ℛ}^{M\times (K+1)}$ is sparse,

$A={\left[A_1^T{A}_2^T\dots A_{K+1}^T\right]}^T$E3

having a component

$A_k=\left[\begin{array}{ccccc}0& \cdots & 1& \cdots & 0\\ 0& \cdots & 1& \cdots & 0\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& \cdots & 1& \cdots & 0\end{array}\right]$E4

in which the $k$th column is filled with unity and all others are zeros. The number of the columns in $A_k$ is exactly $K+1$. However, the number or the row depends on the interval $i_k-i_{k-1}$. Thus, the row‐variant matrix (4) is $A_k\in {ℛ}^{\left(i_k-i_{k-1}\right)\times \left(K+1\right)}$. Additive noise $v$ in Eq. (2) is zero mean, $E\left\{v\right\}=0$, and white Gaussian with the covariance $R={\sigma }_v^{2}I$, where $I\in {ℛ}^{M\times M}$ is an identity matrix and ${\sigma }_v^2$ is a know variance.

The CNA estimation problem is thus to predict the breakpoint locations and evaluate the segmental changes $x=A\left(I\right)a$ with a maximum possible accuracy and precision acceptable for medical applications. The problem is complicated by short number of the probes in each neighboring segment and indistinct edges. Therefore, an analysis of the estimation errors caused by the segmental noise and jitter in the breakpoints is required.

2.2. Jitter probability in the breakpoints

Consider a typical genomic measurement of two neighboring CNA segments in white Gaussian noise with different segmental variances as shown in Figure 1. A constant signal changes from level $a_l$ to level $a_{l+1}$ around the breakpoint$i_l$. In the presence of noise, the location of $i_l$ is not clear owing to commonly large segmental variances ${\sigma }_l^2$ and ${\sigma }_{l+1}^2$. As an example, the Gaussian noise probability density functions (pdfs) $p_l(x)$ and $p_{l+1}(x)$ are shown in Figure 1 for ${\sigma }_l^2>{\sigma }_{l+1}^2$. Let us notice that $p_l(x)$ and $p_{l+1}(x)$ cross each other in two points, ${\alpha }_l$ and ${\beta }_l$, provided that ${\sigma }_l^2\ne {\sigma }_{l+1}^2$.

Now considerer $N$ probes in each segment neighboring to $i_l$ with an average resolution. We thus may assign an event $A_{lj}$ meaning that measurement at point $i_l-N\le j<i_l$ belongs to the $l$th segment. Another event $B_{lj}$ means that measurement at $i_l\le j<i_l+N-1$ belongs to the $(l+1)$th segment. In our approach, we think that a measured value belongs to one segment if the probability is larger than if it belongs to another segment. For example, any measurement point in the interval between ${\alpha }_l$ and ${\beta }_l$ (Figure 1) is supposed to belong to the $(l+1)$th segment. Following Figure 1 and assuming different noise variances ${\sigma }_l^2$ and ${\sigma }_{l+1}^2$, the events $A_{lj}$ and $B_{lj}$ can be specified as follows [16]:

$\begin{array}{c}{A}_{lj}\\ i_l-N_l\le j<i_l\end{array}\text{is }\{\begin{array}{cc}({\alpha }_l<x_j)\vee \left(x_j<{\beta }_l\right),& {\sigma }_l^2>{\sigma }_{l+1}^2,\\ x_j>{\alpha }_l,& {\sigma }_l^2={\sigma }_{l+1}^2,\\ {\alpha }_l<x_j<{\beta }_l,& {\sigma }_l^2<{\sigma }_{l+1}^2,\end{array}$E5

$\begin{array}{c}{B}_{lj}\\ i_l\le j<i_l+N_l-1\end{array}\text{is }\{\begin{array}{cc}{\beta }_l<x_j<{\alpha }_l& {\sigma }_l^2<{\sigma }_{l+1}^2,\\ x_j<{\alpha }_l,& {\sigma }_l^2={\sigma }_{l+1}^2,\\ (x_j<{\alpha }_l)\vee \left(x_j>{\beta }_l\right),& {\sigma }_l^2>{\sigma }_{l+1}^2.\end{array}$E6

Because each point can belong only to one segment, the inverse events are $\overline{A}=1-A_{lj}$ and $\overline{B}=1-B_{lj}$.

Events $A_{lj}$ and $B_{lj}$ can be united into two blocks: $A_l=\left\{A_{l\left(i_l-N\right)}{A}_{l\left(i_l-N+1\right)}\dots A_{l\left(i_l-1\right)}\right\}\text{ and }{B}_l=\left\{B_{l\left(i_l\right)}{B}_{l\left(i_l+1\right)}\dots B_{l\left(i_l+N-1\right)}\right\}$.

If $A_l$ and $B_l$ occur simultaneously with unit probability each, then jitter at $i_l$ will never occur. However, some other events may be found, which do not obligatorily lead to jitter. We ignore such events and define approximately the probability $P(A_l,B_l)$ of the jitter‐free breakpoint as

$P\left(A_l,B_l\right)=P\left({\text{A}}_{{\text{i}}_{\text{l}}-N}\dots {\text{A}}_{{\text{i}}_{\text{l}}-1}{\text{B}}_{{\text{i}}_{\text{l}}}\dots {\text{B}}_{{\text{i}}_{\text{l}}+N-1}\right)$E7

The inverse event $\overline{P}\left(A_l,B_l\right)=1-P(A_l,B_l)$ can be called the approximate jitter probability [17].

2.3. Jitter distribution in the breakpoints

To determine the confidence limits for CNAs using high‐resolution genomic arrays, jitter in the breakpoints must be specified statistically for the segmental Gaussian distribution. This can be done approximately if to employ either the discrete skew Laplace distribution or, more accurately, the modified Bessel function of the second kind and zeroth order.

2.3.1. Approximation with discrete skew Laplace distribution

Following the definition of the jitter probability given in Section 2.1 and taking into consideration that all the events are independent in white Gaussian noise, Eq. (7) can be rewritten as: $P\left(A_l,B_l\right)=P^N\left(A_l\right)P^N\left(B_l\right),$ where, following Eqs. (5) and (6), the probabilities $P(A_l)$ and $P(B_l)$ can be specified as, respectively,

$P\left(A_l\right)=\{\begin{array}{cc}1-\underset{{\beta }_l}{\overset{{\alpha }_l}{{\displaystyle \int }}}{p}_l\left(x\right)dx,& {\sigma }_l^2>{\sigma }_{l+1}^2,\\ \underset{{\alpha }_l}{\overset{\infty }{{\displaystyle \int }}}{p}_l\left(x\right)dx,& {\sigma }_l^2={\sigma }_{l+1}^2,\\ \underset{{\alpha }_l}{\overset{{\beta }_l}{{\displaystyle \int }}}{p}_l\left(x\right)dx,& {\sigma }_l^2<{\sigma }_{l+1}^2,\end{array}$E8

$P\left(B_l\right)=\{\begin{array}{cc}\underset{{\beta }_l}{\overset{{\alpha }_l}{{\displaystyle \int }}}{p}_{l+1}\left(x\right)dx,& {\sigma }_l^2>{\sigma }_{l+1}^2,\\ \underset{-\infty }{\overset{{\alpha }_l}{{\displaystyle \int }}}{p}_{l+1}\left(x\right)dx,& {\sigma }_l^2={\sigma }_{l+1}^2,\\ 1-\underset{{\alpha }_l}{\overset{{\beta }_l}{{\displaystyle \int }}}{p}_{l+1}\left(x\right)dx,& {\sigma }_l^2<{\sigma }_{l+1}^2,\end{array}$E9

where $p_l\left(x\right)=1/\sqrt{2\pi {\sigma }_l^2}{e}^{-({\left(x-a_l\right)}^2)/{\sigma }_l^2}$ is the Gaussian density.

Suppose that jitter occurs at some point $i_l\pm k,0\le k\le N,$ as shown, for example, in Figure 1, and assign two additional blocks of events $A_{lk}=\{A_{i_l-N}\dots A_{i_l-1-k}\}$ and $B_{lk}=\{B_{i_l+k}\dots B_{i_l+N-1}\}$. The probability $P_k^{-}\triangleq P_k^{-}\left(A_{lk}{\overline{A}}_{l\left(i_l-k\right)}\dots {\overline{A}}_{i_l-1}{B}_l\right)$ that jitter occurs at the $k$th point to the left from $i_l$ (left jitter) and the probability $P_k^{+}\triangleq P_k^{+}\left(A_l{\overline{B}}_{l\left(i_l+1\right)}\dots {\overline{B}}_{l\left(i_l+k-1\right)}{B}_k\right)$ that jitter occurs at the $k$th point to the right from $i_l$ (right jitter) can thus be written as, respectively,

$P_k^{-}=P^{N-k}\left(A_l\right){\left[1-P\left(A_l\right)\right]}^k{P}^N\left(B_l\right),$E10

$P_k^{+}=P^N\left(A_l\right){\left[1-P\left(B_l\right)\right]}^k{P}^{N-k}\left(B_l\right),$E11

By normalizing Eqs. (11) and (12) with Eq. (8), one can arrive at a function that turns out to be independent on $N$:

$f_l\left(k\right)=\{\begin{array}{cc}{[P^{-1}\left(A_l\right)-1]}^{\left|k\right|},& k<0,(\text{left})\\ 1& k=0\\ {[P^{-1}\left(B_l-1\right)]}^k,& k<0.(\text{right})\end{array}.$E12

Further normalization of $f_l(k)$ to have a unit area leads to the pdf $p_l\left(k\right)=\frac{1}{{\phi }_l}{f}_l(k)$, where ${\phi }_l$ is the sum of $f_l(k)$ for all $k$,

${\varphi }_l=1+{\displaystyle \sum }_{k=1}^{\infty }[{\phi }_l^A\left(k\right)+{\phi }_l^B\left(k\right)],$E13

where ${\phi }_l^A\left(k\right)={[P^{-1}\left(A_l\right)-1]}^k$ and ${\phi }_l^B\left(k\right)={[P^{-1}\left(B_l\right)-1]}^k$.

It follows from the approximation admitted that $f_l(k)$ converges with $k$ only if $0.5<\tilde{P}=\left\{P\left(A\right),P\left(B\right)\right\}<1$. Otherwise, if $\tilde{P}<0.5$, the sum ${\phi }_l$ is infinite and $f_l(k)$ cannot be transformed to $p_l(k)$. It has been shown in [18] that such a situation is practically rare. It can be observed with extremely small and different segmental SNRs when the probabilities are comparable that the measurement point belongs to one of another segment.

Accepting $0.5<\tilde{P}=\left\{P\left(A\right),P\left(B\right)\right\}<1$, one concludes that $\tilde{P}<0$, $\ln \left(1-\tilde{P}\right)<0$, and $\ln \left(1-\tilde{P}\right)<\text{ln}(\tilde{P})$. Next, using a standard relation ${\displaystyle \sum }_{k=1}^{\infty }{x}^k=1/(x^{-1}-1)$, where $x<1$, and after little transformations, Eq. (14) can be brought to

${\phi }_l=\frac{P\left(A_l\right)+P\left(B_l\right)-1}{[1-2P\left(A_l\right)\left[1-2P\left(B_l\right)\right]}$E14

The jitter pdf $p_l(k)$ associated with the lth breakpoint can finally be found to be

$p_l\left(k\right)=\frac{1}{{\phi }_l}\{\begin{array}{cc}{[P^{-1}\left(A_l\right)-1]}^{\left|k\right|},& k<0,\\ 1,& k=0,\\ {[P^{-1}\left(B_l\right)-1]}^k& k>0,\end{array}$E15

where ${\varphi }_l$ is specified by Eq. (15) and $0.5<P\left(A_l\right),P\left(B_l\right)<1$.

If now to substitute $q_l=P^{-1}\left(A_l\right)-1$ and $d_l=P^{-1}\left(B_l\right)-1$, find $P\left(A_l\right)=1/(1+q_l)$ and $P\left(B_l\right)=1/(1+d_l)$, and provide the transformations, then one may arrive at a conclusion that Eq. (16) is the discrete skew Laplace pdf [19].

$p\left(k\text{|}{d}_l{q}_l\right)=\frac{(1-d_l)(1-q_l)}{1-d_l{q}_l}\{\begin{array}{cc}{p}_l^k,& k\ge 0\\ q_l^{\left|k\right|}& k\le 0\end{array}$E16

where $d_l=e^{-({\kappa }_l/{\nu }_l)}\in [0,1]$ and $q_l=e^{-(1/{\kappa }_l{\nu }_l)}\in [0,1]$ and in which ${\kappa }_l$ and ${\nu }_l>0$ still need to be connected to Eq. (16). With this aim, consider Eqs. (16) and (17) at $k=-1,k=0$, and $k=1$. By equating Eqs. (16) and (17), first obtain $(\left(1-d_l\right)\left(1-q_l\right)d_l)/(1-d_l{q}_l)=1/{\varphi }_l(1-P(B_l))/P(B_l)$ for $k=1$ and $(\left(1-d_l\right)\left(1-q_l\right)q_l)/(1-d_l{q}_l)=1/{\varphi }_l(1-P(A_l))/P(A_l)$ for $k=-1$ that yields

${\nu }_l=\frac{1-{\kappa }_l^2}{{\kappa }_l\ln {\mu }_l},$E17

where ${\mu }_l=(P\left(A_l\right)[1-P(B_l)])/(P(B_l)[1-P(A_l)])$. For $k=0$ we have $(\left(1-d_l\right)\left(1-q_l\right))/(1-d_l{q}_l)=1/{\varphi }_l$ and transform it to the equation $x_l^2-({\varphi }_l(1+{\mu }_l))/(1+{\varphi }_l)x-(1-{\varphi }_l)/\left(1+{\varphi }_l\right){\mu }_l=0$, where a proper solution is

$x=\frac{{\phi }_l(1+{\mu }_l)}{2(1+{\phi }_l)}\left(1-\sqrt{1-\frac{4{\mu }_l(1-{\phi }_l^2)}{{\phi }_l^2{\left(1+{\mu }_l\right)}^2}}\right)$E18

and which $x_l={\mu }_l^{-\left({\kappa }_l^2\right)/(1-{\kappa }_l^2)}$ gives us

${\kappa }_l=\sqrt{\ln \frac{\left(x_l\right)}{\ln (x_l/{\mu }_l)}}.$E19

By combining Eq. (18) with Eq. (20), one may also get a simpler form for ${\nu }_l$, namely ${\nu }_l=-{\kappa }_l/\text{ln}{x}_l$.

Now, introduce the segmental signal‐to‐noise ratios (SNRs): ${\gamma }_l^{-}=\frac{{\Delta }_l^2}{{\sigma }_l^2}$, and ${\gamma }_l^{+}=\frac{{\Delta }_l^2}{{\sigma }_{l+1}^2}$ , where ${\Delta }_l=a_{l+1}-a_l$, substitute the Gaussian pdf to Eqs. (9) and (10), provide the transformations, and rewrite Eqs. (9) and (10) as

$P\left(A_l\right)=\{\begin{array}{cc}1+\frac{1}{2}\left[\text{erf}\left(g_l^{\beta }\right)-\text{erf}\left(g_l^{\alpha }\right)\right],& {\gamma }_l^{-}<{\gamma }_l^{+}\\ \frac{1}{2}\text{erf}c\left(g_l^{\alpha }\right),& {\gamma }_l^{-}={\gamma }_l^{+}\\ \frac{1}{2}\left[\text{erf}\left(g_l^{\beta }\right)-\text{erf}\left(g_l^{\alpha }\right)\right],& {\gamma }_l^{-}>{\gamma }_l^{+}\end{array},$E20

$P\left(B_l\right)=\{\begin{array}{cc}\frac{1}{2}\left[\text{erf}\left(h_l^{\alpha }\right)-\text{erf}\left(h_l^{\beta }\right)\right],& {\gamma }_l^{-}<{\gamma }_l^{+}\\ 1-\frac{1}{2}\text{erfc}\left(h_l^{\alpha }\right),& {\gamma }_l^{-}={\gamma }_l^{+}\\ 1+\frac{1}{2}\left[\text{erf}\left(h_l^{\alpha }\right)-\text{erf}\left(h_l^{\beta }\right)\right],& {\gamma }_l^{-}>{\gamma }_l^{+}\end{array},$E21

where $g_l^{\beta }=({\beta }_l-{\Delta }_l)/\left|{\Delta }_l\right|\sqrt{{\gamma }_l^{-}/2}$, $g_l^{\alpha }=({\alpha }_l-{\Delta }_l)/\left|{\Delta }_l\right|\sqrt{{\gamma }_l^{-}/2}$, $h_l^{\beta }={\beta }_l/\left|{\Delta }_l\right|\sqrt{{\gamma }_l^{+}/2}$, $h_l^{\alpha }={\alpha }_l/\left|{\Delta }_l\right|\sqrt{{\gamma }_l^{+}/2}$, $\text{erf}\left(x\right)$ is the error function, $\text{erfc}\left(x\right)$ is the complementary error function, and

${\alpha }_l,{\beta }_l=\frac{a_l{\gamma }_l^{-}-a_l{\gamma }_l^{+}}{{\gamma }_l^{-}-{\gamma }_l^{+}}\mp \frac{1}{{\gamma }_l^{-}-{\gamma }_l^{+}}\times \sqrt{\left(a_l-a_{l+1}\right){\gamma }_l^{-}{\gamma }_l^{+}+2{\Delta }_l^2\left({\gamma }_l^{-}-{\gamma }_l^{+}\right)\ln \sqrt{\frac{{\gamma }_l^{-}}{{\gamma }_l^{+}}}}.$E22

2.3.2. Approximation of jitter distribution using the modified Bessel functions

An analysis shows that the discrete skew Laplace pdf (17) gives good results only if SNR is >1. Otherwise, real measurements do not fit well, and a more accurate function is required. Below, we show that better approach to real jitter distribution can be provided using the modified Bessel functions.

2.3.2.1. Modified Bessel function

Figure 2 demonstrates the jitter pdf measured experimentally (dotted) for different SNRs. The breakpoint corresponds here to the peak density and the probability of the breakpoint location diminishes to the left and to the right of this point. Note that the discrete skew Laplace pdf (17) behaves linearly in such scales. Therefore, Eq. (17) cannot be applied when SNR is <1 and a more accurate function is required.

Figure 2.

Experimentally defined one‐sided jitter probability densities (dotted) of the breakpoint location for equal segmental SNR $\gamma$ in the range of $M=400$ points with a true breakpoint at $n=200$. The experimental density functions were found using the Maximum Likelihood (ML) estimator. The histogram was plotted over $50\times {10}^3$ runs repeated nine times and average. Approximations (continuous) are provided using the proposed Bessel‐based approximation depicted as MBA.

Among available functions demonstrating the pdf properties, the modified Bessel function of the second kind $K_0\left(x\right)$ and zeroth order is a most good candidate to fit the experimentally measured densities (Figure 2). The following form of $K_0(x)$ can be used:

$\begin{array}{c}{K}_0\left[x\left(k\right)\right]=\underset{0}{\overset{\infty }{{\displaystyle \int }}}\cos \left[x\left(k\right)\sinh t\right]dt\\ =\underset{0}{\overset{\infty }{{\displaystyle \int }}}\frac{\cos \left[x\left(k\right)t\right]}{\sqrt{t^2+1}}dt>0,x\left(k\right)>0,\end{array}$E23

in which a variable $x\left(k\right)$ depends on index $k,$ which represents a discrete departure from the assumed breakpoint location. Because $K_0\left[x\left(k\right)\right]$ is a positive‐valued for $x\left(k\right)>0$ smooth function decreasing with $x$ to zero, it can be used to approximate the probability density.

2.3.2.2. Approximation

In order to use Eq. (24) as an approximating function


$B\left(k\text{|}\gamma \right)=K_0\left[x\left(k\right)\right]$E24

Figure 3.

Simulated CNA with a single breakpoint at n = 200 and segmental standard deviations ${\sigma }_l$ and ${\sigma }_{l+1}$ corresponding to SNRs ${\gamma }_l^{-}={\gamma }_l^{+}=1.37$: (a) measurement and (b) jitter distribution. Here, ML (circled) is the jitter pdf obtained experimentally using an ML estimator through a histogram over $50\times {10}^3$ runs, SkL (solid) is the Laplace distribution, and MBA (dashed) is the Bessel‐based approximation.

conditioned on $\gamma$ for the one‐sided jitter probability densities shown in Figure 2, we represent a variable $x$ via $k$ as $x\left(k,\gamma \right)=\ln [\Phi (k,\gamma )]$ in a way such that small $k\ge 0$ correspond to large values $x$ of and vice versa. Among several candidates, it has been found empirically that the following function $\Phi (k,\gamma )$ fits the histograms with highest accuracy:

$\Phi \left(k,\gamma \right)={\left(\left|k\right|+1\right)}^{\rho +\tau \left|k\right|}\left[\frac{1+\sqrt{\gamma }}{\gamma }-ϵ\right]$E25

if to set $\gamma ={\gamma }_l^{-}$ for $k<0$, $\gamma =\frac{{\gamma }_l^{-}+{\gamma }_l^{+}}{2}$ for $k=0$, and $\gamma ={\gamma }_l^{+}$ for $k>0$, and represent the coefficients and as $\tau (\gamma ),\rho (\gamma )$, and $ϵ(\gamma )$ as

$\tau \left(\gamma \right)=a_0\gamma +a_1$E26

$\rho \left(\gamma \right)=\gamma \left(b_0{\gamma }^{b_1}+a_0\right)+b_2$E27

$ϵ\left(\gamma \right)=c_0{\gamma }^{c_1}+c_2$E28

where $a_0=0.02737$, $a_1=-4.5\times {10}^{-3}$, $b_0=0.3674$, $b_1=-0.3137$, $b_2=0.8066$, $c0=0.8865$, $c_1=-1.033,$ and $c_2=-1.233$ were found in the mean square error (MSE) sense. These values were found in several iterations until the MSE reached a minimum.

In summary, Figure 3 gives a typical example of a simulated CNA, where the modified Bessel function‐based approximation (depicted as MBA) demonstrates better accuracy than the approximation obtained using the skew Laplace distribution (depicted as SkL).

2.4. Probabilistic masks

It follows from Figure 3 that, in view of large noise, estimates of the CNAs may have low confidence, especially with small SNR $\gamma \le 1$. Thus, each estimate requires confidence boundaries within which it may exist with a given probability [20, 21].

Given an estimate ${\widehat{a}}_l$ of the $l$th segmental level in white Gaussian noise, the probabilistic upper boundary (UB) and lower boundary (LB) can be specified for the given confidence probability $P(\vartheta )$ in the $\vartheta$‐sigma sense as [20]

${\widehat{a}}_l^{UB}\cong {\widehat{a}}_l+\epsilon ={\widehat{a}}_l+\vartheta \sqrt{\frac{{\sigma }_j^2}{N_l}}={\widehat{a}}_l+\vartheta {\widehat{\sigma }}_l$E29

${\widehat{a}}_l^{LB}\cong {\widehat{a}}_l-\epsilon ={\widehat{a}}_l+\vartheta \sqrt{\frac{{\sigma }_j^2}{N_l}}={\widehat{a}}_l+\vartheta {\widehat{\sigma }}_l$E30

where $\vartheta$ indicates the boundary wideness in terms of the segmental noise variance ${\widehat{\sigma }}_l$ on an interval $N_l$ points, from ${\widehat{n}}_{l-1}$ to ${\widehat{n}}_l-1$.

Likewise, detected the $l$th breakpoint location ${\widehat{n}}_l$, the jitter probabilistic left boundary $J_l^L$ and right boundary $J_l^R$ can be defined, following [20], as

$J_l^L\cong {\widehat{n}}_l-k_l^R,$E31

$J_l^R\cong {\widehat{n}}_l+k_l^L,$E32

where $k_l^R(\vartheta )$ and $k_l^L(\vartheta )$ are specified by the jitter distribution in the $\vartheta$‐sigma sense.

By combining Eqs. (30) and (31) with Eqs. (32) and (33), the probabilistic masks can be formed as shown in [20] to bound the CNA estimates in the $\vartheta$‐sigma sense for the given confidence probability $P(\vartheta )$. An important property of these masks is that they can be used not only to bound the estimates and show their possible locations on a probabilistic field [20, 21] but also to remove supposedly wrong breakpoints. Such situations occur each time when the masks reveal double UB and LB uniformities in a gap of three neighboring detected breakpoints. If so, then the unlikely existing intermediate breakpoint ought to be removed.

Noticing that the segmental boundaries (30) and (31) remain the same irrespective of the jitter in the breakpoints, below we specify the masks for the jitter represented with the Laplace distribution (17) and Bessel‐based approximation (25).

2.4.1. Masks for Laplace distribution

For the Laplace distribution (17), the jitter left boundary $J_l^L$ (32) and right boundary $J_l^R$(33) can be defined in the $\vartheta$–‐sigma sense if to specify $k_l^R(\vartheta )$ and $k_l^L(\vartheta )$ as shown in [18],

$k_l^R=\frac{\nu }{\kappa }\ln \frac{\left(1-d_l\right)\left(1-q_l\right)}{\xi \left(1-d_l{q}_l\right)},$E33

$k_l^L=\nu \kappa \ln \frac{\left(1-d_l\right)\left(1-q_l\right)}{\xi \left(1-d_l{q}_l\right)},$E34

where [$x$] means a maximum integer lower than or equal to $x$. Note that functions (34) and (35) were obtained in [18] by equating (17) to $\xi \left(N_l\right)=\text{erfc}(\vartheta /\sqrt{2})$ and solving for $k_l$.

The probabilistic UB mask

and LB mask

for the Laplace distribution were formed in [17,20,21] by the segmental upper boundary ${\widehat{a}}_l^{\mathrm{UB}}$ and lower boundary ${\widehat{a}}_l^{\mathrm{LB}}$ and by the jitter left boundary $J_l^{\text{L}}$ and jitter right boundary $J_l^{\text{R}}$. The algorithm for computing

and

masks has been developed and applied to the CNA probes in [22].

2.4.2. Masks for Bessel‐based approximation

The UB mask

and LB mask

for the Bessel‐based approximation can be formed using the same equations as for the Laplace distribution. Suppose that the Laplace pdf (17) is equal to the approximating function ${{\rm B}}_l(k)$ at $k=0$,

$p\left(k=0\text{|}{d}_l,q_l\right)={{\rm B}}_l\left(k=0\right),$E35

that yields ${{\rm B}}_l\left(k=0\right)=\frac{1}{{\varphi }_l}$. Then, define the probabilities $P^B(A_l)$ at $k=-1$ and $P^{{\rm B}}(B_l)$ at $k=1$ as

$P^{{\rm B}}\left(A_l\right)=\frac{B_l\left(k=0\right)}{{{\rm B}}_l\left(k=-1\right)+{{\rm B}}_l\left(k=0\right)},$E36

$P^{{\rm B}}\left(B_l\right)=\frac{B_l\left(k=0\right)}{{{\rm B}}_l\left(k=1\right)+{{\rm B}}_l\left(k=0\right)}.$E37

Next, substitute Eqs. (37) and (38) into Eqs. (19) and (20) to obtain ${\kappa }_l^B$ and ${\nu }_l^B$. The right‐hand jitter $k_l^{B\text{R}}$ and left‐hand jitter $k_l^{B\text{L}}$ can now be specified by, respectively,

$k_l^{B\text{R}}=\frac{{\nu }_l^B}{{\kappa }_l^B}\ln \frac{1}{\xi B_l(k=0)}$E38

$k_l^{B\text{L}}={\nu }_l^B{\kappa }_l^B\ln \frac{1}{\xi B_l(k=0)}$E39

Finally, define the jitter left boundary $J_l^{B\text{L}}$ and right boundary $J_l^{B\text{R}}$ as, respectively,

$J_l^{B\text{L}}\cong {\widehat{n}}_l-k_l^{B\text{R}}$E40

$J_l^{B\text{R}}\cong {\widehat{n}}_l-k_l^{B\text{L}}$E41

and use the algorithm previously designed in the study of Munoz‐Minjares and Shmaliy [22] for the confidence masks based on the Laplace distribution.

3. Results

Figure 4.

and masks for the seventh chromosome taken from “159A–vs–159D–cut” of ROMA: (a) genomic location from 130 to 146 Mb and (b) genomic location from 146 to 156 Mb. Breakpoints î₁, î₆, î₇, î₉, î₁₀, î₁₂, and î₁₃ are well detectable because jitter is moderate. Owing to large jitter the breakpoints î₂, î₃, î₄, î₅, î₈, î₉, and î₁₁ cannot be estimated correctly. There is a probability that the breakpoints î₂, î₃, î₄, î₅, and î₁₁ do not exist. There is a high probability that breakpoint î₅ does not exist.

In this section, we test some CNA measurements and estimates by the algorithm developed in [22] based on the Laplace and Bessel approximations. In order to demonstrate the efficiency of the probabilistic masks and getting practically useful results, we exploit probes obtained by different technologies. First, we employ the results obtained with the HR‐CGH profile and test them by the probabilistic masks using Laplace distribution. We next demonstrate the efficiency of the Bessel‐based probabilistic masks versus the Laplace‐based masks for the probes obtained with the SNP profile.

3.1. HR‐CGH‐based probing

The first test is conducted in the three‐sigma sense suggesting that the CNAs exist between the UB and LB masks with high probability of $P=99.73\%.$ The tested HR‐CGH array data are available from the representational oligonucleotide microarray analysis (ROMA) [23]. The breakpoint locations are also given in [23]. Voluntarily, we select data associated with potentially large jitter and large segmental errors. For clarity, we first compute some characteristics of the detected CNAs and notice that the segmental estimates found by averaging [18] are in a good correspondence with [23]. The database processed is a part of the seventh chromosome in archive “159A–vs–159D–cut” of ROMA a sample of B‐cell chronic lymphocytic leukemia (B‐CLL). It is shown to have 14 segments and 13 breakpoints (Figure 4a and b). Below, we shall show that, owing to large detection noise, there is a high probability that some breakpoints do not exist.

It follows from Figure 4a that the only breakpoint which location can be estimated with high accuracy is $i_1$. Jitter in ${\widehat{i}}_6$ and ${\widehat{i}}_7$ is moderate. All other breakpoints have large jitter. It is seen that the UB mask covering second to sixth segments is almost uniform. Thus, there is a probability that the second to fifth breakpoints do not exist. If to follow the LB mask, the locations of the second to fourth breakpoints can be predicted even with large errors. At least they can be supposed to exist. However, nothing definitive can be said about the fifth breakpoint and one may suppose that it does not exist. It is also hard to distinguish a true location of the eighth breakpoint. In Figure 4b, $i_{10}$, $i_{12}$, and $i_{13}$ are well detectable owing to large segmental SNRs. The breakpoint $i_9$ has a moderate jitter. In turn, the location of $i_{11}$ is unclear. Moreover, there is a probability that $i_{11}$ does not exist.

3.2. SNP‐based probing

Our purpose now is to apply the probabilistic mask with SNP profile that represents the CNA with low levels of SNR. Specifically, we employ the probes of the first chromosome available from “BLC_B1_T45.txt” a sample of primary breast carcinoma.

Inherently, the more accurate Bessel‐based approximation extends the jitter probabilistic boundaries with respect to the Laplace‐based ones, especially for low SNRs. We illustrate it in Figure 5, where the estimates of the first chromosome were tested by

,

,

, and

for $\vartheta =3$ (confidence probability $P=99.73\%).$

In Figure 6, the masks

and

are placed in the vicinity of segment ${\widehat{a}}_{18}$ for several confidence probabilities: $\vartheta =0.6745(P=50\%)$, $\vartheta =1(P=68.27\%)$, $\vartheta =2(P=95.45\%)$, and $\vartheta =3(P=99.73\%)$. What the masks suggest here is that the CNA evidently exists with high probability, but the segmental levels and the breakpoint locations cannot be estimated with high accuracy, owing to low SNRs.

Figure 5.

Jitter left boundaries ${ℬ}_l^{ℬ\text{L}}$, $J_l^{\text{L}}$ and right boundaries $J_l^{ℬ\text{R}}$, $J_l^{\text{R}}$ for the breakpoint ${\widehat{i}}_2$ of first chromosome from sample BLC_B1_T45.txt (primary breast carcinoma). The probabilistic masks detect a breakpoint with a confidence probability $\vartheta =3(P=99.73\%)$.

Figure 6.

The and masks placed around the segmental level $a_{18}$ for several confidence probabilities [20]. Here, the CNA exists with high probability, but the segmental levels and the breakpoint locations cannot estimate with high accuracy.

Figure 7.

The confidence masks placed around $a_{10}$ for $\vartheta =0.6745(P=50\%)$ and $\vartheta =3(P=99.73\%)$. Masks and do not confirm an existence of segmental changes while and indicate a small change.

Figure 8.

The confidence masks , ${ℒ}_l^{\text{LB}},{ℬ}_l^{\text{UB}}$ and . placed around the breakpoint ${\widehat{i}}_{20}$ for confidence probabilities $\vartheta =0.6745$ and $\vartheta =3$ of first chromosome from sample BLC_B1_T45.txt. The confidence masks based on Laplace distribution cannot detect the breakpoint ${\widehat{i}}_{20}$.

A special case can also be noticed when the masks

and

are not able to confirm or deny an existence of segmental changes with high probability, owing to the inability of computing the Laplace‐based masks for extremely low SNRs. Figures 7 and 8 illustrate such situations. Just on the contrary, masks

and

can be computed for any reasonable SNR.

A conclusion that can be made based on the results illustrated in Figures 5–8 is that the Bessel‐based probabilistic masks can be used to improve estimates of the chromosomal changes for the required probability.

We finally notice that the computation time required by the masks to process the first chromosome from sample “BLC B1 T45.txt” with a length of n = 905215 was 2.634599 s using MATLAB software on a personal computer with a processor Intel Core i5, 2.5 GHz.

4. Discussion

We evaluate the breakpoints obtained by the projects representational oligonucleotide microarray analysis [23] and GAP [14] with the confidence masks. As has been shown before, not all of the detected chromosomal changes have the same confidence to mean that there is a probability that some breakpoints do not exist. In order to improve the CNA estimates for the required confidence, the following process can be used:

1. Obtain estimates of the CNA using the standard CBS algorithm [24, 25] or any other algorithm.

2. Compute masks

   
   and
   
   for the given confidence probability $P,\%$ and bound the estimates.

3. If the masks reveal double uniformities, in $\text{UB}$ and $\text{LB}$, in a GAP of any three neighboring breakpoints, then remove the intermediate breakpoint and estimate the segmental level between the survived breakpoints by simple averaging. The CNAs estimated in such a way will be valid for the given confidence $P,\%.$

Application of this methodology to the CNA structure detected in frames of the Project GAP is shown in Figure 9. Its special feature is a number of hardly recognized small chromosomal changes (Figure 9a). We test them by the proposed masks ${ℬ}_l^{\text{UB}}$ and ${ℬ}_l^{\mathrm{LB}}$. To this end, we first start with equal confidence probabilities of $P=50\%$ for each estimate to exist or not exist and find out that three breakpoints demonstrate no detectability. We remove these breakpoints and depict their locations with “×”. Reasoning similarly, we remove four breakpoints to retain only probable changes, by $P=75\%,$ nine breakpoints to show a picture combined with almost certain changes, by $P=93\%$, and 10 breakpoints in the three‐sigma sense, $P=99.73\%$. Observing the results, we infer that the masks are able to correct only the estimates obtained under the low SNRs. The relevant chromosomal sections S1–S7 are circled in Figure 9. It is not surprising because changes existing with high SNRs are seen visually. An estimator thus can easily detect them with high confidence.

Figure 9.

Improving estimates of the CNAs obtained in Project GAP [25] by removing some unlikely existing breakpoints: (a) original estimates, (b) even changes, $P=50\%$, (c) probable changes, $P=75\%,$ (d) almost certain changes, $P=93\%$, and (e) three‐sigma sense, $P=99.73\%.$

5. Conclusions

Modern technologies developed to produce the CNA profiles with high resolution still admit intensive white Gaussian noise. Accordingly, not one estimator even ideal is able to provide jitter‐free estimation of segmental changes. Thus, in order to avoid wrong decisions, the estimates must be bounded for the confidence probability. Jitter exists in the CNA's breakpoints fundamentally. When SNR is >1, it can statistically be described using the discrete skew Laplace distribution. Otherwise, if SNR is <1, the Bessel‐based approximation produces more accuracy. By the jitter distribution, it is easy to find a region within which the breakpoint exists for the required probability. Of practical importance are the confidence UB and LB masks, which can be created based on the segmental and jitter distributions for the given confidence probability. The masks can serve as an auxiliary tool for medical experts to make decisions about the CNA structures. Applications to probes obtained using the HR‐CGH and SNP technologies confirm efficiency of the confidence masks.