DTDHM: detection of tandem duplications based on hybrid methods using next-generation sequencing data

Workflow of DTDHM

DTDHM analyzes TDs without the need for control samples. The workflow of DTDHM is shown in Fig. 1 and consists of three main steps.

In the first step, RD and MQ are extracted from BAM files and preprocessed, accounting for missing values and the “N” positions in the reference genome. The sliding window splits the whole gene sequence into several continuous and non-overlapping bins. The RD and MQ are calculated for each bin and any bias caused by uneven GC content is corrected. The TV model is used to smooth the two signals and reduce noise. Then, the circular binary segmentation (CBS) algorithm is used to segment the bins into genome segments according to RD value fluctuations. Finally, RD and MQ are standardized.

In the second step, the candidate TD regions are obtained. RD and MQ are input into the DTDHM model as two features. The model calculates the outliers of each genome segment based on the KNN algorithm and detects some regions that may contain TDs as noise. The boxplot procedure sets an appropriate threshold for outliers. If the outlier in a region is greater than the threshold, it is considered a candidate TD region.

In the third step, the SR strategy is used to obtain split reads from BAM files. Then, CIGAR fields are used to infer the breakpoint position of the TDs. Finally, PEM is used to find and analyze discordant reads to refine the boundaries of TDs.

The DTDHM software is implemented primarily in Python. The following subsections describe each step of DTDHM in detail.

• 1) Dealing with missing values and “N” positions

• The reference genome usually consists of four bases: “A,” “T,” “G,” and “C.” The character “N” is used to indicate that the position of the base has not been determined and could be any of the four bases mentioned above. There are missing values and “N” positions in any current version of the reference genome. These problems lead to mapping errors, where the read count deviates from the true read count. To provide a more reasonable read count, DTDHM fills the missing position with zero and adopts the removal strategy for the “N” position, as explained by Liu et al. (2020) and Yuan et al. (2019).

• 2) Obtaining RD and MQ

• Due to the influence of noise from sequencing errors, mapping errors, and GC content bias, the RD signal does not follow a linear relationship with the TD region. Therefore, the variance of the RD signal is insufficient for detecting regions of variation in the genome sequence. MQ can indicate how reliably the read aligns to a position in the reference sequence and can be used in combination with other features for mutation detection (Zhao et al., 2020). For these reasons, DTDHM introduces MQ signals as one of the characteristics used for detecting TDs.

  The genome sequence is divided into contiguous non-overlapping bins using a sliding window of fixed width. Because the “N” position in the reference sequence does not indicate whether tandem duplication events have occurred in the region, bins containing “N” are filtered. Then, the read pair count (RC) of each position in the bin is counted. The average value of RC in each bin is taken as the read depth signal of this bin, denoted as RD. The Formula (1) for calculating RD is as follows: (1) $$R{D}_i={\displaystyle \frac{{\sum }_{j=1}^{L_b}R{C}_{ij}}{L_b},iϵ\left[1,N_r\right],}$$

• where $R{D}_i$ represents the RD value of the i-th bin; $R{C}_{ij}$ represents the RC value of the j-th position in the i-th bin; $L_b$ represents the length of the genome bins, such as 1,000 bp; and $N_r$ is the number of genomic bins generated from all regions.

  MQ value is the average Mq value at each bin position. Its calculation method is similar to RD, and is calculated using the following Formula (2): (2) $$M{Q}_i={\displaystyle \frac{{\sum }_{j=1}^{L_b}M{q}_{ij}}{L_b},iϵ\left[1,N_r\right],}$$

• where $M{Q}_i$ indicates the MQ value of the i-th bin, and $M{q}_{ij}$ indicates the Mq value at the j-th position in the i-th bin.

• 3) Correction of GC content bias

• NGS sequencing can lead to GC content bias during PCR amplification (Dohm et al., 2008). GC content bias affects the RD signal and then interferes with the prediction of tandem duplications. To eliminate the GC content bias, the GC fraction of each bin is calculated and then GC correction is performed on the RD signal of each bin (Yoon et al., 2009). The Formula (3) for GC correction is as follows: (3) $$R{D}_i^{\mathrm{\prime }}={\displaystyle \frac{\overline{RD}}{{\overline{RD}}_{GC}}\cdot R{D}_i,iϵ\left[1,N_r\right],}$$

