
==== Front
Brief Bioinform
Brief Bioinform
bib
Briefings in Bioinformatics
1467-5463
1477-4054
Oxford University Press

10.1093/bib/bbae022
bbae022
Review
AcademicSubjects/SCI01060
On the core segmentation algorithms of copy number variation detection tools
https://orcid.org/0009-0004-7363-6509
Zhang Yibo Key Laboratory of Biomedical Information Engineering of Ministry of Education and Department of Biomedical Engineering, School of Life Science and Technology, Xi’an Jiaotong University, Xi’an, China

Liu Wenyu Key Laboratory of Biomedical Information Engineering of Ministry of Education and Department of Biomedical Engineering, School of Life Science and Technology, Xi’an Jiaotong University, Xi’an, China

Duan Junbo Key Laboratory of Biomedical Information Engineering of Ministry of Education and Department of Biomedical Engineering, School of Life Science and Technology, Xi’an Jiaotong University, Xi’an, China

Corresponding author: Junbo Duan, School of Life Science and Technology, Xi’an Jiaotong University, 28 Xianning West Road, Xi’an 710049, China. Tel.: 029-82668668; E-mail: junbo.duan@mail.xjtu.edu.cn
Yibo Zhang and Wenyu Liu contributed equally to this work.

3 2024
10 2 2024
10 2 2024
25 2 bbae02217 7 2023
26 10 2023
11 1 2024
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Shotgun sequencing is a high-throughput method used to detect copy number variants (CNVs). Although there are numerous CNV detection tools based on shotgun sequencing, their quality varies significantly, leading to performance discrepancies. Therefore, we conducted a comprehensive analysis of next-generation sequencing-based CNV detection tools over the past decade. Our findings revealed that the majority of mainstream tools employ similar detection rationale: calculates the so-called read depth signal from aligned sequencing reads and then segments the signal by utilizing either circular binary segmentation (CBS) or hidden Markov model (HMM). Hence, we compared the performance of those two core segmentation algorithms in CNV detection, considering varying sequencing depths, segment lengths and complex types of CNVs. To ensure a fair comparison, we designed a parametrical model using mainstream statistical distributions, which allows for pre-excluding bias correction such as guanine-cytosine (GC) content during the preprocessing step. The results indicate the following key points: (1) Under ideal conditions, CBS demonstrates high precision, while HMM exhibits a high recall rate. (2) For practical conditions, HMM is advantageous at lower sequencing depths, while CBS is more competitive in detecting small variant segments compared to HMM. (3) In case involving complex CNVs resembling real sequencing, HMM demonstrates more robustness compared with CBS. (4) When facing large-scale sequencing data, HMM costs less time compared with the CBS, while their memory usage is approximately equal. This can provide an important guidance and reference for researchers to develop new tools for CNV detection.

copy number variants
next generation sequencing
circular binary segmentation
hidden Markov model
National Natural Science Foundation of China 10.13039/501100001809 61771381
==== Body
pmcINTRODUCTION

Copy number variants (CNVs) [1] refer to alterations in copy number (CN) (usually >50 base pairs (bp) [2]) of a segment of deoxyribonucleic acid (DNA) sequence in a genome, including amplification or reduction in CN. They are often associated with biological phenomena such as genetic disease [3], adaptive evolution [4] and gene expression [2]. The development of high-throughput sequencing technologies [or next-generation sequencing (NGS)] has provided an important approach for the rapid detection of individual genomic CNVs [5–8].

In NGS-based CNV detection, the core task is to segment genomic regions and CN typing of sequences based on the observed changes in the data [9]. For paired-end whole-genome sequencing data, there are mainly five methods to detect CNVs [2], namely (1) paired-end mapping [8], (2) split read [1], (3) reads depth (RD) [10], (4) assembly [11] and (5) integration-based detection methods [12].

Among the five detection methods mentioned above, the RD-based detection method is currently the most widely used in the field of CNV detection [13]. The RD-based method uses the read counts information obtained from NGS data after alignment to detect CNV. The underlying assumption is that the RD values in genomic regions are CN-related, with CN duplication regions having higher RD values than normal regions and CN deletion regions having lower RD values than normal regions [14]. The advantage is that larger segments of complex CNV regions can be detected [12].

In the RD methods, CN detection usually starts with segmentation, where the genome is partitioned into non-overlapping regions, with the assumption that each region shares the same CN. This process aims to identify the boundaries of segments where CNV events occur, often referred to as breakpoints. As a result, the efficacy of CNV detection depends heavily on the segmentation algorithm employed [12].

Based on this, we have conducted a comprehensive study on the RD-based CNV detection tools over the past decades (Table 1), and we found that the most widely used segmentation algorithms are circular binary segmentation (CBS) and hidden Markov models (HMMs) [15–50]. To fairly compare the algorithm performance, we introduced a parametric modeling approach to simulate experimental data [51]. Finally, we systematically studied the advantages, disadvantages and applicability of core segmentation algorithms, such as CBS and HMM, in the hope to provide important guidance and reference for the development of new tools.

Table 1 List of CNV detection tools

Method	Segmentation algorithm	Statistical distribution	Language	Publishing year	
SegSeq [15]	CBS	Poisson	Matlab	2009	
ReadDepth [16]	NB	R	2011	
CNAnorm [17]	Normal	R	2012	
cn.MOPS [18]	Poisson	R	2012	
CNVrd2 [19]	Poisson	R	2014	
CNVkit [20]	NB	R	2016	
Canvas [21]	Normal	C#	2016	
iCopyDAV [22]	NB	R	2018	
CNV-LOF [23]	\	Python/R	2019	
ClinCNV [24]	Normal	R	2019	
CamCNV [25]	Poisson	R	2020	
IhybCNV [26]	\	Python/R	2021	
CNAseg [27]	HMM	NB	R	2010	
JointSLM [28]	Normal	R	2011	
GENSENG [29]	Mixture	C++	2012	
ERDs [30]	Poisson	Perl/C	2012	
PennCNV-Seq [31]	Mixture	Perl/C	2017	
CLImAT-HET [32]	NB	C/Matlab	2017	
CNV-HMM [33]	NB	R	2019	
JAX-CNV [34]	Poisson	C++	2021	
Coverage-	Normal	Python	2022	
Master [35]	
SCONCE [36]	NB	C++	2022	
CNV-seq [37]	Robust statistical model	Poisson	R	2009	
RDXplorer [38]	EWT	Normal	Python	2009	
mrCaNaVar [39]	mrFAST	\	C	2009	
CNVnator [40]	Mean-shift	Normal	C++	2011	
BIC-Seq [41]	MLE	\	Perl/R	2011	
Control-FREEC [42]	LASSO	Normal	C++	2012	
RSICNV [43]	LRS	NB	C++	2014	
GROM-RD [44]	Sliding	Beta	C	2015	
window	
CNV-RF [45]	Random	\	Perl/R	2016	
forest	
XCAVATOR [46]	SLM	Poisson	R	2017	
DCC [47]	DCC	Normal	C++	2018	
CNV_IFTV [48]	Isolation	\	Python	2019	
forest	
DpGMM [49]	cghFLasso	Normal	Python	2020	
ifCNV [50]	Isolation	\	Python/R	2022	
forest	

