SPCNV is an RD-based CNV detection method that is suitable for the detection of a single sample. Figure 1 shows the workflow of SPCNV in detail, which consists of the following main five steps: 1) The sequenced donor samples (Fastq) and reference genome (Fasta) are prepared for input; 2) The reads are aligned to the reference genome using BWA (Li and Durbin, 2010) to generate sequence alignment files (SAM), which are converted to BAM format using SAMtools (Li et al., 2009); 3) The data is preprocessed. This step mainly includes read count (RC) profile extraction with SAMtools, bin definition (Yuan et al., 2018), anomalous bin removal, obtaining the read depth (RD) profiles, GC bias calibration, noise removal, and the dimension transformation of the RD profiles; 4) The relative shortest path score is built and assigned for each RD; 5) Based on the score profile, a boxplot is utilized to predict CNVs. The SPCNV software is developed in R and Python languages. It can be downloaded from https://github.com/gj-123/SPCNV/releases and is easy to install and use after reading the user manual. In the following section, each step in the workflow of SPCNV is analyzed and discussed in detail.

The sequenced donor samples are aligned to the reference genome using the BWA tool, which generates sequence alignment files in the SAM format. The SAM files are further converted into binary sequence alignment files in BAM format using SAMtools. The read count (RC) profiles are then extracted with SAMtools from the BAM files. We define a sliding window procedure (bin) (Yuan et al., 2018), with which the RC profiles are divided continuously and are non-overlapping to generate the RD profiles. This process is described using Eq. 1. R D = R D 1 , R D 2 , R D 3 , ⋅ ⋅ ⋅ , R D n , (1) where R D n represents the RD value of the n-th bin, which is equal to the mean RC in a bin. Since the reference genome contains a large number of “N" positions, reads aligning to these positions will result in RCs equal to 0 that will be mistaken for a loss at that position. Therefore, if a bin contains “N" positions, we remove the bin from the genome sequence (Yuan et al., 2021a). Due to the complexity of the human genome, the distribution of GC-content is uneven, which can lead to the misidentification of copy number deletions. The GC content of each bin is calibrated using the median method (Yoon et al., 2009). Factors such as alignment errors and biases can cause the resulting sequencing data to be noisy. RD signal noise will seriously affect the detection accuracy, which is a key step in the detection of CNV. Here, the total variation model (Condat, 2013; Duan et al., 2013) is used to smooth and segment the RD profile to generate an RD segment (RDS) profile, which is represented by Eq. 2. R D S = R D S 1 , R D S 2 , R D S 3 , ⋅ ⋅ ⋅ , R D S n , (2) where R D S n represents the value of the n-th read depth segment, which is equal to the mean of all RDs contained in this segment. The RDS profile is converted to two-dimensional space to generate the R D S ′ profile (Liu et al., 2020), which is composed of the RDS ratio and differences between adjacent RDS ratios and is expressed by Eq. 3. R D S ′ = R D S X i , R D S Y i | i ∈ N , 1 ≤ i ≤ n , (3) where R D S X i represents the value of the i-th RD ratio, R D S Y i represents the difference between the i-th RD ratio and its adjacent RD ratios.

This transformation process can detect differences in RD from two perspectives. The first dimension can approximately reflect the copy number status corresponding to each RD from a global perspective. The second dimension can approximately reflect the difference between an RDS and its adjacent RDSs from a local perspective. By extracting two features of RD, the proposed method can more easily discover globally isolated and local small cluster variants. At the same time, this step also provides an effective data platform for constructing the relative shortest path score in the next section.

Based on R D S ′ profile, we construct a relative shortest path score (RSPS) to evaluate the degree of anomaly of each RDS. Here, we regard each element in R D S ′ as an object represented by o. The RSPS fully reflects the closeness between an object and its surrounding objects and is highly suitable for application in CNV detection scenarios. The RSPS of an object depends on the ratio between the object’s shortest path and the mean of its k nearest neighbors’ shortest paths. Some related basic concepts and definitions must be introduced before giving the definition of RSPS, which mainly include the k-distance of an object, the k-distance neighborhood of an object (Breunig et al., 2000), the shortest path set, the shortest path relation set, and the shortest path cost set.

Definition 1