• where $R{D}_i^{\mathrm{\prime }}$ and $R{D}_i$ represent the revised value and original value of the RD of the i-th bin, respectively; $\overline{RD}$ is the mean of all RD’s; and ${\overline{RD}}_{GC}$ represents the average RD that is similar to the GC score of the ith bin. Bins with similar GC scores are those bins with GC score differences less than or equal to 0.002.

• 4) Denoising using total variation

• The RD and MQ between neighboring bins have a natural correlation. However, due to sequencing and mapping errors, the RD and MQ of adjacent bins may fluctuate randomly, resulting in noise. Therefore, a regular term is added to keep RD and MQ smooth. Unlike simple denoising methods, the TV model can remove noise while preserving the boundary information of the bins. Because of this, DTDHM uses the TV model to smooth and de-noise RD and MQ. The TV model is expressed as the solution to the minimization problem, and the Formula (4) is as follows: (4) $$\underset{x\in R^{N_r}}{min}\left\{{\displaystyle \frac{1}{2}{\Vert y-x\Vert }_2^2+\lambda \sum _{i=1}^{N_r-1}\left|x_{i+1}-x_i\right|}\right\},$$

• where x represents the signal containing noise; y represents the signal obtained after smoothing; both x and y are vectors, i.e., $x={\left[x_1,\dots ,x_{N_r}\right]}^T$ and $y={\left[y_1,\dots ,y_{N_r}\right]}^T$; $N_r$ represents the total number of bins; T is the matrix transpose; the first term represents the fitting error between x and y; the second term represents the total variation penalty; and λ represents the penalty parameter that controls the tradeoff between the first and second terms. Because choosing parameter λ can be a problem (Condat, 2013; Vaiter et al., 2013), DTDHM allows the user to specify the value of parameter λ, which is recommended to be selected from the range [0.15, 0.30]; the default is 0.15 (Yuan et al., 2021).

• 5) Region segmentation and standardization processing

• Because the KNN algorithm is unsuitable for large sample sizes, DTDHM uses CBS to segment the bin-based genome into segments of different sizes, which reduces the computational complexity of KNN. RD and MQ are input into the KNN algorithm as two features. However, MQ and RD have a data imbalance problem because the value of MQ is much larger than the value of RD. When the KNN algorithm finds the nearest neighbor samples in the feature space, the contribution of MQ is greater than that of RD in the algorithm. To balance the influence of the two factors in the KNN algorithm, RD and MQ are normalized, as shown in the following Formula (5): (5) $$r^{\mathrm{\prime }}=\left(r-\mu \right)/\sigma ,$$

• where r indicates RD or MQ, $r^{\mathrm{\prime }}$ is the standardized sample, $\mu$ represents the mean value of all sample data, and $\sigma$ represents the standard deviation of all sample data.

Calling KNN to get candidate TDs

After preprocessing, the normalized RD and MQ are input into the KNN algorithm as two features. RD and MQ reflect magnitude of variation and mapping quality, respectively. The RD signal in one-dimensional space is transformed into a two-dimensional profile, D, combining RD and MQ (Yuan et al., 2019), as shown in Formula (6). The variation interval is analyzed from two perspectives to improve detection accuracy further. (6) $$D=\left\{\left(R{D}_i,M{Q}_i\right)|,iϵ\left[1,N_r\right]\right\}.$$

DTDHM adopts a KNN algorithm based on anomaly detection, which uses the basic principles of KNN but calculates outlier scores rather than predicting labels or values. The KNN anomaly detection algorithm is a distance-based method. Its core idea is to compute the distances between a genomic fragment point on a two-dimensional profile and its nearest k neighboring fragment points, evaluate the degree of anomaly of the sample through the size of the distance between the sample points, and judge whether it is a TD sample point. The characteristics of TDs cause the RD value and MQ value of the variant region to generally deviate from the value range of the normal region, and the k-nearest neighbor distance of a TD is usually much larger than that of a normal point. When using this rule to further calculate the outlier scores of genomic fragments, the larger the outlier scores of a genomic fragment, the more likely it is a TD. In this article, a genomic fragment is regarded as an object p.