MATERIALS AND METHODS

The core step in CN variation detection is to detect the location in the chromosome where the CN changes. Various segmentation methods for CN variation have been developed, mainly based on breakpoint detection and/or smoothing techniques. For RD-based whole-genome CNV detection tools, there are two main segmentation algorithms: (1) CBS [52] and (2) HMMs [53]. As shown in Table 1, CBS and HMM are dominant segmentation algorithms in CNV detection and are used by many detection tools. Other minor segmentation algorithms used by other CNV detection tools differ from one another, making it difficult to make a fair comparison under a unified framework. Therefore, in this study, only CBS and HMM are considered for comparison. Table 1 shows the RD-based whole-genome CNV detection tools over the last decades, and it should be noted that some newly developed tools are not listed in Table 1; e.g. GATK-gCNV [54] was designed for whole-exome sequencing data analysis in which the window size varies according to exons; mCNA [55] employs feature other than RD signals to detect CNV, etc.

Circular binary segmentation

CBS is a further refinement of binary segmentation by Olshen et al. [52] in 2004. The algorithm focuses on detecting breakpoints, which are the genomic locus where the CN changes. The algorithm starts with the whole chromosome and recursively segments it by detecting breakpoints and stops when none are found in any segment. The goal of the algorithm is to identify the breakpoints that segment a chromosome into segments of constant CN. Once the chromosomes have been segmented, additional information such as chromosome ploidy can be used to estimate the CN of each segment. The process is shown below.

(i) Assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_1,\dots, {X}_n$\end{document} represents the data for n markers in the given genomic segment. The test statistic is defined as the maximum T-statistic, given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $T=\underset{i<j\le n}{\max}\kern0.1em \left|{T}_{ij}\right|$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${T}_{ij}$\end{document} represents the two-sample T-statistic that compares the mean of the observations with indices from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i+1$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document} with the mean of the remaining observations.

(ii) If the P-value is less than the threshold \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\alpha$\end{document} (usually 0.01), we declare the change to be statistically significant and estimate the location of the breakpoint as maximizing the test statistic for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document} and\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document} as well as applying the algorithm recursively to the three resulting segments: from 0 to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document}; from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document} and from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $n$\end{document} until no further break points are found.

(iii) After finding partition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${I}_1,\dots, {I}_m$\end{document} of the data \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_1,\dots, {X}_n$\end{document} with CBS, the original single partition test is applied to each triple \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${I}_r$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${I}_{r+1}$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${I}_{r+2}$\end{document} to confirm if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${I}_{r+1}$\end{document} is a breakpoint. If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${I}_{r+1}$\end{document} is not a breakpoint, it is dropped.

Considering the algorithm performance as well as optimization issues, in the subsequent comparison experiments, we employed the migration of the relevant part of the segmentation algorithm from the existing software package to implement the target algorithm. Based on existing research, we found that the majority of CN variation detection tools based on CBS rely on the R package DNAcopy. This package was developed by Olshen [56], the author of the CBS algorithm, and is regularly maintained on the software website, providing easy access to the original code. This allows researchers to use the package within a programming environment, facilitating customization of their analysis and generation of publication-quality figures. Therefore, in the subsequent experiments in this paper, we utilized the software package DNAcopy for invoking the CBS algorithm.

The DNAcopy package comprises three primary functions, namely ‘CAN’, ‘segment’ and ‘plot.DNAcopy’. The CAN function generates a ‘CN array’ data object for analyzing DNA CNs in programs such as CBS. The segment function segments DNA CN data into regions of estimated equal CN using CBS. The plot.DNAcopy function plots the data and results from segmenting a CNA object. Regarding setting parameters, we focus on adjusting the parameters of the segment function. The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\alpha$\end{document} value was adjusted to 0.01 to control the significance threshold of segment breakpoints and improve the detection rate of CNV, while undo.splits was adjusted to true, which means that the segmentation was post-processed to remove small segments potentially stemming from noise.

Hidden Markov model

The HMM proposed by Baum et al. [57] described a statistical model of Markov processes containing hidden states, which was originally applied to the field of speech recognition. The model assumes that there is an underlying sequence of hidden states that generates the observed sequence of events. Hidden states cannot be directly observed, but observable events are assumed to depend on them. HMM consists of three main components: an initial probability distribution of hidden states, a transition probability matrix describing the probability of transitioning from one hidden state to another and an emission probability distribution describing the probability of observing a particular event given a hidden state.

HMM has been applied to CNV detection studies [27–36]. In these studies, the genome is generally divided into fixed length, non-overlapping windows, and the CN status within each window is generally positively correlated with RD in that window. Therefore, when the HMM is applied to CNV detection, the CN state within each window is the hidden state and RD within each window is the observation. The transfer probability is the conditional probability that the CN state of one window corresponds to a certain state, while the CN state of another window corresponds to a different state. As the distant between the two windows increases, the likelihood that the CN states of the two windows do not match also increases. The emission probability distribution is the conditional probability distribution of sequenced segments based on the CN status. The process is shown below.

(i) Defining the hidden states: The state represents the underlying CN. We can define a set of hidden states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $q=\left({q}_1,..,{q}_K\right)$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $K$\end{document} is the total number of all possible CNs.

(ii) Defining the initial state probabilities_\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\pi$\end{document}: Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\pi}_i$\end{document} be the initial state probability, the probability that the state of the first window is state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document}.