The k-distance of an object o, which is defined using Eq. 4. k − d i s t o = d i s t o , o ′ , (4) where d i s t o , o ′ represents euclidean distance between object o and object o ′ ∈ R D S ′ \ o , o ′ indicates that the k-th object closest to o is sorted in ascending order, k − d i s t o represents the k-distance of an object o. Here, k is a positive integer.

Definition 2

The k-distance neighborhood of an object o is a collection of objects whose distance from o is less than or equal to k − d i s t o , which is defined using Eq. 5. N k − d i s t o = a | a ∈ R D S ′ \ o , d i s t o , a ≤ k − d i s t o , (5) where N k − d i s t o represents the k-distance neighborhood of an object o and a collection of objects, and the distance between each object in the collection and o is not greater than k − d i s t o .

The shortest path set (SPS) of an object o is composed of object o and the k-distance neighborhood of object o, which are connected to form a path with the shortest distance. The shortest path relation set (SPRS) of an object o is defined as the edge between two objects on the shortest path. The shortest path cost set (SPCS) of an object o is defined as the distance between two objects on the shortest path. Algorithm 1 describes the calculation process of SPS(o), SPRS(o), and SPCS(o) in detail.

The following example is used to clearly explain the SPS, SPRS, and SPCS calculation process of an object. As shown in Figure 2, there are a total of 10 objects {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, and their corresponding values are {(0.5, 0.5), (0.13, 0.78), (0.18, 0.13), (0.2, 0.4), (0.32, 0.38), (0.47, 0.44), (0.56, 0.6), (0.6, 0.56), (0.87, 0.5), and (0.91, 0.93)}. We calculate the SPS, SPRS, and SPCS of object one here, where the value of k is set to 5. Eq. 4 and 5 are employed to obtain the five nearest neighbors of object 1 (6, 7, 8, 5, 4). The distance from object one to object six is the smallest. According to Algorithm 1, SPS (1) is equal to {1, 6}, SPRS (1) is equal to {(1, 6)}, and SPC S(1) is equal to {0.07}. The distance between objects 7, 8, 5, 4, and objects 1, six is calculated to obtain a minimum distance, and the distance between objects one and seven is the smallest. Similarly, SPS (1) is equal to {1, 6, 7}, SPRS (1) is equal to {(1, 6), (1, 7)}, and SPCS (1) is equal to {0.07, 0.12}. The procedure ends when each neighbor of object one finds its closest object from SPS(1). Finally, SPS (1) is equal to {1, 6, 7, 8, 5, 4}, SPRS (1) is equal to {(1, 6), (1, 7), (7, 8), (6, 5), (5, 4)}, and SPCS (1) is equal to {0.07, 0.12, 0.06, 0.16, 0.12}. As shown in Figure 2, the red line segments form the final shortest path.

Algorithm 1

Calculate SPS, SPRS and SPCS of object o.

Based on the above definitions, the mean shortest path cost of an object o ( S P C m o ) is defined by Eq. 6. S P C m o = 1 k ∑ i = 1 k S P C S i o , (6) where S P C S i o represents the i-th element in S P C S o .If the value of S P C m o is larger, the distance between o and its k nearest neighbors is sparser; if it is smaller, they are closer together.

After estimating the mean shortest path cost of all objects, we further construct the relative shortest path score (RSPS) to measure the degree of deviation between object o and its k nearest neighbors, which is defined by Eq. 7. R S P S o = N k − d i s t o ⋅ S P C m o ∑ a ∈ N k − d i s t o S P C m a , (7) where S P C m o represents the mean shortest path cost of o, S P C m a represents the mean shortest path cost of one of its nearest neighbors, N k − d i s t o represents the number of elements in N k − d i s t o , R S P S o represents the ratio between S P C m o and the mean of mean shortest path cost of its k nearest neighbors. If the value of R S P S o is larger, the distance between o and its k nearest neighbors is sparser; if it is smaller, they are closer together. This means that the higher the RSPS of an object, the more likely the object is a CNV.

Although we have evaluated the relative shortest path score for each object, it is not yet possible to distinguish abnormal objects from normal objects. This step is critical in CNV detection, and a reasonable threshold selection strategy will significantly improve the accuracy of the detection results. Traditional threshold selection strategies mainly include: 1) Fitting a statistical model using the score profile to evaluate the significance of each object and using hypothesis testing to predict CNVs; 2) Selecting an empirical value as a threshold to identify CNVs. The limitation of the first strategy is that due to the bias and noise of the sequencing data, the actual distribution of the data and the fitted model are quite different, resulting in inaccurate detection results. The limitation of the second strategy is that a fixed threshold can effectively identify abnormal objects in some scenarios. However, the performance of the method may drop significantly in certain application scenarios. Considering the above two scenarios, we use a boxplot to determine thresholds based on the score files. This method does not require an assumption that the score profile obeys a certain distribution in advance and is able to dynamically determine thresholds according to different score files. Here, we use Eq. 8 to estimate a threshold for judging anomalies. τ = R S P S Q 3 + λ ⋅ R S P S Q 3 − R S P S Q 1 , (8) where R S P S Q 3 represents upper quartile of RSPS, R S P S Q 1 represents lower quartile of RSPS, λ represents multiple of interquartile range of RSPS, τ represents the maximum value of the inner limit of RSPS, which is used as the threshold. If an object’s RSPS is greater than the threshold τ , it is considered to be a CNV. After predicting CNVs, we further differentiate the types of CNV (gain and loss). If the RD of an object is greater than or equal to the average RD of all normal objects, it is regarded as a gain; if it is less, it is regarded as a loss.

