SOSCNV: a stochastic outlier selection algorithm-based method for copy number variation detection using hybrid strategy

CNV detection method based on Stochastic Outlier Selection (SOS)

The SOS method only requires an unlabeled data set as input. It quantifies the relationship between data points through correlation, and then calculates the outlier probability of each data point. The core idea is that the correlation is proportional to the similarity between two data points. That is, the smaller the similarity between two data points, the smaller the correlation. Once a data point has a relatively small correlation with other data points, this data point will be judged as an outlier. In the CNV detection progress, the RD value and Guanine and Cytosine (GC) content of the CNV region often have significant differences compared with the normal region, and the correlation between the copy number at a specific gene locus and that of the entire genomic data is also appropriate. Therefore, according to the magnitude of the correlation between the collection of genomic data, the CNV regions in the measured genome can be effectively identified. The variation detection process of SOSCNV is shown in Fig. 1.

In the workflow diagram of the SOSCNV method, the alignment file in Binary Alignment Map (BAM) format and the reference genome file in Fasta format are regarded as input files (Li & Durbin, 2009). Among them, the BAM file can be generated by aligning sequencing samples (Bonfield & Mahoney, 2013) and the reference genome through the BWA-MEM algorithm, and then sorted by the SAMTools software (Li et al., 2009). In the preprocessing stage, the RD value is extracted as a feature signal, and according to the alignment information in SR strategy, the breakpoints of CNV regions in the genome sequence are determined. In addition to the RD signal, the mapping quality (MQ) of the reads in the window is also counted. MQ is employed to measure the quality of read alignments and to exclude low-quality bases, thereby reducing noise in the sequencing data. In the segmentation stage, a two-stage segmentation strategy is given. First, the circular binary segmentation (CBS) algorithm is used for global segmentation (Venkatraman & Olshen, 2007). Then, the large fragments obtained by the CBS algorithm segmentation are re-segmented into small sub-fragments, Then, the SOS algorithm models the feature signal and calculate the anomaly score. Finally, the OTSU algorithm determines the threshold and identifies the final CNV (Goh et al., 2018).

Preprocessing

In the preprocessing stage, it covers extracting alignment information, dividing the bins, correcting GC deviations and searching for breakpoints (Du et al., 2025a). The RD signal is regarded as the input feature of the model, which is the most important feature information in the RD-based detection strategy. After calculating the RD value at each position on the genome sequence, further processing is needed to meet subsequent detection requirements.

(1) Dividing windows and calculating RD and MQ signals

The reference genome is divided into consecutive, equal-sized, non-overlapping bins with a width of 1,000 bp. For each bin, the sequencing read depth information from the BAM file is calculated and marked as the initial RD. If a bin in the reference genome contains “N” bases, the initial RD of this bin is set to a negative value and the bin is marked as an abnormal bin; the RD value of such bins is excluded from further statistical analysis. For normal bins that do not contain “N” bases, the average read count across all bases within the bin is calculated as the RD value for that bin and is used in subsequent detection analysis, as shown in Formula (1). In this formula, $R{D}_i$ represents the average read depth of the i-th bin, $r{c}_j$ is the read count at the j-th position within the bin, and $L_i$ is the size of the i-th bin. This approach effectively reduces noise-induced fluctuations and lowers computational cost.

(1) $$R{D}_i={\displaystyle \frac{{\sum }_1^{L_i}r{c}_j}{L_i}.}$$

In addition to the RD signal, the MQ of the read segments within each bin is also calculated, as shown in Formula (2). In this formula, $M{Q}_i$ represents the average mapping quality of the i-th bin, $m{q}_j$ denotes the mapping quality score at the j-th position within the bin, and $L_i$ is the size of the i-th bin. The MQ signal reflects the overall alignment reliability of reads in the bin, where higher MQ values indicate more confident alignments. By filtering out low-quality mapping positions and replacing unreliable base-level signals with the average MQ value, this approach effectively reduces the impact of noise on CNV detection results.