(iii) Defining the state transition probabilities_\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A$\end{document}: The underlying hidden Markov chain is defined by state transitions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $P\left({q}_t\mid{q}_{t-1}\right)$\end{document} and is represented by a stochastic transition matrix. The specific setting of the state transition matrix is described in detail later.

(iv) Defining the emission probabilities_\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $B$\end{document}: There are various ways to set the emission probability of RD signal, and the common ones are Normal distribution [28, 35], Poisson distribution [30, 34], negative binomial distribution [27, 32, 33, 36] or mixed distribution [29, 31]. In general, the negative binomial distribution is widely used for HMM emission probabilities in CN variation detection. The negative binomial distribution has the characteristics of high data volatility and low probability of variation in the process of deep variation [16], which will be described in detail in the Parametric modeling of RD data section.

(v) Training: The Baum–Welch algorithm [58] is the most commonly used algorithm in HMM parameter optimization, which implements the EM algorithm by executing two steps alternately. First, the expectation step calculates the expectations of the data hidden variables (hidden states), and then the M step (Maximization step) uses these expectations to calculate the maximum likelihood estimates of the model parameters.

(vi) Calling CN: The Viterbi algorithm [59] is usually the last step of the HMM algorithm, which calculates the most likely sequence of hidden states in the model given the HMM model \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\lambda =\left(A,B,\pi \right)$\end{document} and the sequence of observations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $O=\left({O}_1,{O}_2,\dots, {O}_n\right)$\end{document}. The main concept is to obtain the optimal sequence of states by recursively extrapolating the most probable sequence of states at each moment.

Due to the absence of a standardized HMM package for CNV detection in sequencing data, we extracted the HMM segmentation algorithm from an existing HMM package. HMMcopy is an R package developed by Lawrence [60] for unbiased CN detection based on high-throughput whole-genome sequencing data, which is also regularly maintained by the authors. It facilitates the partitioning of sequencing data into windows, guanine-cytosine (GC) content correction, CN segmentation and visualization of results. We have mainly migrated the HMMsegment function in the HMMcopy package to implement our own HMM segmentation algorithm and use the Baum–Welch algorithm for parameter optimization.

To improve the segmentation performance of the HMM, we made two changes to the algorithm. First, we made the number of hidden states adjustable. Second, we rewrote the transition matrix to make it more realistic, as CNs in whole-genome sequencing data do not typically jump directly from amplification to deletion or vice versa. If the number of hidden states is three, then the transition matrix is defined as follows:

(1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} A=\left[\begin{array}{ccc}1-d& d& 0\\{}d& 1-2d& d\\{}0& d& 1-d\end{array}\right], \end{equation*}\end{document}

where d is the probability that the CN state transits. Here, we set the transition probability of jumping from deletion to amplification and from amplification to deletion between two consecutive windows as 0, given that this scenario is almost impossible in genomes.

Parametric modeling of RD data

Traditional sequencing data simulation usually follows: downloading standard human reference genome (e.g. hg38) [61], simulating the occurrence of CNVs to generate test genome with CNVs, simulating reads to mimic sequencing, mapping data to reference genome with alignment tools such as MAQ [62] or Bowtie [63], bias correction such as GC content and mappability and generating RD simulation signals. The process of data alignment and bias correction inevitably introduces systematic errors into the RD signal, making a fair comparison of the performance of segmentation algorithms more challenging. Meanwhile, the alignment of short reads is often time-consuming and requires extensive computing resources. Therefore, we opted to use a parametric model to simulate the genome-wide RD signals directly. This approach avoids the challenges associated with data alignment and bias correction, making a fairer comparison of segmentation algorithm performance.

In previous literature, Jeng et al. [64] considered the CNV detection problem as a general problem of detecting short segments from a long string of Gaussian noise data. These short segments are characterized by their boundaries known as breakpoints, which are distributed sparsely along the whole genomic locations. Cai et al. [51] assumed additivity and symmetric noise with constant variance across the genomic data. However, this assumption may be compromised by the variation in GC content, DNA repeats and sequencing errors. Through further studies, when assuming a constant read sampling rate across the genome, Xie et al. [37] showed that RD data follows a Poisson distribution, and Miller et al. [16] further showed that the negative binomial distribution can be used to model the RD data in CNV detection. Therefore, Saran et al. [43] integrated a natural exponential family with a quadratic variance function to model the RD data, including the Poisson and negative binomial distribution as special cases. Table 1 also shows the statistical distributions used in literature, revealing that the normal, Poisson and negative binomial distributions are the most frequently used.

Based on these studies, we assume that the genome of a given individual is sequenced and short reads are mapped to a reference genome. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_1,\dots, {X}_n$\end{document} be the genomic serial RD data observed at window \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i=1,\dots, n$\end{document}. Suppose that there are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $q$\end{document} CNVs with q possibly increasing with n, and let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathbb{I}=\left\{{I}_1,\dots, {I}_q\right\}$\end{document} be the collection of the corresponding true CNV segments. We assume that the number, locations and mean values of the CNV segments are unknown. Finally, we modeled the RD data parametrically using Poisson distribution and negative binomial distribution as shown in the following text.

Poisson distribution

When the RD data \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_1,\dots, {X}_n$\end{document} are assumed to follow a Poisson distribution, they can be expressed as in Equation (2):

(2) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {X}_i\sim \mathrm{Poisson}\left({\lambda}_i\right),{\lambda}_i={\lambda}_0\frac{c_i}{2}, \end{equation*}\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\lambda}_0$\end{document} is the average RD of the normal genomic region (with CN \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $c=2$\end{document}), which is related to the sequencing depth; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${c}_i$\end{document} is a non-negative integer that represents the value of the CN state of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_i$\end{document} at the location \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document}.

Negative binomial distribution

When the RD data \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_1,\dots, {X}_n$\end{document} are assumed to follow a negative binomial distribution, they can be expressed as in Equation (3):

(3) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {X}_i\sim \mathrm{Negative}\ \mathrm{Binomial}\left({r}_i,p\right),{r}_i={r}_0\frac{c_i}{2}, \end{equation*}\end{document}

where the parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${r}_i$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p$\end{document} can be expressed as the mean \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $E\left[{X}_i\right]={\mu}_i$\end{document} and the variance \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $Var\left[{X}_i\right]={\sigma}_i^2$\end{document} of the negative binomial distribution by Equation (4):