IntSIM (Yuan et al., 2017) simulation software was adopted to generate the simulation data sets. Before using the software, the two key parameters of sample tumor purity (TP) and sequencing coverage (SC) were set from 0.2 to 0.8 and 5x, respectively. To ensure the reliability of the test results, 50 samples were generated under each set of configuration conditions and the average of which was used as the test result. There were six gains and eight losses embedded in each sample, whose lengths range from 10,000 to 50,000 bp.

Based on the simulated datasets, the performance of SPCNV and four other alignment methods (CNV-LOF, FREEC, CNVnator, and BIC-seq2) were tested by calculating their recall, precision, and F1-scores. Recall is defined as the number of correctly detected CNV events divided by the total number of simulated CNV events, which can be calculated by the ground truth file (Magi et al., 2013). Precision is defined as the number of correctly detected CNV events divided by the total number of detected CNV events (Magi et al., 2013). The F1-score is defined as the harmonic mean of recall and precision. The experimental results of each method are depicted in Figure 3, where the three-performance metrics (recall, precision, and F1-score) of each method are compared in the four simulation sample sets. According to the overall trend, the performance of the majority of methods improves with increasing tumor purity. For example, the recall of CNVnator is close to 0.2 when the tumor purity is equal to 0.2, and its recall exceeds 0.7 when the tumor purity is equal to 0.8. Correspondingly, its F1-score increases from 0.18 to 0.58. Among the five methods, SPCNV achieves the best F1-score in each dataset. BIC-seq2 obtains the lowest F1-score at a purity equal to 0.8, but its F1-scores are better than other three methods (CNV-LOF, FREEC, and CNVnator) at a purity equal to 0.6. The above situations indicate that BIC-seq2 can provide base level resolution and detect a large number of CNVs, but its precision is very low when detecting high-purity samples. BIC-seq2 is not suitable for the detection of high-purity samples. The F1-score of FREEC is better than other three methods (CNV-LOF, BIC-seq2, and CNVnator) when the tumor purity is equal to 0.8. When detecting low and medium purity samples, its performance is relatively balanced between recall and precision. When FREEC detects high-purity samples, its F1-score is superior to the other three comparison methods, but its recall is significantly higher than precision, which indicates that its performance in detecting high-purity samples is uneven. The F1-scores of CNV-LOF are better than other three methods (BIC-seq2, FREEC, and CNVnator) when the tumor purity is equal to 0.2 and 0.4, which indicates that it is suitable for detecting low and medium purity samples. Its advantage is to obtain high precision when detecting low and medium purity samples. In contrast, its recall rate is low. CNVnator obtains the lowest F1-scores at purities equal to 0.2, 0.4, and 0.6, because the precision of CNVnator is the lowest in all sample sets. CNVnator detects a large number of long CNVs, most of which are false-positive positions. In terms of recall, SPCNV achieves the best recall, except when the tumor purity is equal to 0.6. In terms of precision, SPCNV obtains the best precision in each sample set. Overall, the SPCNV gets the best trade-off between recall and precision in each simulation sample set.