(2) $$M{Q}_i={\displaystyle \frac{{\sum }_1^{L_i}m{q}_j}{L_i}.}$$

(2) GC correction

GC bias (Benjamini & Speed, 2012) is a key factor that causes the inconsistency between RD signal and sequence depth. In regions with low or high GC content, the RD value will deviate, which will affect the analysis of subsequent detection results. The GC content also has a certain correlation with the RD value. Specifically, it is manifested that the GC content in the CNV region may have a large difference compared with that in the normal region. Therefore, in the SOSCNV method, the GC content is also used as the input feature signal of the anomaly detection model. The following method can effectively correct the RD value of each window. as shown in Formula (3). Among Formula (3), $R{D}_i^{\mathrm{\prime }}$ and $R{D}_i$ respectively represent the corrected RD value and the original RD value of the i-th bin, $R{D}_{avg}$ represents the average RD value of all bins, and $R{D}_{gc}$ represents the average RD value of all bins with the same GC content as this bin.

(3) $$R{D}_i^{\mathrm{\prime }}={\displaystyle \frac{R{D}_{avg}}{R{D}_{gc}}R{D}_i.}$$

(3) Searching for possible breakpoints

SOSCNV employs a split-read strategy based on the CIGAR string to identify breakpoints. For alignment information in BAM files, it prioritizes read records that contain soft clipping or hard clipping operations. The presence of clipping indicates that part of the sequencing read could not be successfully aligned to the reference genome during the mapping process, and these unmapped segments may suggest potential CNV. Formula (4) summarizes the possible combinations of matching and clipping operations in the Compact Idiosyncratic Gapped Alignment Report (CIGAR) string. “M” denotes the continuous segment of the sequencing read that aligns to the reference genome, while “S” and “H” represent soft clipping and hard clipping in the CIGAR string, respectively. $xM$ indicates that x bases are aligned to the reference genome; $yS$ means that y bases remain unaligned but are kept in the read; and $yH$ means that y bases are unaligned and removed from the read. When the clipped segment occurs before the aligned region ( $y\mathrm{S}x\mathrm{M}$ and $yHxM$), the breakpoint is assigned to the start position of the aligned block, i.e., $a=R.pos$. Conversely, when the clipped segment occurs after the aligned region ( $x\mathrm{M}y\mathrm{S}$ and $xMyH$), the breakpoint is assigned to the end position of the aligned block, i.e., $b=R.pos+x-1$.