(4) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {\mu}_i=\frac{r_ip}{1-p},{\sigma}_i^2=\frac{\mu_i^2}{r_i}+{\mu}_i. \end{equation*}\end{document}

Similar to the Poisson distribution, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mu}_i$\end{document} denotes the average RD of the normal genomic region; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\sigma}_i^2$\end{document} represents the data variance, which increases with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${\mu}_i$\end{document} in a quadratic form. The parameters of the negative binomial distribution offer more flexibility compared to those of the Poisson distribution, allowing it to have a larger variance, which better fit the observed sequencing data.

In general, the simulation data based on Poisson distribution are closer to the ideal sequencing situation, while the simulation data based on negative binomial distribution are closer to the real sequencing situation.

SIMULATION STUDIES

In order to compare the performance of the two types of core CNV detection algorithms comprehensively and accurately, we performed simulation experiments in three aspects, including different sequencing coverage, different CNV segment lengths and complex CNV cases. Similar to Yang et al. [33], we introduced Precision and Recall as the performance evaluation criteria. In addition, we evaluated the time and space complexity of the two algorithms.

Performance evaluation criteria

We introduce Precision and Recall to evaluate the segmentation algorithm performance, where Precision measures the ability of CNVs predicted correctly, and Recall measures the ability of the algorithm to identify all true CNVs. Specifically, Precision is the ratio of the number of correct CNVs detected to the total number of all CNVs detected, and Recall is the ratio of the number of correct CNVs detected to the total number of all true CNVs.

To further ensure the reliability of the results, we repeated the simulation experiment 100 times for each configuration, calculated Precision and Recall separately and the mean and SD for all results.

Comparison under different sequencing depths

Sequencing depth is the ratio of the total number of sequenced bp to the number of whole-genome bp, which is one of the important indicators to evaluate the sequencing quality [13]. In general, high sequencing depth means better data quality and lower error rates. However, high sequencing depth also entails heavy budgets and data processing cost. Therefore, decision of the sequencing depth needs to balance between the cost and data quality. In this simulation experiment, we generate data based on Poisson and negative binomial distributions and simulate different sequencing depths by setting different baseline RDs to compare the performance of the two segmentation algorithms at different sequencing depths. We set the RD signal length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N$\end{document} as 10 000 and generated CNV segments randomly. The CNV events were set as two groups, i.e. amplification and deletion (both heterozygous), and the number of CNV in each group was three, and the lengths were set as 30. The baseline RD was set to 6, 10 and 30, indicating sequencing depths of 3×, 5× and 15×, respectively. As for the HMM algorithm, we modify the number of hidden states to three (normal, amplification and deletion) based on the simulation data. The performance results of the experiment are shown in Table 2.

Table 2 Performance comparison with different sequencing depths

Depth	CBS	HMM	
Precision/%	Recall/%	Precision/%	Recall/%	
6	79.1 (32.0)	41.9 (20.9)	69.2 (12.1)	80.8 (10.6)	
9.8 (28.0)	3.6 (13.0)	63.9 (37.8)	55.9 (20.4)	
10	95.6 (4.4)	78.5 (15.0)	88.4 (6.6)	96.1 (4.2)	
89.6 (21.9)	41.0 (14.7)	84.2 (11.8)	77.1 (11.6)	
30	98.8 (1.3)	99.1 (0.8)	94.7 (3.8)	99.7 (0.6)	
97.3 (2.3)	97.1 (3.5)	92.2 (5.1)	98.9 (1.3)	
Precision and Recall are the average measurements based on 100 replications with the standard errors given in the parentheses. At each sequencing depth, the first line is the result of Poisson distribution, and the second line is the result of NB distribution.

From Table 2, it can be seen that the values of Precision and Recall increase with the increase of sequencing depth. When the sequencing depth increases to 15×, the algorithm segmentation performance gradually reaches a plateau. It is well confirmed that the data quality and algorithm segmentation performance increase with increasing sequencing depth, and this also indicates that the algorithm segmentation performance is affected at low sequencing depth. Besides, we can find that the segmentation metrics of the Poisson distribution outperform the negative binomial distribution, especially at low sequencing depths. This is because the simulated data based on the negative binomial distribution are more volatile than the Poisson distribution and its data are more discrete, making it challenging for the algorithm to perform the segmentation step.

Overall, the Precision of CBS is higher than that of HMM, which indicates that CBS detects more accurate CNV segments. In contrast, the Recall of HMM is consistently better than that of CBS, which implies that HMM detects more comprehensive CNV segments. However, CBS barely detected the variant segments when the data are simulated based on the negative binomial distribution and the sequencing depth is 3×. And, the SDs of Precision and Recall of both algorithms were high at the sequencing depth of 3×, indicating the increase of outliers. So, we further analyzed the data with a sequencing depth of 3×, as shown in Figure 1.

Figure 1 Algorithm performance metrics at sequencing depth of 3×. (A) Poisson distribution performance metrics, (B) Negative binomial distribution performance metrics.

Figure 1 shows the data distribution of the performance metrics of the two algorithms for 100 sets of simulation experiments at a sequencing depth of 3×. From Figure 1, it can be found that under simulated data with negative binomial distribution (i.e. increased data fluctuation), the recall of CBS decreases sharply and almost fails to detect CNV segments, which also affects the precision of CBS. Although the recall of HMM also decreases, it is more robust than CBS. This suggests that the segmentation performance of CBS is significantly degraded by its own recall at low sequencing depths, while HMM is more robust. Figures S2 and S3 in the Supplement shows the example figure of segmentation based on Poisson distribution for the two algorithms. The original data, CBS segmentation result and HMM segmentation result are shown from left to right in the figure.

Comparison under different CNV segment lengths

CNV segment length is an important measure of genomic variation and can be used to find the genomic variants associated with disease [65]. In general, longer CNV lengths may be associated with more severe phenotypes, but some studies have also shown that short CNV segments may likewise affect the gene expression or protein function, leading to the development and progression of some diseases [66]. For CNV detection, the length of the CNV segments is also an important factor to consider. In general, shorter CNV segments are more likely to be masked by noise or other bias, making them difficult to detect, while longer CNV segments are usually more stable and their detection is easier, providing more reliable CNV detection results. In this simulation experiment, we also simulated two CNV cases: amplification and deletion. Each contains three CNV segments with lengths from 1 to 50 (×1000 bp), respectively, and the sequencing depth is set to 15×; i.e. baseline RD is 30, while the rest of the parameters are set as in the above experiment. This simulation allows us to compare the performance of the two core algorithms in segmenting CNV segments of different lengths. The performance results of the experiment are shown in Figure 2. Afterward, we extracted results with lengths of 5, 10 and 30 (×1000 bp) for performance comparison, as shown in Table 3.