Based on the above analysis and discussion, we further analyze the performance of each method in detecting gains and losses. Figure 4 details the performance of each method in detecting gains and losses in the four datasets. In general, the number of gains detected by each method is more than losses. The performance of SPCNV in detecting the gain and loss is relatively balanced compared to the other four methods. While CNV-LOF detects the most gains in each sample set, it obtains the least losses in the vast majority of cases. In most cases, FREEC detects far more gains than losses. CNVnator detects the least gains in each sample set and obtains fewer losses than most methods. BIC-seq2 detects the most losses in three sets of samples, indicating that it is suitable for detecting losses. In summary, SPCNV detects more gains and losses than most methods in each simulation sample set, which shows that its performance is robust in gain and loss detection.

As a supplement to the above experiments, we further analyzed the sensitivity of each method to detect CNV of different lengths. Figure 5 details the sensitivity of five methods to detect three CNV lengths (10 kb, 20 kb and 50 kb) in four sets of samples. The sensitivity is defined as the ratio between the number of correctly detected CNVs and the number of simulated CNVs. The proposed method achieves the best sensitivity when the TP is equal to 0.2 and the CNV length is equal to 50 kb, and ranks second when CNV sizes are equal to 10kb and 20 kb. BIC-seq2 and CNV-LOF get the best sensitivity when the TP is equal to 0.2 and CNV sizes are equal to 10 kb and 20 kb. CNV-LOF achieves lower sensitivity at CNV sizes equal to 10 kb and 50 kb than other three comparison methods (SPCNV, FREEC and BIC-seq2). BIC-seq2 achieves lower sensitivity at CNV sizes equal to 20 kb and 50 kb than SPCNV. The proposed method gets the best sensitivity when the TP is equal to 0.4 and CNV sizes are equal to 20 kb and 50 kb, and ranks second when CNV size is equal to 10 kb. CNV-LOF and BIC-seq2 get the best sensitivity when CNV size is equal to 10 kb, and achieves lower sensitivity at CNV sizes equal to 20 kb and 50 kb than the proposed method. The proposed method gets the best sensitivity when the TP is equal to 0.6 and CNV size is equal to 20 kb, and ranks second when CNV sizes are equal to 10 kb and 50 kb. At the same time, CNV-LOF and BIC-seq2 get the best sensitivity when CNV size is equal to 10 kb and 50 kb, respectively. CNV-LOF achieves lower sensitivity at CNV sizes equal to 20 kb and 50 kb than SPCNV. BIC-seq2 achieves lower sensitivity at CNV sizes equal to 10 kb and 20 kb than SPCNV. The proposed method gets the best sensitivity when the TP is equal to 0.8 and CNV size are equal to 20 kb and 50 kb, and ranks second when CNV size is equal to 10 kb. CNV-LOF gets the best sensitivity when CNV size is equal to 10 kb, and achieves the lowest sensitivity at CNV sizes equal to 20 kb and 50 kb. In general, the proposed method performs better than the other four comparison methods in detecting CNVS of different sizes.

In order to verify the performance of the proposed method in real application scenarios, we use six real data samples (NA12878, NA12891, NA12892, NA19238, NA19239, and NA19240) from the 1,000 Genomes Project, which can be downloaded for free from http://www.internationalgenome.org/. Some of the test results for these samples are recorded in the DGV, which is considered the ground truth for calculating the recall, precision, and F1-scores for each method. Similarly, we select the four methods (CNV-LOF, FREEC, CNVnator, and BIC-seq2) of the above experiments to compare with the proposed method. The experimental results are shown in Figure 6. SPCNV obtains the highest F1-scores among the five samples and ranks second in the NA19238 sample. CNV-LOF obtains the highest F1-score in the NA19238 sample and ranks second in F1-score among the other five samples. BIC-seq2 does not detect the correct CNVs in NA12878, NA12891, and NA12892, and its F1-scores rank third in NA19238, NA19239, and NA19240. The performance of FREEC and CNVnator is relatively close, and their F1-scores are ranked third and fourth in NA12878, NA12891, and NA12892 and in NA12878, NA12891, and NA12892, respectively. In terms of recall, FREEC achieves the best recall three times, CNVnator has the best recall twice, and CNV-LOF achieves the best recall once. In terms of precision, SPCNV obtains the best precision in five of the six samples, and CNV-LOF has the best precision once. Overall, SPCNV has obvious advantages over the other four methods in terms of precision and F1-score.