(4) $$CIGAR=\{\begin{array}{c}xMyS\hfill \\ xMyH\hfill \\ ySxM\hfill \\ yHxM\hfill \end{array}.$$

Global segmentation and local segmentation

After the preprocessing is completed, the segmentation step will be executed, which includes two stages: global segmentation and local segmentation, with the aim of identifying consecutive regions with the same RD signal. Global segmentation is performed using the circular binary segmentation (CBS) algorithm. The CBS algorithm first searches for a continuous sequence of genomic windows. Based on the maximum statistic, the method compares the average RD value within a continuous window region to the average RD value of the remaining windows. If the difference between the two exceeds a predefined threshold, the position is considered a potential breakpoint. The region is then partitioned into new segments using this breakpoint as the boundary. This process is recursively repeated across the entire genome until no further significant change points can be detected, ultimately achieving multiple segmentations of the genome. Based on the divided segments, the RD value and GC content are recalculated, as shown in Formulas (5) and (6). Among formulas, $RD{s}_i$ and $GC{s}_i$ respectively represent the RD value and GC content in the i-th fragment, $size\mathrm{\_}{s}_i$ represents the size of the i-th fragment, and $R{D}_i$ and $G{C}_i$ respectively represent the RD value and GC content of each window contained in the fragment. After segmentation is completed, the total number of samples changes from the number of windows to the number of segments divided.

(5) $$RD{s}_i={\displaystyle \frac{{\sum }_{j=1}^{size\mathrm{\_}{s}_i}R{D}_j}{size\mathrm{\_}{s}_i}.}$$

(6) $$GC{s}_i={\displaystyle \frac{{\sum }_{j=1}^{size\mathrm{\_}{s}_i}G{C}_j}{size\mathrm{\_}{s}_i}.}$$

After the global segmentation is completed, proceed with local segmentation. Due to the substantial variation in the lengths of candidate segments produced by CBS segmentation, SOSCNV employs a local re-segmentation strategy to reduce the influence of length heterogeneity on subsequent analyses and to improve the local resolution of CNV detection. In this process, the candidate segments are further partitioned into uniformly sized sub-segments with fixed lengths. First, the length of the sub-segments needs to be defined. The size of sub-segment is related to the resolution of CNV detection. A smaller sub-segment can provide a higher breakpoint resolution, but it will lead to an increase in false positive events; while a larger sub-segment can improve the recognition accuracy, but it will cause additional false negative events. The length threshold of sub-segments is defined as $L_{bin}$. The parameter $L_{bin}$ is flexible and can be selected based on the characteristics of the sequencing data. In practice, $L_{bin}$ may be configured to 500, 1,000, 3,000 bp, or larger values to balance breakpoint resolution and noise robustness across different datasets. Subsequently, the segments with a length greater than $L_{bin}$ are re-segmented into multiple consecutive sub-segments with a length of $L_{bin}$. If the length of the sub-segment does not reach $L_{bin}$, this fragment is merged into the previous sub-segment, and the segmentation continues on the concatenated segment. For sub-segments that do not require merging, their sequencing signals (RD and GC) are computed directly from the bins they cover. For sub-segments generated by merging a short residual fragment with the previous segment, the sequencing signals are obtained by averaging the RD and GC values of all bins contained in the merged interval.

Abnormal score calculation and confirmation threshold

After local segmentation is completed, the next step is to calculate the abnormal score. The purpose of this step is to quantify the degree of abnormality of each sub-segment, which helps in identifying potential CNVs more accurately. In anomaly detection algorithms, anomaly scores can be used to represent the degree of anomaly of each data sample. The anomaly score in the SOS algorithm is the probability value of anomalies between 0 and 1, so there is no need for additional normalization processing, avoiding the problem of large differences in anomaly score ranges on genomic data with different coverage levels. The outlier score represents the probability of a sample point being an outlier, and the higher the score, the greater the likelihood of it being an outlier. The SOSCNV algorithm primarily uses the preprocessed read depth signal and GC-content signal as the initial inputs. The abnormal score calculation process of the SOSCNV model can be seen as a transformation process of four matrices, and the specific steps are as follows:

• (1)

  Let $M=[M_{ij}]\in R^{N\times 2}$ denote the input segmentation matrix, where n is the number of segments, and the two columns represent the preprocessed RD and GC signals;

• (2)

  Calculate the dissimilarity matrix D based on the distance between each segmentation point;

• (3)

  Use the matrix D to calculate the correlation between segmentations and represent it with the matrix P;

• (4)

  Calculate the correlation probability matrix Q based on matrix P;

• (5)

  Calculate the anomaly score for each segmentation from the correlation probability matrix Q.

In Step 1, the input data sample set includes two features, namely the RD signal and GC content of each sub-segment, which are represented by a matrix M of size $n\times 2$, where n denotes the number of sub-segments, as shown in Formula (7). Each row of the matrix M can be regarded as a sample point and is represented by a two-dimensional feature vector x_i = [ $r{d}_i$, $g{c}_i$].

(7) $$M=\left(\begin{array}{c}r{d}_{1}g{c}_1\\ r{d}_{2}g{c}_2\\ r{d}_{3}g{c}_3\\ \dots \\ r{d}_{n}g{c}_n\end{array}\right).$$

In Step 2, the two features of the data points are used for the calculation of the dissimilarity between the data points. The dissimilarity is obtained through the Euclidean distance between two data points in the two-dimensional Euclidean space and is a non-negative scalar. The calculation process of the dissimilarity $d_{ij}$ between the data points $m_i$ and $m_j$ is given by Formula (8).

(8) $$d_{ij}=\sqrt{{\left(r{d}_j-r{d}_i\right)}^2+{\left(g{c}_j-g{c}_i\right)}^2}.$$

It can be learned from the above formula that the dissimilarity matrix D is a symmetric matrix, that is, $d_{ij}$ = $d_{ji}$. At the same time, the distance from a data point to itself is 0. The dissimilarity matrix D ( $n\times n$) can be expressed in the form of Formula (9).

(9) $$D=(\begin{array}{cccc}0& d_{12}& \dots & d_{1n}\\ d_{21}& 0& \dots & d_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ d_{n1}& d_{n2}& \cdots & 0\end{array}).$$

In Step 3, the similarity relationship between two data points is quantified through the degree of association. The higher the similarity between the two data points, the greater the degree of association. Conversely, the greater the dissimilarity, the smaller the degree of association. The calculation formula for the degree of association $p_{ij}$ between the two data points is shown in Formula (10).

(10) $$p_{ij}=\{\begin{array}{c}\exp \left(-{d_{ij}}^2/2{{\alpha }_i}^2\right)\\ 0\end{array}\begin{array}{c}\\ \end{array}\begin{array}{c}\\ \end{array}\begin{array}{c}if\\ if\end{array}\begin{array}{c}\\ \end{array}\begin{array}{c}i\ne j\\ i=j\end{array}.$$

Here, ${{\alpha }_i}^2$ represents the variance of the data point m_i. In the SOS algorithm, each data point corresponds to a variance. The variance value is related to the density of the data point. Points with high density correspond to a smaller variance. The smaller the variance, the faster the attenuation of the degree of association; conversely, the larger the variance, the smaller the influence of dissimilarity on the degree of association. The SOS algorithm uses an adaptive method to determine the variance of each data point so that each data point has the same number of effective neighbors. The number of neighbors of a data point is called the complexity and is also the only parameter of the SOS algorithm. Each data point can have 1 to $n-1$ effective neighbors. By calculating the degree of association between data points, the association matrix P ( $n\times n$) is obtained, as shown in Formula (11). The diagonal elements in matrix P are 0 because there is no association between a data point and itself. Different from the dissimilarity matrix D, matrix P is asymmetric. This is because the variance of each data point is generated through an adaptive method, and the probability that two data points have the same variance is extremely low, thereby causing the asymmetry of the association degree ( $p_{ij}\ne p_{ij}$).

(11) $$P=(\begin{array}{cccc}0& p_{12}& \dots & p_{1n}\\ p_{21}& 0& \dots & p_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ p_{n1}& p_{n2}& \cdots & 0\end{array}).$$

In Step 4, the association probability matrix Q is obtained by normalizing the row vectors in the association matrix P. The calculation of the association probability value $q_{ij}$ between two data points is shown in Formula (12).

(12) $$q_{ij}={\displaystyle \frac{p_{ij}}{{\sum }_{k=1}^n{p}_{ik}}.}$$

Similar to matrix P, the correlation probability value $q_{ij}$ between the data points themselves is 0. Based on the correlation probabilities between the data points, the anomaly score of each sample is calculated, as shown in Formula (13). Among Formula (13), $C_o$ represents the anomaly class, and the anomaly $score\left(m_i\right)$ is a probability value ranging from 0 to 1. When the correlation probability value between a sample point and all other sample points is low, it is highly likely to be an abnormal point. Generally speaking, the segment represented by sample points with larger anomaly scores are potential CNV regions.

(13) $$score\left(m_i\right)=p\left(m_i\in C_o=\prod _{j\ne i}1-q_{ji}\right).$$

For different genomic data, the abnormal scores calculated by the SOSCNV model are also different. If a threshold of fixed size is used to classify the abnormal scores, it will produce certain deviations and cause inaccurate subsequent detection results. So, this article uses the OTSU to automatically obtain the thresholds applicable to various different genomic samples.

The OTSU algorithm is a global binary segmentation algorithm that dynamically obtains the threshold by traversing all the scores in the interval to maximize the variance between the normal class and the abnormal class. The abnormal score obtained by the SOS model is a probability value between 0 and 1. The abnormal score values of large abnormal samples are usually larger, and their number is much smaller than that of normal samples. Therefore, when using the OTSU algorithm to select the threshold, a range interval can be specified instead of starting from the minimum abnormal score to traverse until the threshold T that maximizes the variance between the normal class and the abnormal class is found. To improve computational efficiency, the OTSU algorithm can search for the threshold T that maximizes the between-class variance within a specified score interval (e.g., [0.5, 0.7]) where abnormal scores are more likely to fall. This avoids unnecessary traversal of very low or high scores while still adapting dynamically to the sample distribution. After T is determined, the abnormal score value of each fragment is compared with T, and the fragment with a score higher than T is determined as the CNV region. Then, the RD value corresponding to the CNV region is compared with the average RD value of the normal region to distinguish whether each variation region is copy number amplification or deletion.

Simulation experiment analysis

The results are shown in Figs. 2 and 3. The indicator values in the figure are the average values of the evaluation quantities in 50 samples for each configuration. The gray curve represents the F1-score. To ensure the fairness of the experiment, each method uses its default parameters.

Regardless of whether the simulated data are from Chromosome 13 or Chromosome 21, the SOSCNV method consistently achieves the highest F1-scores across all combinations of sequencing depths and tumor purities, demonstrating superior overall detection performance. Overall, SOSCNV, Manta, and TARDIS exhibit strong robustness, as their detection performance is relatively unaffected by changes in sequencing depths, tumor purities, and CNV lengths. Specifically, on the simulated data from Chromosomes 13 and 21, Manta’s F1-scores consistently range between 0.6 and 0.7. For TARDIS, the F1-scores on Chromosome 13 also range from 0.6 to 0.7, whereas on Chromosome 21, they slightly increase to a range of 0.6–0.8 with greater fluctuation, suggesting that TARDIS performs more stably when detecting longer CNVs. In contrast, SOSCNV achieves F1-scores between 0.7 and 0.8 on Chromosome 13, and an even higher range of 0.75–0.85 on Chromosome 21, indicating better detection performance. Across both chromosomes, SOSCNV outperforms Manta and TARDIS. Meanwhile, RSICNV and Matchclips2 show opposite performance trends: RSICNV achieves higher F1-scores under low sequencing coverage than at higher depths, whereas Matchclips2 performs better at higher coverage levels, with noticeably lower F1-scores at 10 $\times$.

In the simulated data, Chromosome 13 is used to construct longer CNVs, whereas Chromosome 21 is used to construct shorter CNVs. SOSCNV consistently maintains high and stable F1-scores in both scenarios, achieving scores of 0.7–0.8 on Chromosome 13 and even higher scores of 0.75–0.85 on Chromosome 21. This indicates that SOSCNV is not only effective for long CNVs but also demonstrates enhanced sensitivity when detecting shorter CNVs, a scenario where many traditional methods typically struggle.

Moreover, the F1-scores of SOSCNV across different sequencing depths and tumor purities exhibit a narrower variance compared to other tools, suggesting that its performance is less affected by these common experimental variables. The method achieves a balanced trade-off between precision and recall, as reflected by its favorable positioning on the F1-score contour lines, without sacrificing one for the other. In contrast, other tools such as Manta and TARDIS tend to exhibit slightly lower F1-scores and greater fluctuations depending on CNV length and coverage, while RSICNV and Matchclips2 show more pronounced performance changes when sequencing depth varies.

Collectively, these observations highlight that SOSCNV offers not only superior average detection accuracy but also stronger robustness across different CNV lengths, sequencing coverages, and tumor purity conditions, making it particularly suitable for real-world applications where CNV size and sequencing quality can vary widely.

Real experimental analysis

To compare the detection performance of SOSCNV with several other tools on real data, samples from the Yoruba family trio in the 1,000 Genomes Project were selected, including NA19238 (mother), NA19239 (father), and NA19240 (daughter). In addition, the HG002 sample was downloaded from NCBI as another real dataset for evaluation. SOSCNV and four other tools were evaluated based on these four real samples. For performance assessment, the Database of Genomic Variants (DGV, https://dgv.tcag.ca/dgv/app/home) and Genome in a Bottle (GIAB, https://www.nist.gov/programs-projects/genome-bottle) were used as ground truth references. In the evaluation, the ground truth was derived from the GRCh38 reference genome, using the version released on February 25, 2020 (https://dgv.tcag.ca/dgv/app/downloads?ref=GRCh37/hg19; accessed February 2026). For the ground truth of the real data HG002, we used the GIAB benchmark released in November 2024 (https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/mosaic_v1.00/GRCh38/SNV/).

Figure 4 illustrates the performance of five CNV detection tools across four real human samples. The three bar plots respectively show the average precision, recall, and F1-score over the 22 chromosomes. As shown in the figure, SOSCNV achieves substantially higher precision in the NA19238 and NA19240 samples compared to other tools, indicating superior detection accuracy. Although Manta exhibits relatively high precision in the NA19240 sample, the number of detected CNVs is too small, resulting in a very low recall and thus lowering its overall F1-score. Due to the low sequencing depth of the three real samples, Mathclips2 is almost unable to effectively detect copy number variations. In the real data sample HG002, SOSCNV not only achieved the highest precision but also the highest F1-score, further demonstrating its robustness across different datasets. Overall, SOSCNV achieves the highest F1-scores across all four real samples, demonstrating that it outperforms the other four tools not only in simulated datasets with varying CNV lengths, sequencing depths, and tumor purities, but also in real-world data. Due to the large number of variants included in the DGV ground truth, although some detection tools achieve relatively high precision on certain variants, the number of true positives detected is still small compared with the total number of ground truth variants. Consequently, the recall of all five tools remains relatively low.

NA19238, NA19239, and NA19240 belong to the same family, making the Mendelian inconsistency rate a reliable metric for CNV evaluation. If a CNV detected in the offspring overlaps (50% reciprocal overlap) with a CNV of the same type in either parent, then the CNV is defined as the Mendelian consistent CNV, else the CNV is defined as the Mendelian inconsistent CNV. The rate of Mendelian inconsistent CNVs in the offspring is defined as Mendelian inconsistency rate (Liu, Wu & Wang, 2021). Table 1 shows the Mendelian inconsistency rates of the five CNV detection tools on the real data samples. RSICNV exhibited the lowest Mendelian inconsistency rate among all tools, which may be related to the larger number of CNVs detected by RSICNV: it identified the highest number of CNVs in NA19238, NA19239, and NA19240. In contrast, Manta, Matchclips2, and TARDIS showed relatively high Mendelian inconsistency rates, all exceeding 30%, while SOSCNV’s inconsistency rate was lower than that of these three tools.

Evaluation of detection time

Figure 5 compares the runtime of the five CNV detection tools across the three real data samples. All five detection tools were executed on the same computer, which was running Ubuntu 22.04.5 LTS, equipped with 16.0 GiB of RAM and a 13th Gen Intel® Core™ i7-13700F processor with 24 cores. As shown in the figure, Manta requires the longest runtime, whereas Matchclips2 is the fastest among the five tools. SOSCNV’s runtime is the second longest, just below Manta. The overall runtime of SOSCNV includes both the preprocessing and detection stages. While it is slightly longer than some other tools, this also represents an area for future improvement. Further optimization of runtime, particularly for large-scale datasets, would enhance its applicability in research and clinical settings.