Figure 2 Algorithm performance metrics for variant segment length of 1–50. From left to right are precision and recall. (A) and (B) indicate Poisson distribution and NB distribution.

Table 3 Performance comparison with different segment lengths

Length/×1000 bp	CBS	HMM	
Precision/%	Recall/%	Precision/%	Recall/%	
5	96.0 (6.6)	52.5 (13.5)	9.0 (16.0)	6.7 (11.8)	
75.7 (39.2)	20.3 (14.8)	7.2 (16.6)	3.5 (7.2)	
10	96.2 (3.6)	93.9 (7.4)	81.9 (9.5)	95.7 (7.3)	
93.3 (7.0)	58.1 (16.0)	74.6 (11.9)	76.1 (16.1)	
30	98.8 (1.3)	99.1 (0.8)	94.7 (3.8)	99.7 (0.6)	
97.3 (2.3)	97.1 (3.5)	92.2 (5.1)	98.9 (1.3)	
Precision and Recall are the average measurements based on 100 replications with the standard errors given in the parentheses. At each segment length, the first line is the result of Poisson distribution, and the second line is the result of NB distribution.

From Table 3, we can clearly see that the values of Precision and Recall increase as the length of the variant segments increases. The algorithm segmentation performance gradually reaches a plateau when the length of variant segments increases to 30. This observation agrees with the actual situation because CNVs with small segments tend to be easily overwhelmed by noise.

The Precision and Recall of CBS are significantly higher than those of HMM in the face of tiny variant segments, i.e. when the length of the variant segment is 5. This indicates that CBS is better than HMM in detecting tiny segments. To further explore the performance of CBS and HMM in tiny segment detection, we investigated their data distribution.

From Figure 3, it can be found that when faced with tiny segments, HMM shows a large number of missed detections in both Poisson distribution and negative binomial distribution-based simulation data, indicating the poor performance of HMM in detecting such extreme conditions as tiny segments. This phenomenon improves rapidly with the increase of the variant segment lengths. On the other hand, CBS performs well in the case of detecting tiny fragments.

Figure 3 Algorithm performance metrics for a variant segment length of 5. (A) Poisson distribution performance metrics, (B) Negative binomial distribution performance metrics.

In addition, the segmentation based on Poisson distribution is still better than the segmentation based on NB distribution; especially when the CNV segment length is 5, this phenomenon is most obvious. This is consistent with the conclusions in the experiments in Comparison under different sequencing depths section. The specific segmentation can be seen in Figures S4 and S5 in the Supplement.

Comparison under complex CNV cases

CNVs can be mainly classified into two types: CN amplification and CN deletion. According to the degree of variation, they can also be classified into heterozygous variants and pure variants. Take the human genome as an example, where it is diploid; if one of the two homologous chromosome sequences are different from the sequence of the corresponding region of the reference genome, the variation is called heterozygous variation; if both homologous chromosome sequences are different from the sequence of the corresponding region of the reference genome, the variation is called homozygous (or pure) variation. This makes the situation of CNVs quite complicated, which poses a great challenge to accurately segment the CNV regions.

In this simulation experiment, we simulated various complex CNVs, including pure amplification, heterozygous amplification, pure deletion and heterozygous deletion, to evaluate the performance of the two core algorithms under complex CNVs. Based on the previous experimental results, each CNV case consists of two segments, and the segment length of each CNV segment is 30 and 50, respectively, and the sequencing depth is set to 15×. The rest of the parameters are set as in the previous model. Considering that the number of CNV species is changed, the number of hidden states of HMM is adjusted to 5 in this experiment according to the simulation data, and the transition matrix of HMM is rewritten. In addition, due to the limitations of sequencing technology, low sequencing coverage and sequence matching errors, some regions of the gene sequences may be missing in the actual sequencing situation. Therefore, we simulated the loss of gene sequencing data by randomly inserting two empty sequences of length 200 from the simulated data. The performance results of the experiment are shown in Table 4.

Table 4 Performance comparison under complex CNVs

	CBS	HMM	
	Precision/%	Recall/%	Precision/%	Recall/%	
No	98.7 (6.4)	99.2 (2.7)	97.4 (2.8)	99.5 (0.5)	
missing	96.6 (9.9)	98.9 (0.9)	96.7 (3.0)	99.3 (0.7)	
Missing	91.1 (18.1)	91.8 (14.5)	96.8 (12.0)	92.7 (13.4)	
89.2 (22.3)	90.5 (20.7)	93.8 (15.8)	92.5 (16.4)	
Precision and Recall are the average measurements based on 100 replications with the standard errors given in the parentheses. In each situation, the first line is the result of Poisson distribution, and the second line is the result of NB distribution.

As can be seen in Table 4, when the sequencing depth and variant segment length are moderate, increasing the CNV types and CN has minimal effect on the segmentation performance of the algorithm. However, by inserting the empty sequence, i.e. when the phenomenon of missing data occurs, the segmentation performance of both algorithms degrades. The CBS performance degrades significantly, and the SD of CBS Precision and Recall increases rapidly, indicating an increase in data outliers. This is because when the empty sequence is inserted, the simulation data are divided into multiple segments. When CBS performs the segmentation step, the simulated data are considered as several independent subsets at the interval of empty sequences, resulting in shorter segment lengths and posing challenges for segmentation. On the contrary, HMM treats all simulated data as a whole and performs an overall segmentation. This gives HMM an advantage over CBS in the case of complex and discontinuous sequencing data. In real sequencing scenarios, sequencing data are often discontinuous. The specific segmentation can be seen in Figures S6 and S7 in the Supplement.

Comparison of computational resources

We also analyzed the computation time and memory consumption of both algorithms when processing complex CNVs. The lengths of the test data are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${10}^3$\end{document},\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${10}^4$\end{document},\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${10}^5$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${10}^6$\end{document}, and the rest of the parameters are the same as those mentioned in the Comparison under different CNV segment lengths section. Similarly, we inserted empty sequences in the simulation data to increase the complexity of the data. The test environment is windows 10, R version x64 4.1.3, CPU is intel i9-10900F and the memory is 32G. We use the memory.size function in R to calculate the memory usage of the two algorithms, and we use the sys.time function to calculate the calculation time of the two algorithms. The results are shown in the Figure 4.