In the KNN algorithm, in order to reduce the computational cost and speed up the search of k-nearest neighbors, DTDHM uses k-dimensional tree (KDTree), a high-dimensional index tree data structure. Using KDTree search, the k-nearest neighbor of any object in data space has a computational cost of $O\left(N\lg N\right)$. For the calculation of distance, DTDHM uses Euclidean distance. There are three common outlier scoring schemes: maximum, median, and average. The maximum and median outlier scoring schemes ignore the information from other nearest neighbors and exhibit instability. Because of this, DTDHM uses the more stable average outlier scoring scheme, which is the average distance between object p and all objects within its neighborhood.

When identifying tandem duplications, adjacent regions in the genome are positionally correlated (Yuan et al., 2017). After executing KNN and KDTree algorithms, the returned noise points can be regarded as a set of candidate TDs containing adjacent regions. Successive fragments are then merged to obtain candidate regions. However, because the KNN algorithm is sensitive to outliers, it can produce many false positives in the detection process, so the candidate regions then need to be filtered.

To accurately declare TDs based on outliers, a reasonable statistical model is needed to provide a cutoff point for outliers. The classical method produces an uneven distribution of sequencing reads, affecting the detection of low-amplitude TDs, so DTDHM uses the boxplot procedure to mark outliers (Sim & Chang, 2005) without any restrictions on the data (such as data following a certain distribution). This method uses interquartile range (IQR) to detect outliers. The interquartile range is the difference between the upper quartile (UQ) and the lower quartile (LQ) and contains half of all the data. Using parameter θ, the upper and lower values are defined as $UQ+\theta \cdot IQR$ and $LQ-\theta \cdot IQR$. Objects that exceed the upper or lower limit of the boxplot are considered outliers. In tandem duplications, a higher outlier on the boxplot indicates a higher probability of an outlier, so the objects above the upper limit of the boxplot are considered candidate regions.

For an object p, Table 1 describes the steps KNN takes to declare a candidate TD. The KNN algorithm takes RD and MQ as features. It selects any sample point in the feature space and searches for its k-nearest neighbor domains using KDTree. Then, the average outlier score scheme is used to calculate the average distance between the sample point and its k-nearest neighbors, or the outlier score of the sample point. These outlier scores are obtained for all sample points, then a statistical analysis is carried out using the boxplot procedure. When the sample points fall above the upper limit of the boxplot, they are classified as outliers. The candidate TDs are obtained by counting all outliers.

Filtering candidate regions

Candidate TD regions then have to be filtered for other types of structural variations, such as interspersed duplication. The SR strategy is used to identify and filter out the split reads that do not belong to the target tandem duplication. The SR strategy uses the local alignment algorithm to align the pair-end read pair with the reference genome and split the read segments that do not completely match the reference genome. Tandem and interspersed duplication can be distinguished by analyzing the characteristics and alignment of split reads. The CIGAR field in the BAM file indicates the alignment of the read segment with the reference genome, usually consisting of the characters “H,” “M,” and “S.” Split reads are determined by whether the CIGAR field contains an “S” cut symbol. In this article, CIGAR mainly analyzes two types of pre-alignment and post-alignment reads, as represented by Formula (7):