Figure 4 Calculation time and memory usage of two algorithms under different signal lengths.

As can be seen in Figure 4, the computation time of CBS is always longer than that of HMM, especially in the case of small-scale data, such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N={10}^3$\end{document}, where the computation time of CBS exceeds that of HMM by two orders of magnitude. This difference decreases with the increasing RD data length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N$\end{document}. This is because the computation time of HMM grows almost linearly with the RD data length compared to the computation time of CBS, which grows slowly with the RD data length. However, the computation time of CBS is always longer than that of HMM, and when the signal length is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${10}^6$\end{document}, the computation time of CBS still exceeds two times that of HMM. This means that CBS will take far more time than HMM when processing a large genome with high resolution.

In terms of memory usage, the two algorithms are not very different, with CBS consuming less memory compared to HMM. The memory usage of both algorithms increases with the increase of signal length, and the memory usage increases faster when the signal length is larger than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${10}^4$\end{document}. Experimental testing indicates that different distributions of the data have negligible impacts on the time and space cost of the algorithms.

DISCUSSION

In this paper, we comprehensively analyzed and summarized the segmentation algorithms employed by CN variation detection tools based on high-throughput sequencing over the last decade. We focused on extracting the core segmentation algorithms, namely, CBS and HMM, and conducted a systematic and comprehensive comparison of their performance. Meanwhile, through data parametric modeling simulation, we were able to eliminate the errors as associated with data alignment and GC content bias, which makes the comparison of the algorithm performance more equitable. These findings can provide important an guidance and reference for the research and development of new CNV detection tools.

The experimental results illustrate that under most ideal conditions, where only the heterozygous variation is included and the CN depth and variation segmentation length are distinct, CBS has a higher detection accuracy, while HMM has a higher sensitivity. However, in more practical cases (complex CNV), HMM outperforms the CBS, providing more accurate and reliable detections that are less affected by breakpoints. It is worth mentioning that the computation time of HMM is much smaller than that of CBS. When performing large-scale sequencing, HMM has a significant time-saving advantage over CBS. In addition, compared to CBS, HMM is more flexible in that it can be modified accordingly to meet specific data requirements to achieve better results. Therefore, we recommend researchers to use HMM in selecting segmentation algorithms in CNV detection. Nevertheless, CBS has its own special application scenario in which the variant to be detected is small.

Last, there is still room for improvement. (1) In our experiments, the method we used to approximate the baseline depth by using sample-averaged RD is only applicable to whole-genome sequencing, but which is not straightforwardly applicable to exome and other targeted sequencing. (2) CN amplification and deletion usually have different magnitudes, as amplification can vary from 3 to tens, while deletion is limited to no more than 2. This unbalanced in signal feature can potentially introduce bias into the statistical model and lead to a high false discovery rate, which needs attentions in future works. (3) We used the Baum–Welch algorithm to tune the model parameters of the HMM, but since the Baum–Welch algorithm uses the EM algorithm, local optimal solution is usually encountered, and global optimality is not guaranteed. (4) We compared the performance of two core CNV detection algorithms using real sequencing data from the 1000 genomes project [67] and the database of genomic variants [68] as the ground truth. However, since real sequencing data are complicated, ground-true CNVs are insufficient for most samples, e.g. only deletion without amplification, and it is hard to comprehensively compare and discuss the two segmentation algorithms based on those limited samples. Therefore, we mainly used simulated sequencing data for performance evaluation. To complement our study, an example of real data segmentation is shown in the Supplementary Fig. S1. More real examples are needed.

Key Points

For ideal conditions, CBS demonstrates high precision, while HMM exhibits a high recall rate.

For real conditions, HMM is advantageous at low sequencing depths, while CBS is more competitive in detecting variants with small segments.

In the case of complex CNVs close to real sequencing, HMM is more robust compared to CBS.

HMM is more time-efficient compared with CBS when processing large-scale sequencing data, while their memory usage is approximately equal.

Supplementary Material

supplement(10_bbae022_26)_bbae022 Click here for additional data file.

FUNDING

This work was supported by the National Natural Science Foundation of China (61771381).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author Biographies

Yibo Zhang received the BS degree in biomedical engineering from the Beijing Institute of Technology in 2022. He is currently pursuing the master's degree from the School of Life Science and Technology, Xi'an Jiaotong University, Xi'an, Shannxi, China. His research interests include signal processing and bioinformatics.

Wenyu Liu received the BS degree in biomedical engineering from the Central South University, in 2020, and the MS degree in biomedical engineering from the Xi’an Jiaotong University, Xi'an, Shannxi, China, in 2023. His research interests include signal processing and bioinformatics.

Junbo Duan received the BS degree in information engineering and the MS degree in communication and information system from the Xi'an Jiaotong University, Xi'an, China, in 2004 and 2007, respectively, and the PhD degree in signal processing from the Université Henry Poincaré, Nancy, France, in 2010. After graduation, he was a postdoctoral fellow at the Department of Biomedical Engineering and Bio- statistics and Bioinformatics, Tulane University, USA, until 2013. He is currently an associate professor at the Department of Biomedical Engineering, Xi'an Jiaotong University.
==== Refs
References

1. Mills RE , WalterK, StewartC, et al.  Mapping copy number variation by population-scale genome sequencing. Nature  2011;470 (7332 ):59–65.21293372
2. Alkan C , CoeBP, EichlerEE. Genome structural variation discovery and genotyping. Nat Rev Genet  2011;12 (5 ):363–76.21358748
3. Fellermann K , StangeDE, SchaeffelerE, et al.  A chromosome 8 gene-cluster polymorphism with low human beta-defensin 2 gene copy number predisposes to Crohn disease of the colon. Am J Hum Genet  2006;79 (3 ):439–48.16909382
4. Zhang F , GuW, HurlesME, LupskiJR. Copy number variation in human health, disease, and evolution. Annu Rev Genomics Hum Genet  2009;10 :451–81.19715442
5. Wheeler DA , SrinivasanM, EgholmM, et al.  The complete genome of an individual by massively parallel DNA sequencing. Nature  2008;452 (7189 ):872–6.18421352
6. Bentley DR , BalasubramanianS, SwerdlowHP, et al.  Accurate whole human genome sequencing using reversible terminator chemistry. Nature  2008;456 (7218 ):53–9.18987734
7. McKernan KJ , PeckhamHE, CostaGL, et al.  Sequence and structural variation in a human genome uncovered by short-read, massively parallel ligation sequencing using two-base encoding. Genome Res  2009;19 (9 ):1527–41.19546169
8. Korbel JO , UrbanAE, AffourtitJP, et al.  Paired-end mapping reveals extensive structural variation in the human genome. Science  2007;318 (5849 ):420–6.17901297
9. Singh AK , OlsenMF, LavikLAS, et al.  Detecting copy number variation in next generation sequencing data from diagnostic gene panels. BMC Med Genomics  2021;14 (1 ):214.34465341
10. Campbell PJ , StephensPJ, PleasanceED, et al.  Identification of somatically acquired rearrangements in cancer using genome-wide massively parallel paired-end sequencing. Nat Genet  2008;40 (6 ):722–9.18438408
11. Kidd JM , CooperGM, DonahueWF, et al.  Mapping and sequencing of structural variation from eight human genomes. Nature  2008;453 (7191 ):56–64.18451855
12. Zhao M , WangQ, WangQ, et al.  Computational tools for copy number variation (CNV) detection using next-generation sequencing data: features and perspectives. BMC Bioinform  2013;14 (11 ):1–16.
13. Zare F , DowM, MonteleoneN, et al.  An evaluation of copy number variation detection tools for cancer using whole exome sequencing data. BMC Bioinform  2017;18 (1 ):1–13.
14. Quenez O , CassinariK, CoutantS, et al.  Detection of copy-number variations from NGS data using read depth information: a diagnostic performance evaluation. Eur J Hum Genet  2021;29 (1 ):99–109.32591635
15. Chiang DY , GetzG, JaffeDB, et al.  High-resolution mapping of copy-number alterations with massively parallel sequencing. Nat Methods  2009;6 (1 ):99–103.19043412
16. Miller CA , HamptonO, CoarfaC, MilosavljevicA. ReadDepth: a parallel R package for detecting copy number alterations from short sequencing reads. PloS One  2011;6 (1 ):1–7.
17. Gusnanto A , WoodHM, PawitanY, et al.  Correcting for cancer genome size and tumour cell content enables better estimation of copy number alterations from next-generation sequence data. Bioinformatics  2012;28 (1 ):40–7.22039209
18. Klambauer G , SchwarzbauerK, MayrA, et al.  cn. MOPS: mixture of Poissons for discovering copy number variations in next-generation sequencing data with a low false discovery rate. Nucleic Acids Res  2012;40 (9 ):e69–9.22302147
19. Nguyen HT , MerrimanTR, BlackMA. The CNVrd2 package: measurement of copy number at complex loci using high-throughput sequencing data. Front Genet  2014;5 :248.25136349
20. Talevich E , ShainAH, BottonT, BastianBC. CNVkit: genome-wide copy number detection and visualization from targeted DNA sequencing. PLoS Comput Biol  2016;12 (4 ):e1004873.27100738
21. Roller E , IvakhnoS, LeeS, et al.  Canvas: versatile and scalable detection of copy number variants. Bioinformatics  2016;32 (15 ):2375–7.27153601
22. Dharanipragada P , VogetiS, ParekhN. iCopyDAV: integrated platform for copy number variations—detection, annotation and visualization. PloS One  2018;13 (4 ):e0195334.29621297
23. Yuan X , LiJ, BaiJ, et al.  A local outlier factor-based detection of copy number variations from NGS data. IEEE/ACM Trans Comput Biol Bioinform  2019;18 (5 ):811–1820.
24. Demidov G , OssowskiS. ClinCNV: novel method for allele-specific somatic copy-number alterations detection. BioRxiv. 2019;837971.
25. Dennis J , WalkerL, TyrerJ, et al.  Detecting rare copy number variants from Illumina genotyping arrays with the CamCNV pipeline: segmentation of z-scores improves detection and reliability. Genet Epidemiol  2021;45 (3 ):237–48.33020983
26. Xie K , LiuK, AlviHAK, et al.  IhybCNV: an intra-hybrid approach for CNV detection from next-generation sequencing data. Digit Signal Process  2022;121 :103304.
27. Ivakhno S , RoyceT, CoxAJ, et al.  CNAseg—a novel framework for identification of copy number changes in cancer from second-generation sequencing data. Bioinformatics  2010;26 (24 ):3051–8.20966003
28. Magi A , BenelliM, YoonS, et al.  Detecting common copy number variants in high-throughput sequencing data by using JointSLM algorithm. Nucleic Acids Res  2011;39 (10 ):e65.21321017
29. Szatkiewicz JP , WangWB, SullivanPF, et al.  Improving detection of copy-number variation by simultaneous bias correction and read-depth segmentation. Nucleic Acids Res  2013;41 (3 ):1519–32.23275535
30. Zhu M , NeedAC, HanY, et al.  Using ERDS to infer copy-number variants in high-coverage genomes. Am J Hum Genet  2012;91 (3 ):408–21.22939633
31. de  AraújoLL, WangK. PennCNV in whole-genome sequencing data. BMC Bioinform  2017;18 (11 ):49–56.
32. Yu Z , LiA, WangM. CLImAT-HET: detecting subclonal copy number alterations and loss of heterozygosity in heterogeneous tumor samples from whole-genome sequencing data. BMC Med Genomics  2017;10 (1 ):1–11.28057009
33. Yang H , ZhuD. Improved detection algorithm for copy number variations based on hidden Markov model. Multimed Tools Appl  2020;79 :9237–53.
34. Lee WP , ZhuQ, YangX, et al.  JAX-CNV: a whole genome sequencing-based algorithm for copy number detection at clinical grade level. Genom Proteom Bioinform  2022;20 (6 ):1197–206.
35. Melivoia R , YassineZ, JennyM, et al.  CoverageMaster: comprehensive CNV detection and visualization from NGS short reads for genetic medicine applications. Brief Bioinform  2022;23 (2 ):1–8.
36. Hui S , NielsenR. SCONCE: a method for profiling copy number alterations in cancer evolution using single-cell whole genome sequencing. Bioinformatics  2022;38 (7 ):1801–8.35080614
37. Xie C , TammiMT. CNV-seq, a new method to detect copy number variation using high-throughput sequencing. BMC Bioinform  2009;10 :1–9.
38. Yoon S , XuanZ, MakarovV, et al.  Sensitive and accurate detection of copy number variants using read depth of coverage. Genome Res  2009;19 (9 ):1586–92.19657104
39. Alkan C , KiddJM, Marques-BonetT, et al.  Personalized copy number and segmental duplication maps using next-generation sequencing. Nat Genet  2009;41 (10 ):1061–7.19718026
40. Abyzov A , UrbanAE, SnyderM, GersteinM. CNVnator: an approach to discover, genotype, and characterize typical and atypical CNVs from family and population genome sequencing. Genome Res  2011;21 (6 ):974–84.21324876
41. Xi R , HadjipanayisAG, LuquetteLJ, et al.  Copy number variation detection in whole-genome sequencing data using the Bayesian information criterion. Proc Natl Acad Sci  2011;108 (46 ):E1128–36.22065754
42. Boeva V , PopovaT, BleakleyK, et al.  Control-FREEC: a tool for assessing copy number and allelic content using next-generation sequencing data. Bioinformatics  2012;28 (3 ):423–5.22155870
43. Vardhanabhuti S , JengXJ, WuY, LiH. Parametric modeling of whole-genome sequencing data for CNV identification. Biostatistics  2014;15 (3 ):427–41.24478395
44. Smith SD , KawashJK, GrigorievA. GROM-RD: resolving genomic biases to improve read depth detection of copy number variants. PeerJ  2015;3 :e836.25802807
45. Onsongo G , BaughnLB, BowerM, et al.  CNV-RF is a random forest–based copy number variation detection method using next-generation sequencing. J Mol Diagn  2016;18 (6 ):872–81.27597741
46. Magi A , PippucciT, SidoreC. XCAVATOR: accurate detection and genotyping of copy number variants from second and third generation whole-genome sequencing experiments. BMC Genom  2017;18 (1 ):1–15.
47. Yuan X , ZhangJ, YangL, et al.  Detection of significant copy number variations from multiple samples in next-generation sequencing data. IEEE Trans Nanobioscience  2017;17 (1 ):12–20.
48. Yuan X , YuJ, XiJ, et al.  CNV_IFTV: an isolation forest and total variation-based detection of CNVs from short-read sequencing data. IEEE/ACM Trans Comput Biol Bioinform  2019;18 (2 ):539–49.
49. Li Y , ZhangJ, YuanX, LiJ. dpGMM: a Dirichlet process Gaussian mixture model for copy number variation detection in low-coverage whole-genome sequencing data. IEEE Access  2020;8 :27973–85.
50. Cabello-Aguilar S , VendrellJA, van GoethemC, et al.  ifCNV: a novel isolation-forest-based package to detect copy-number variations from various targeted NGS datasets. Mol Ther Nucleic Acids  2022;30 :174–83.36250203
51. Tony Cai T , Jessie JengX, LiH. Robust detection and identification of sparse segments in ultrahigh dimensional data analysis. J R Stat Soc Series B Stat Methodology  2012;74 (5 ):773–97.
52. Olshen AB , VenkatramanES, LucitoR, WiglerM. Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics  2004;5 (4 ):557–72.15475419
53. Wang K , LiM, HadleyD, et al.  PennCNV: an integrated hidden Markov model designed for high-resolution copy number variation detection in whole-genome SNP genotyping data. Genome Res  2007;17 (11 ):1665–74.17921354
54. Babadi M , FuJM, LeeSK, et al.  GATK-gCNV enables the discovery of rare copy number variants from exome sequencing data. Nat Genet  2023;55 (9 ):1589–97.37604963
55. Viailly PJ , SaterV, ViennotM, et al.  Improving high-resolution copy number variation analysis from next generation sequencing using unique molecular identifiers. BMC Bioinform  2021;22 (1 ):1–15.
56. Venkatraman ES , OlshenAB. A faster circular binary segmentation algorithm for the analysis of array CGH data. Bioinformatics  2007;23 (6 ):657–63.17234643
57. Baum LE , PetrieT. Statistical inference for probabilistic functions of finite state Markov chains. Ann Math Stat  1966;37 (6 ):1554–63.
58. Baum LE . An inequality and associated maximization technique in statistical estimation for probabilistic functions of Markov processes. Inequalities  1972;3 (1 ):1–8.
59. Viterbi A . Error bounds for convolutional codes and an asymptotically optimum decoding algorithm. IEEE Trans Inf Theory  1967;13 (2 ):260–9.
60. Ha G , RothA, LaiD, et al.  Integrative analysis of genome-wide loss of heterozygosity and monoallelic expression at nucleotide resolution reveals disrupted pathways in triple-negative breast cancer. Genome Res  2012;22 (10 ):1995–2007.22637570
61. Pan B , KuskoR, XiaoW, et al.  Similarities and differences between variants called with human reference genome HG19 or HG38. BMC Bioinform  2019;20 (2 ):17–29.
62. Kumps C , Van RoyN, HeyrmanL, et al.  Multiplex amplicon quantification (MAQ), a fast and efficient method for the simultaneous detection of copy number alterations in neuroblastoma. BMC Genom  2010;11 (1 ):1–8.
63. Langmead B , SalzbergSL. Fast gapped-read alignment with Bowtie 2. Nat Methods  2012;9 (4 ):357–9.22388286
64. Jeng XJ , CaiTT, LiH. Optimal sparse segment identification with application in copy number variation analysis. J Am Stat Assoc  2010;105 (491 ):1156–66.23543902
65. Cooper GM , CoeBP, GirirajanS, et al.  A copy number variation morbidity map of developmental delay. Nat Genet  2011;43 (9 ):838–46.21841781
66. Malhotra D , SebatJ. CNVs: harbingers of a rare variant revolution in psychiatric genetics. Cell  2012;148 (6 ):1223–41.22424231
67. 1000 Genomes Project Consortium, AbecasisGR, AltshulerD, et al.  A map of human genome variation from population-scale sequencing. Nature  2010;467 (7319 ):1061–73.20981092
68. MacDonald JR , ZimanR, YuenRKC, et al.  The database of genomic variants: a curated collection of structural variation in the human genome. Nucleic Acids Res  2014;42 (D1 ):D986–92.24174537