(7) $$CIGAR=\{\begin{array}{c}xMyS\\ ySxM\end{array}x,y\in N,x+y=R_l,$$where “M” means that the digit length in front of the symbol matches the reference position, and “S” means that it does not match and is soft cut. “xMyS” indicates pre-alignment; “ySxM” indicates post-alignment; and $R_l$ represents the length of the short sequencing read segment, which is generally set to 100 bp. For example, if the length of the pair-end reads segment is 100 bp, the CIGAR field in one of the reads is “20S80M,” indicating that the read is post-alignment, the first 20 bp of the read does not match the reference position, but the remaining 80 bp does match the reference position.

Split reads in the alignment results are then scanned to find the set of potential breakpoints and summarize the SR characteristics of tandem and interspersed duplications. Figure 2 shows the SR characteristics of tandem and interspersed duplications. In Fig. 2A, R1 and R2 are two reads that span the junction of the repeat sequence. When the reads are mapped to the reference genome, both R1 and R2 are multiply mapped to points A and B. Point A represents the starting position of the variation, while point B represents the ending position of the variation. R1 and R2 are post-alignment (ySxM) at point A and pre-alignment (xMyS) at point B. This indicates that split reads generated by tandem duplications can map multiply to two breakpoints, which are, respectively, the starting position and the ending position of the tandem duplication. Furthermore, the alignment pattern at the breakpoints is consistent, being post-alignment at the starting position and pre-alignment at the ending position. In Fig. 2B, R3 and R4 are two split reads generated by interspersed duplication. R3 is multiply mapped to points A and C, while R4 is mapped to points B and C. Point C represents the insertion site of the interspersed duplication, with two alignment patterns (xMyS + ySxM). This shows that split reads generated by interspersed duplication can map to multiple breakpoints, with all breakpoints other than the starting and ending position including both pre-alignment and post-alignment.

DTDHM uses the features of split reads to extract tandem duplications and filter out interspersed duplications from candidate regions. Post-alignment split reads are filtered out at the starting position of the candidate region, then other mapping positions of these reads are identified. If all split reads mapped to these other positions are pre-alignment and are located at the ending position of the candidate region, the candidate region is labeled as a tandem duplication; otherwise, it may be an interspersed duplication. Pre-alignment split reads are then filtered out at the ending position of the candidate region and the breakpoints of these reads are recorded. If all reads at the breakpoint position of this read are pre-alignment, and all reads at the other mapping positions of this read are post-alignment, the candidate region is labeled as a tandem duplication.

Refining boundaries

Because the KNN algorithm only uses RD and MQ signals, its detected candidate TD region boundary is inaccurate. To improve TD boundary accuracy to the bp level rather than the bin level, DTDHM optimizes the boundaries of TD regions based on the breakpoint positions of tandem duplications. There are many split reads in whole sequencing data and most split reads are useless for boundary detection, so DTDHM searches for split reads within the valid range. Figure 3 depicts an example of DTDHM exploring the effective split reads in a candidate TD region [a, b]. The length of m bins is used as the effective search step (S = m∙ $L_b$). Prior experience has indicated that most of the effective split reads gather in the range of two bins, and expanding that range may lead to noisy data, so m was set to 2. The starting position a and ending position b of the TD region detected are used as the search center. The detection range at both ends of the candidate TD region is 2S. Within this range, SRs containing a large amount of breakpoint information are extracted to correct the boundary of the candidate TD region [a, b]. Ideally, the boundary [a, b] can be precisely [A, B]. Using the SR strategy, the boundary of the TD region can be expressed by the following Formula (8):

(8) $$A=POS(Rpost(i)),B=POS(Rpre(j))+x_{Rpre(j)}-1,$$where A represents the exact starting position of the candidate TD region; B represents the exact ending position of the candidate TD region; Rpost(i) represents a post-aligned split read in the range [a – S, a + S]; Rpre(j) represents a pre-aligned split read in the range [b – S, b + S]; POS() indicates the position where the read aligns with the reference genome; and $x_{Rpre\left(j\right)}$ represents the length of the Rpre(j) that matches the reference genome. If R1 is post-alignment near position a, the exact starting position of the TD region can be determined by $A=POS\left(Rpost\left(1\right)\right)$. If R2 is pre-alignment near position b, the exact ending position of the TD region can be determined by $B=POS\left(Rpre\left(2\right)\right)+x_{Rpre\left(2\right)}-1$. Finally, the rough boundary [a, b] of the TD region is exactly [A, B].

However, SR can only refine a limited area at low coverage depth and tumor purity. After processing with SR, some boundaries of TD regions [a, b] remain rough, such as [a, B], [A, b], and [a, b]. Therefore, the PEM strategy is then used to further deal with the boundaries not covered by SR. Figure 4 shows an example of PEM in tandem duplication. This strategy focuses on pair-end reads spanning tandem duplication junctions that align to the reference sequence to form discordant read pairs. Among these discordant read pairs, the read (R2′) at the front of the alignment position is closest to the starting position of the TD region. Similarly, the read (R5′) at the back of the alignment position is closest to the ending position of the TD region. Using this rule, a more accurate TD region is obtained.

Figure 5 depicts an example of using PEM signals to refine candidate TD boundaries. Since SR signals have been previously used to refine the boundaries of TDs, only the reads without SR signals are included when processing PEM signals. First, the location of the extracted discordant read pair must be within the search radius S of the starting position (which may contain the true breakpoint information). Second, the absolute value of the insert size of discordant reads must be within a certain range of the length of the detected candidate TD region (the default range is 2S). DTDHM sorts discordant reads (such as R4′ and R6′) that meet these two conditions by the size of the alignment position. The alignment position of read R4′ is the most forward, and the alignment position of read R6′ is the most backward. The position of R5′, another pair-end read of R6′, is closest to the TD ending position. According to these rules, the alignment position of read R4′ can approximate the TD starting position, and the position of R5′ plus the length of the read segment ( $R_l$) can approximate the TD ending position. Using the pair-end mapping strategy, the boundary of the TD region can be expressed by the following Formula (9):

(9) $$A=POS(Rfirst(i)),B=POS(Pnex{t}_{Rlast(j)})+R_l,$$where $Rfirst\left(i\right)$ represents the read pair with the smallest comparison position among the discordant read pairs that meet two conditions; $Rlast\left(j\right)$ represents the read pair with the largest comparison position among the discordant read pairs that meet all conditions; and $Pnex{t}_{Rlast\left(j\right)}$ represents the other read pair corresponding to $Rlast\left(j\right)$ in paired-end sequencing; $R_l$ is defined in Formula (7).

Table 2 details the process of declaring TDs and the precise location of TD boundaries.

Simulation studies

In generating simulated data, Python was used to generate a FASTA file containing variation information. SinC (Pattnaik et al., 2013) and seqtk (https://GitHub.com/lh3/seqtk) software were combined to generate datasets of different levels of tumor purity and coverage depth. The SInC_readGen program in the SInC software can generate two types of FASTQ files based on FASTA files: one is normal and the other contains mutation information. The “−R” and “−C” parameters in the SInC_readGen program can specify read length and coverage depth. “−R” was set to 100. It is worth noting that due to the paired-end sequencing used in this study, merging a pair of read FASTQ files doubled the coverage depth of the merged FASTQ file; therefore, the value of “−C” should be half of the actual coverage depth. Next, the sub-sampling feature in Seqtk was used to extract data from both the normal and mutant FASTQ files to generate data with different levels of tumor purity. Then, the sequencing reads were mapped to the reference genome using the BWA-MEM (Li & Durbin, 2010) tool to generate a SAM file. Finally, SAMtools (Li et al., 2009) software was used to convert the files from SAM format to BAM format and sort the BAM files.

In the simulated dataset, tumor purity was set to 0.2, 0.4, and 0.6, and coverage depth was set to 4, 6, and 8X. In each simulation configuration, the size of tandem and interspersed duplications ranged from 10 to 50 kbp. Chromosome 21 in GRCh38 was selected as the reference genome. Each configuration generated 50 samples and their average sensitivity, precision, and F1-score were calculated to reduce the randomness of the experiment. DTDHM and three existing methods were tested on the simulated dataset. A true positive was recorded when the detection result covered more than half of a true TD area. The calculation methods for sensitivity, precision, and F1-score are shown in the following Formulas (10)–(12): (10) $$Sensitivity={\displaystyle \frac{TP}{P}}$$

(11) $$Precision={\displaystyle \frac{TP}{TP+FP}}$$

(12) $$F1-score=2\times {\displaystyle \frac{Precision\times Sensitivity}{Precision+Sensitivity}}$$where sensitivity, also known as recall, refers to the ratio of correctly predicted tandem duplication count (TP) to the total actual tandem duplication count (P); FP is the count of tandem duplications incorrectly predicted by the method; precision refers to the ratio of correctly predicted tandem duplication count (TP) to the total predicted tandem duplication count (TP + FP); and F1-score represents the harmonic mean of sensitivity and precision, two mutually balancing metrics. A higher F1-score indicates better performance of the detection method.

Table 3 shows the F1-score of the four methods in different configurations. Among all coverage depth and tumor purity, DTDHM consistently maintained the highest F1-score. In 450 simulated data samples, the average F1-score of DTDHM, SVIM, TARDIS, and TIDDIT were 80.0%, 56.2%, 43.4%, and 67.1%, respectively. The F1-score of DTDHM was 1.2 times higher than that of the second-best method.

Figure 6 depicts the F1-score comparison of the four algorithms in different configurations. The F1-score of these four methods improved significantly as coverage depth and tumor purity increased. The DTDHM method exhibited a relatively small variation in F1-score, indicating stable detection performance. At the coverage depth of 4X, with the change of tumor purity, the F1-score of DTDHM changed from 85.9% to 69.5%, increasing by 16.4%, while the F1-score of SVIM increased by 17.2%, the F1-score of TARDIS increased by 43.5%, and the F1-score of TIDDIT increased by 31.5%. At 6X and 8X coverage depths, the F1-scores of DTDHM improved by 14.7% and 13.8%, respectively. TARDIS was the most sensitive to changes in tumor purity, followed by TIDDIT. All methods were most effective when the tumor purity was 0.6. The results of the four methods were more stable at different coverage depths than they were with changes to tumor purity.

Figure 7 compares the sensitivity, precision, and F1-score between DTDHM and the three methods on simulated datasets. Under each set of conditions, DTDHM achieved the best F1-score, followed by TIDDIT, SVIM, and TARDIS. The average sensitivity of DTDHM was the highest at 79.9%, which was 21.9%, 47.7%, and 13.5% higher than SVIM, TARDIS, and TIDDIT, respectively. The average precision of DTDHM was 80.7%, and the average precision of SVIM, TARDIS, and TIDDIT were 57%, 82.7%, and 73.2%, respectively. SVIM achieved the lowest precision under each set of conditions, meaning it detected many TDs, most of which were false positives. TARDIS achieved the lowest sensitivity but had a higher precision, indicating this method was more conservative in its detection than the other methods. The performance of TIDDIT was better than that of both SVIM and TARDIS, ranking second in terms of F1-score. Of all the tested methods, DTDHM achieved the best balance of precision and sensitivity.

For the detection accuracy of TD boundaries, the boundary bias of 50 sets of experiments under each configuration were plotted to better show the distribution of data. The boundary bias was defined as the average base pair (bp) of the recalled TDs from the actual TD boundary deviation. The smaller the boundary bias, the more accurate the TD detection method. The comparison of boundary bias of the four methods is shown in Fig. 8. As seen in Fig. 8, the boundary bias of the four methods decreased as tumor purity and coverage depth increased. This is due to the fact that there are fewer effective SR and PEM signals in low tumor purity samples, leading to inaccurate breakpoint identification. Overall, the SVIM method exhibited smaller boundary biases at 4X coverage depth with boundary biases consistently around 36 bp, followed by the DTDHM method, with most boundary biases around 50 bp. For 6X and 8X coverage depth, DTDHM demonstrated smaller boundary biases with an average bias of 34 bp, followed by SVIM. The TARDIS method had the worst and most volatile results. The detection of these methods was also analyzed in the 10X_0.6 configuration, and the results are presented in the Supplemental File. The results of these comparisons of simulated data indicate that DTDHM is an effective and reliable TD detection method.

Application to the real samples

Because tests on simulated data cannot reflect the real detection effect of the algorithm on biological data, DTDHM was also used to analyze five real sequencing samples (NA19238, NA19239, NA19240, HG00266, and NA12891) from the 1,000 Genomes Project. The first three samples were from the YRI trio of the Yoruba family, the fourth sample was a woman of European descent from Finland, and the fifth sample was from the CEU trio of the European ancestry. TDs were predicted from a single tumor sample using whole genome sequencing (WGS) data.

DTDHM was compared with three other methods (SVIM, TARDIS, and TIDDIT). Due to the lack of ground truth files for the first four samples, it was not possible to calculate sensitivity, precision, and F1-score. Therefore, the detection effects of the four methods were analyzed by calculating the ODS of each method (Liu et al., 2020; Yuan et al., 2020, 2021). The calculation of ODS is shown in Formula (13):

(13) $$ODS=Ave{r}_{Olap}\cdot Rati{o}_{Olap},$$where $Ave{r}_{Olap}$ represents the average number of overlapping events between an algorithm and other algorithms, or the total overlapping TDs divided by the number of compared algorithms, which is roughly the number of overlapping TDs detected by the algorithm; and $Rati{o}_{Olap}$ represents the ratio of the average number of overlapping events of the algorithm to the number of predicted events. Overlapping events refer to pairwise intersections of prediction results between any algorithms. In this article, overlapping events between different algorithms were assumed to be true positives, so $Rati{o}_{Olap}$ was considered an indicator of precision.

Figure 9 shows the number of overlapping TDs predicted by the algorithm on each chromosome, in kbp. In the chord diagram, the upper half of the circle is divided into four sections denoting DTDHM, SVIM, TARDIS, and TIDDIT. The lower half of the circle is divided into 22 sections, representing autosomes 1 through 22. As shown in Fig. 9, TIDDIT detected the highest total number of overlapping TDs, followed by DTDHM, TARDIS, and SVIM. Tandem duplication structures are usually different for different chromosomes. DTDHM detected the most overlapping TDs on chromosome 1, followed by chromosome 10. On chromosome 10, the total number of overlapping events between DTDHM and other algorithms was 141.7kbp, while the total number of overlapping events of SVIM, TARDIS, and TIDDIT compared to the other algorithms were 84.36, 104.73, and 89.84 kbp, respectively. These results indicate that the four algorithms achieved good prediction results on chromosome 10.

Figure 10 records the average number of overlapping events ( $Ave{r}_{Olap}$) per chromosome in four samples for the four methods. Table 4 shows the ODS calculation results for these four methods. As shown in Fig. 10, TIDDIT detected the most number of overlapping TDs in all four samples due to the fact that TIDDIT predicts the largest number of total events. However, the precision of TIDDIT was not high. In sample NA19240, TIDDIT detected 15,059 bp overlapping TDs, but its ODS was only 945.38, which is much lower than DTDHM’s ODS of 2,503.12. This indicates that detecting a large number of TDs may cause too many false positives. TARDIS detected fewer TDs, but its ODS was lower than that of DTDHM. DTDHM detected a modest number of TDs in these samples and had a relatively higher overlapping density than the other methods.

Table 4 shows that DTDHM achieved the highest ODS in three of the real samples, followed by TARDIS, TIDDIT, and SVIM. In the four real samples, the ODS of DTDHM was 1.02, 0.94, 1.13, and 1.44 times the ODS of TARDIS, respectively. TARDIS had a relatively high ODS in each real sample, but its low $Ave{r}_{Olap}$ resulted in poorer detection performance compared to DTDHM. DTDHM maintained a high $Ave{r}_{Olap}$ and $Rati{o}_{Olap}$ and had the best ODS among the real samples.

The long-read sequencing data results of chromosome 21 for NA12891 was downloaded from the NCBI database to provide a more intuitive display of method performance (https://www.ncbi.nlm.nih.gov/). Each method was repeated 20 times for each sample, and the average performance metrics of precision, sensitivity, and F1-score were recorded for each method. As shown in Fig. 11, DTDHM demonstrated the best performance on NA12891, achieving the highest F1-score. DTDHM was able to balance precision and sensitivity to provide relatively accurate results. On NA12891, DTDHM’s F1-score was 42.1%, which was 17.1%, 20.7%, and 22.1% higher than SVIM, TARDIS, and TIDDIT, respectively. The results for TARDIS and TIDDIT were very similar. Furthermore, the precision and sensitivity of DTDHM were superior to other methods. It is worth noting that, in contrast to the results from the simulated datasets (Fig. 7), the results of many methods were not satisfactory on chromosome 21 for NA12891. This is because the types of variations in the real genome are more complex, leading to deviations in method performance between the simulated dataset and the real dataset. These results indicate that the application of the DTDHM method in actual data is relatively reliable.