VACmap: an accurate long-read aligner for unraveling complex genomic rearrangements

Domain: physics
Published: 2026-01-05
Editorial outcome: minor_revision

Abstract
Sequence alignment is essential for genomic research and clinical diagnostics, yet detecting complex rearrangements such as inversions, duplications, and gene conversions remains challenging due to allele complexity and limitations of current methods. We introduce VACmap, a non-linear mapping approach to enhance the detection and representation of all genetic variations. VACmap improves duplication detection from 20% to 90% in the Challenging Medically-Relevant Genes (CMRG) benchmark and improves characterization of complex inversions in repetitive regions and gene conversion events. It improves resolving clinically significant loci, including the LPA gene (with repetitive KIV-2 units linked to coronary heart disease), GBA1 and STRC genes (risk factors for Parkinson’s disease and hearing loss, respectively, affected by pseudogene recombination with GBAP1 and STRCP1). Here, we show that VACmap delivers better alignment accuracy and SV detection, providing a robust tool for genomic analysis and clinical insights, with potential to advance understanding of genetic diversity and disease mechanisms.

Introduction
Sequence alignment is a fundamental starting point for genomic research and clinical diagnostics, serving as the critical bridge between raw sequencing data and meaningful biological insights. Accurate alignment directly influences the quality of downstream analyses, such as variant detection, genome assembly, comparative genomics, and personalized medicine. In particular, the detection of structural variations (SVs)—genomic rearrangements such as inversions, duplications, translocations, and complex clustered events—relies heavily on precise alignment1. However, the inability of existing linear aligners to adequately represent complex SVs remains a barrier to progress in genomic science. Misaligned or misrepresented SVs hinder downstream analyses, leading to gaps in our understanding of genetic variation and its impact on human health and disease.

SVs, defined as genomic alterations of 50 base pairs or larger, are among the most impactful sources of genetic variation1–4. They affect more nucleotides in the genome than smaller variants, such as single-nucleotide polymorphisms (SNPs) or small insertions and deletions (indels)3,5–7. Consequently, their influence is being recognized across evolutionary processes, human health, and contribute to both Mendelian and complex diseases as well as cancer development8–11. Despite their critical importance, our understanding of SVs remains limited12. This is in part due to the complexities of complex SV, but also the inability of the current state of art aligners to represent them. Research has predominantly focused on simpler SV classes, such as deletions and insertions, while more intricate variations—like duplications, inversions, and complex clustered events—are often misrepresented and thus missed. This gap in our analytical approach hinders a comprehensive understanding of these important genomic phenomena13,14, despite their importance being highlighted in several studies already7–11,14–16.

These studies have mainly been driven by long-read sequencing that indeed enable the characterization of tandem repeats and thus regions where SV are predominantly observed17. Despite their advancements, long-reads are challenging to align due to their length and generally higher sequencing errors. In the past, others have highlighted that specialized aligners are required to accurately align them using predominantly linear alignments to the reference, where multiple linear subalignments represent a potential SV. These subalignments are identified through a process commonly known as the seed-chain-extend algorithm, tailored for long-read mapping. In this approach, seeds (exact matches, such as k-mers) are identified, and a co-linear subset of these seeds is selected to form chains, which are then extended bidirectionally until significant differences between the read and reference sequences are encountered. For reads without SVs, a single chain (alignment) can represent the entire read. However, the presence of SVs requires multiple chains (subalignments) to represent the read fully—for example, a read spanning an inversion typically requires three subalignments to capture the inversion and its flanking regions. Due to the complexity of SVs and their tendency to occur in repetitive regions, the seed-chain-extend approach often generates a pool of redundant subalignments, necessitating an additional step to determine the optimal subset of subalignments. For instance, minimap2 employs a co-linear chaining algorithm to identify the set of all possible co-linear chains. It then uses a greedy strategy during primary chain selection18, which determines the optimal subset of subalignments to represent the read. This process initializes an empty set Q and iteratively processes subalignments from the highest to the lowest chaining scores: if a subalignment overlaps with a chain in Q by 50% or more of the shorter subalignment’s length, it is marked as secondary; otherwise, it is added to Q, ultimately representing the read with the subalignments in Q18. In contrast, YAHA adopts a graph-based approach, leveraging Optimal Query Coverage algorithm to finds the optimal set of subalignments that cover the length of the query19. NGMLR builds on a similar strategy, enhancing it with a refined scoring function to identify the optimal combination of subalignments with the highest joint score20. Despite these advancements, post-alignment processes for determining the ideal set of sub-alignments are often inadequate because of the complexity of the allele and the underlying repetitive regions. They fail especially for duplication, inversion, and translocation, which are often missed or falsely identified by the linear alignment algorithms. Duplications, for example, are often misaligned as insertion because linear alignment algorithms prefer a single continuous alignment and treat duplication, which would require splitting reads, as insertion. Similarly, inversion and translocation, are often misaligned as splitting the reads is penalized and thus avoided.

To overcome the challenges posed by existing alignment methods, we introduce VACmap, a long-read mapping tool developed to improve the representation of all types of SVs. VACmap uses a non-linear alignment algorithm that captures an entire read as a unified, non-linear alignment. This approach streamlines the traditional alignment process by eliminating the need for splitting reads and selecting from multiple linear alignments. We demonstrate that this approach improves the representation of complex alleles, providing a more accurate and comprehensive view of SVs.

Results
The workflow of the non-linear alignment algorithm
Figure 1 gives an overview of VACmap’s non-linear mapping approach. The key important differentiation between VACmap and other approaches is implemented after initial matches between reference and read sequences have been identified. Here, existing methods try to conserve the order of all subalignments by heavily penalizing splits when searching chains of matches maintained. The linear alignment approach can efficiently model genomic alterations such as deletions, insertions, and substitutions since these don’t break the co-linearity of a chain. However, the linear approach penalizes the detection of complex SV such as duplication, inversion, translocation, or combinations of SV. In VACmap, we propose a hybrid alignment algorithm, which combines both linear and non-linear linkage approaches in a chain. In detail, VACmap represents matches as quadruples called ‘anchors’, which include the start positions in the long read and reference sequences, the strand match, and the anchor’s length. They are ordered by their end positions in the long read. The VACmap’s non-linear chaining algorithm then promotes the extension of the chain to subsequent anchors that preserve a strictly linear relationship with the preceding anchor, enhancing this connection with a positive score. Conversely, it penalizes extension to anchors that disrupt this linearity by assigning negative scores to such connections. Then the optimal non-linear alignment of the entire sequence is the chain with the highest aggregate score (the longest path). Each of these linear subalignments can be extracted by dividing the non-linear alignment at the non-linear junction, eliminating the traditional necessity of additional post-alignment steps that reconstruct genomics rearrangement from a pool of error-prone independent subalignments (See “Methods” for details).

Assessing the impact of VACmap on detecting complex SVs in synthetic data
To assess the impact of VACmap on variant detection in downstream applications, we conducted a series of tests using synthetic datasets. We generated synthetic long-read datasets containing a wide range of SVs, both simple and complex, using a custom tool we developed called VACsim, addressing the absence of simulation tools for complex SVs. VACsim introduced 30,000 SVs, each composed of 1 to 20 basic SV events, including deletions, insertions, duplications, inversions, and translocations.

As illustrated in Fig. 2a, VACmap’s alignment data enhanced SVIM’s20 ability to detect complex SVs in simulated data from PacBio CLR, PacBio HiFi, and ONT, with F1 score improvements ranging between 29.5 and 73.2 percent (refer to Supplementary Table 1). For complex SVs located within repetitive sequences, the use of VACmap-produced alignments provided gains in precision and recall, performing better than other methods by about 35.2 to 64.6 percent in F1 score (see Supplementary Table 1). Figure 2b shows the recall rates of SVs detection under different SV complexity and sequencing technology. NGMLR, minimap2, Winnowmap2, and LRA21 were shown to be adequate only for identifying complex SVs comprising up to two simple SV events. Beyond this complexity, the recall rates of SVIM decreased when using alignments from these tools. Conversely, SVIM with VACmap alignments consistently displayed sensitive and reliable SV detection across the full spectrum of SV complexities.

For precise gene copy number quantification, particularly in tandem duplications that might influence protein levels, accurate mapping is essential. To investigate the performance of copy number estimations with alignments produced by different aligners, we generated 10,000 tandem duplications on chromosome 1 using VACsim. These duplications had repeat unit sizes ranging from 100 to 500 base pairs and repeat counts between 1 and 20. SVIM was employed to estimate the copy number for each tandem duplication from the various aligners’ alignments. According to the results depicted in Fig. 2c and Supplementary Table 2, alignments from current mapping methods led to a bias in copy number estimation. There was a decline in the linear correlation between the actual and estimated copy numbers as the repeat count grew, especially within repetitive areas. In contrast, alignments from VACmap resulted in more precise copy number estimates across diverse copy number intervals and within repetitive regions. This underscores VACmap’s capability in accurately ascertaining the copy number of tandem duplications, indicating its effectiveness and accuracy in dealing with complex genomic structures.

Evaluation using genome in a bottle benchmark
We evaluated the SV detection performance of VACmap, NGMLR, Winnowmap2, minimap2, and LRA alignments using SVIM and cuteSV with the GIAB benchmark set4,20–25. Truvari26 was used to assess precision, recall, and F1 scores. Before evaluation, SVIM’s and cuteSV’s tandem duplication calls were relabeled as insertions to allow for comparability to the GIAB assembly-derived benchmark. As expected, all five alignment approaches demonstrated similar performance in detecting deletions and insertions in both GIAB tier 1 and CMRG regions (Fig. 3a, b). And the runtime of VACmap is faster than NGMLR and comparable with Winnowmap2 and LRA, but slower than minimap2. However, VACmap requires lower memory usage than the other aligners (Supplementary Table 3). It should be noted that NGMLR is no longer actively maintained, which may contribute to its performance limitations compared to more actively developed tools.

To evaluate SVIM’s and cuteSV’s sensitivity in detecting duplications using alignments from different tools, we isolated tandem duplication calls within the GIAB benchmark set using REPTYPE annotation. The results (Fig. 3c, d) showed that SVIM, using VACmap-produced alignments, exhibited high sensitivity for duplication detection, identifying approximately 70% to 80% more duplications compared to other alignment approaches in the GIAB tier 1 and CRMG regions, respectively. This is highly important for the interpretability of the impact of SV. Additionally, the SV distribution detected with VACmap alignments showed notable differences compared to other aligners (Fig. 3e and Supplementary Fig. 1). VACmap indicated that more than 67% of the sequence gain was due to duplications, consistent with previous findings14. In contrast, minimap2 attributed only 1% of the total sequence gain to duplications. This discrepancy in SV classification is critical for interpreting the biological impact of SVs, underscoring the importance of accurate SV detection.

VACmap’s ability to accurately map duplicated segments also enabled us to characterize a previously reported de novo variation27 (Fig. 4a–c and Supplementary Fig. 2). This variation, located within a Tandem Repeat (TR) region at chr14:23,280,711 (GRCh38), was originally labeled as a de novo insertion, as different insertion sizes were observed in the child (HG002: 537 bp) and the parents (HG003: 214 bp and HG004: 15 bp). However, with VACmap’s alignment, what was initially thought to be an insertion was revealed to be a 109-bp Variable Number Tandem Repeat (VNTR), with varying repeat counts in the child (five repeats) and the paternal parent (two repeats). TR regions are known to be variable in the number of repeats, often changing between generations due to mechanisms like replication slippage and unequal crossing over during meiosis. These processes can lead to differences in repeat counts, which explains the variation observed between the child and the father in this case.

This example highlights a limitation of conventional alignment algorithms, which often misinterpret duplications as insertions. Traditional aligners rely on maintaining the relative order of sequences when aligning them. However, duplications disrupt this order, making it difficult for linear aligners to correctly map such regions. As a result, duplications are often misaligned as insertions or entirely ignored. In contrast, VACmap’s non-linear alignment approach accurately handles these complex repeat structures, providing a more precise representation of the true genetic variation.

Enhance the characterization of complex inversions in repetitive regions
We then analyzed the inversion call sets generated by five different SV detection pipelines. The VACmap-SVIM call sets captured nearly all of the inversions (105 out of 116) identified by the combined call sets of minimap2, Winnowmap2, NGMLR and LRA, and additionally uncovered 97 inversions not detected by these approaches (Fig. 3f). When comparing inversions that overlapped with a previously reported call set28, the VACmap-SVIM pipeline identified nearly all the inversions (48 out of 49) detected by the other three pipelines, while also discovering 14 inversions that were missed by the other methods (Supplementary Fig. 3). Upon manual inspection of an inversion missed by VACmap-SVIM, we found a more complex structure—an inversion flanked by an inverted duplication and deletion. While VACmap could resolve this complex structure, SVIM failed to detect it because the intricate structure did not align with SVIM’s predefined rules for identifying inversions (Supplementary Fig. 4).

Thus, highlighting that inversions remain challenging to resolve because their locations are often surrounded by large segmental duplications. To further investigate this, we analyze the combined call set of 213 inversion regions from five aligners. Across all inversions, 32% (68/213) of them overlap with segmental duplications, and half of them (39/68) are only detectable through VACmap alignment. For instance, VACmap alignment enables accurate identification of a homozygous 16-kb inversion located in the SPIDR gene (Fig. 4d), a gene involved in DNA repair and associated with gonadal dysgenesis diseases29. On the contrary, other aligners’ alignments are less reliable, as they showed more mismatch bases (i.e., signal of wrongly mapping of reads20) and inconstant breakpoints across different read alignments. The standard deviation of inversion sizes called by SVIM is 291.4 for VACmap alignments and 2066.4 for NGMLR alignments, respectively. A higher variance will be considered an unreliable SV prediction and assign a lower quality score (The SVIM quality score for this inversion is 14 and 0 for VACmap and NGMLR alignments, respectively, and will be discarded).

Figure 4e demonstrates why minimap2 and other linear aligners fail to accurately pinpoint inversion breakpoints. Linear alignment methods, such as minimap2, rely on heuristic strategies like the Z-drop heuristic to infer breakpoints18. These methods monitor the alignment score and split the alignment when the score drops below a predefined threshold (indicated by the red dashed line in the figure). However, this approach often fails to identify the precise breakpoint because after the inversion, the sequence in the read is not significantly divergent from the reference. As shown in the figure, the alignment score continues to increase slowly rather than showing a sharp drop, leading minimap2 to incorrectly place the breakpoint upstream (marked by the red dashed line).

In contrast, VACmap’s non-linear alignment algorithm can simultaneously evaluate both forward and reverse strands (blue and orange curves, respectively) and automatically switch between them to maximize the alignment score. This allows VACmap to correctly identify the true inversion breakpoint, as it can seamlessly align both strands and captures subtle changes in the alignment score. The result is a more accurate alignment and a precise breakpoint, as reflected in the figure, where VACmap’s breakpoint (blue dashed line) aligns with the actual inversion. Supplementary Figs. 5–9 provide further examples of how VACmap performs better than traditional aligners in mapping complex inversions.

Improve identification on SIGLEC11::SIGLEC16 and RHCE::RHD gene conversion
Gene conversion is a challenging form of SV that is difficult to capture accurately using current alignment algorithms and SV detection tools. Figure 5a and Supplementary Fig. 10 illustrate an inversion initially misidentified by SVIM, which was actually a gene conversion event between the SIGLEC11 and SIGLEC16 genes on the maternal haplotype. These two genes share highly similar sequences in the regions encoding their extracellular domains, due to past gene conversion events30. The most recent conversion, which occurred approximately one million years ago, involved regions A in SIGLEC11 and A* in SIGLEC1630 (Fig. 5b). However, VACmap’s alignment revealed a gene conversion event involving different regions, B and B*, in these two genes.

Notably, the B* region in SIGLEC16 had previously been flagged by the GIAB consortium due to a cluster of heterozygous small variants23. However, GIAB’s alignment methods, which rely on minimap2, were unable to resolve this gene conversion, resulting in numerous false-positive SNP calls in both the GIAB CMRG benchmark set and the draft release of the GIAB T2T SV benchmark (Fig. 5d). This outcome is not surprising given minimap2’s limitations in handling complex rearrangements, as it struggles to split reads or assemblies appropriately to represent gene conversion events, leading to misalignments and erroneous variant calls.

Additionally, VACmap successfully resolved a homozygous gene conversion event between the RHCE and RHD genes, which had been inaccurately represented by existing aligners (Fig. 5c). This correction reduced over a hundred false-positive SNP and indel calls in the GIAB benchmark sets (Fig. 5e and Supplementary Figs. 11 and 12). This highlights VACmap’s ability to detect and accurately characterize gene conversion events that are typically missed or misclassified by conventional linear alignment methods.

Evaluation using the LPA, GBA1, and STRC genes
We next assessed the LPA gene to highlight a medically important region that is further improved using VACmap. The complexity of this region raises due to high diversity in the population which represents 5–40 copies of the KIV-2 repeat in the LPA gene10. This copy number is inversely correlated with human lipoprotein(a) levels, which are strongly linked to coronary heart disease10. However, quantifying the KIV-2 copy number accurately poses challenges due to repetitiveness and thus the low mappability of sequencing reads in the LPA gene region31. We assessed the performance of five mapping methods by aligning PacBio HiFi and ONT sequencing data from human samples (CHM13 and HG002) against the GRCh38 reference genome. IGV visualizations revealed that NGMLR, Winnowmap2, minimap2, and LRA produced alignments with more mismatches and less informative coverage information compared to VACmap (Fig. 6a). VACmap demonstrated an ability to accurately represent KIV-2 repeats, showing clear and distinct coverage boundaries (Supplementary Figs. 13–15).

To simplify KIV-2 copy number determination, we modified the GRCh38 reference by removing the second to sixth KIV-2 repeat units and including the first 1000 bp sequence of the follow-up KIV-1 unit (Fig. 6c). We then realigned the PacBio HiFi and ONT data to the modified reference. The IGV visualizations indicated that VACmap-produced alignments (Fig. 6b and Supplementary Figs. 16–18) showed the expected alignment scheme of both type A and type B KIV-2 units (Fig. 6e, f. Other mapping methods struggled to produce correct alignments despite the reduced complexity of the modified reference. Furthermore, the ONT reads facilitated the resolution of all 23 copies of the KIV-2 repeat unit in the CHM13 sample due to its longer read length compared to PacBio HiFi data (Supplementary Fig. 19).

Then, we aligned the GRCh38 assembly, CHM13 assembly, and HG002 assembly to the modified GRCh38 reference. The non-linear alignment of these three assemblies is shown in Fig. 6d. Consistent with previous findings10, we found the GRCh38 assembly consisted of six copies of KIV-2 repeat units with the pattern “AAABAA” where “A” indicates the type A KIV-2 repeat unit, and “B” indicates the type B KIV-2 repeat unit. In the CHM13 assembly, 23 KIV-2 repeat units were identified, following the pattern: “BBBBBBAABAAAAAAAAAAAAAA”. Similarly, the HG002 paternal assembly contains 24 KIV-2 repeat units with the pattern “BBBBBBAAABAAAAAAAAAAAAAA”, while the HG002 maternal assembly consists of 14 KIV-2 repeat units arranged as “AAAAAAAAAAAAAA”.

To further demonstrate the clinical utility of VACmap, we chose GBA1. This is a major risk factor for Parkinson’s disease32, a challenging gene to analyze33, which is prone to structural variants caused by recombination with a nearby highly homologous pseudogene (GBAP1). We previously detected using ONT long reads with adaptive sampling a pathogenic deletion which could not be correctly called after minimap2 or NGMLR alignment34. In contrast, VACmap allowed SVIM and cuteSV to correctly report the breakpoints (Fig. 7a), which is crucial in determining whether a deletion is pathogenic. Similarly, the STRC gene is a known deafness-associated gene causing mild-to-moderate hearing loss34 and is inherited in an autosomal recessive manner. However, it’s hard to detect due to its location in tandem duplication region and the presence of a highly homologous (>99%) pseudogene (STRCP1)35. By examining the GIAB samples using VACmap-produced alignments (Fig. 7b), we identified a heterozygous deletion in NA1924036 involving the loss of one copy of the CKMT1B-STRC-CATSPER2 gene cluster. However, all four other aforementioned aligners failed to pinpoint the deletion. In addition, Duplomap37, a specialized aligner for remapping reads in tandem duplications, cannot detect the deletion, since it uses minimap2 internally, which often hesitates to split reads.

Methods
Ethics statement
Ethics approval for the GBA1 carrier was provided by the National Research Ethics Service London—Hampstead Ethics Committee as part of the RAPSODI study (www.rapsodi.com)40. Informed consent was provided.

Variant-aware chaining algorithm
Algorithm overview. Traditional sequence alignment algorithms13,18–20,22 rely on linear edits—insertions, deletions, and substitutions—that preserve sequence order and orientation. While effective for simple genomic alterations, these methods struggle to represent complex rearrangements (e.g., duplications, inversions, and translocations) and are often indistinguishable from block substitutions. In cases of clustered rearrangements, linear methods attempt to reconstruct the structure by selecting from a pool of local subalignments, but this approach rarely yields a globally optimal alignment, especially in repetitive or highly rearranged regions41.

VACmap introduces a hybrid approach to address these limitations. Unlike linear alignment methods, VACmap employs the Variant-aware Chaining (VAC) algorithm, which integrates linear and non-linear edits within a weighted directed acyclic graph (DAG) framework. The algorithm identifies exact k-mer matches (anchors) between a reference and a long-read sequence, constructs a DAG where anchors are nodes, and connects them with edges representing possible alignments. Edges are classified as normal (for co-linear alignments) or variation (for rearrangements), with weights adjusted by penalties to account for indels and rearrangement complexity. By finding the longest path in the DAG, VACmap produces a unified, globally optimal alignment that accurately captures genomic rearrangements without relying on post-alignment selection.

Anchor identification. VACmap begins by identifying exact k-mer matches between the reference and long-read sequence (x, y, s, k), where:
x: Start position in the long-read sequence.
y: Start position in the reference sequence.
s: Orientation (1 for forward strand, −1 for reverse strand).
k: Length of the matched k-mer.

Graph construction. Anchors are modeled as nodes in a weighted DAG. A directed edge from node i on node j is created if the nodes satisfy:
x i+k i ≤x j+k j (1)

This condition tests for all node j that satisfy it and ensures that the alignment respects the read sequence order, allowing overlaps between nodes. The overlap size is defined as:
overlapsize=max(x i+k i −x j,0) (2)

This accounts for cases where anchors partially overlap in the long-read sequence.

Edge classification and weighting. Edges are classified as normal or variation based on the relative positions and orientations of the connected nodes. The initial edge weight, or bonus, is calculated as:
bonus i,j=min(x j+k j −x i −k i,k j ) (3)

This bonus reflects the length of the aligned region between nodes. To determine the edge type, we compute the sequence gain or loss (diff) as:
diff= |readgap−refgap| (4)

where readgap and refgap are the distances between node j’s first not overlapped base pair to node i and node i’s last base pair in the query and reference, respectively:
readgap=max(x j−x i−k i,0) (5)

refgap= |y j+y i−y j−y i| (6)

An edge is classified as normal if it satisfies:
diff<maxdiff,0≤refgap<maxgap,s i=s j (7)

Here, maxdiff (default: 50 for first round chaining, 30 for second) regulates small indels, maxgap (default: 1000 for first round chaining, 100 for second) prevents misalignments. Edges not meeting these criteria are classified as variation edges, representing rearrangements.

Penalty calculations. For normal edges, an additional penalty (NP) account for small indels:
a= 0:5 if diff<10
2 if diff≥10 (8)

NPi,j=0:01×k j*diff+α×log 2(max(diff,1))+0:1×log 2(max(readgap,1)) (9)

For a large read gap (default: >30), we use a larger penalty:
NPi,j=0:01×k j*diff+α×log 2(max(diff,1))+0:5*readgap (10)

For variation edges, the penalty (VP) accounts for rearrangements:
VPi,j=β+min(30, 1/γ * ln(max(1,diff))) (11)

Where the default settings for β is 59, 40, 30, and 30, for mode L, H, S, and R, respectively. For parameter γ, the default settings are 36. After applying the additional penalty, the final weight for normal and variation edges are defined in Eqs. (12) and (13), respectively.

w(i,j)=bonus i,j −NPi,j (12)
w(i,j)=bonus i,j −VPi,j (13)

Optimization. Here, we consider the optimal alignment of the read and reference sequence to be the longest path in the weighted directed graph. To find the longest path among N nodes can be computed in O(N2) time using dynamic programming. In detail, for a given node j, we can determine its maximum score S(j) and its best predecessor node i by utilizing Eq. (14):
S(j)=max maxi<j (S(i)+w(i,j)),k j (14)

However, we don’t need to test all node i because node j’s best predecessor score S(i) is always higher than the current highest score S(k) minus (β+γ+20). So, for node j, we define a set M j which includes all nodes that score higher than S(k) minus (β+γ+20). Then we only need to test nodes in set M j to identify the best predecessor of node j.

This further reduces the time complexity to O(hN), where h is equal to the average size of set M.

We observed that the best predecessor of node j is usually among the top-ranked nodes in set M j when sorted in descending order by their scores S(i). To exploit this property and enhance computational efficiency, we initialize the score of node j as S(j)=k j and then iterate over all nodes i in M j (where i<j) in descending order of their scores S(i) to identify the best predecessor. The iteration stops early if S(i)+k j<S(j), as this condition indicates that further predecessors cannot improve the score. In practice, we maintain an array of nodes in descending order by their scores using binary search. Then the overall time complexity of the algorithm is O(kN+NlogN). The average value of k for HiFi, ONT, and CLR data in the first round of chaining is 21, 33, and 44, respectively.

Once the maximum score of all nodes is computed, we can identify the highest scores non-colinear chain by backtracking. The optimal set of colinear subalignments is recovered by discarding variation edges in the highest score chain.

Map quality calculations. For the map quality of the highest score chain, VACmap uses minimap2’s equation to compute the map quality:
mapquality=40×(1−f2/f1)×min(1, m/10)×ln(f1) (15)

where m is the number of anchors on the highest scores chain, f1 is the chaining score, and f2 is the score of the second high chain.

Local index and anchor reduction. When sequencing errors or SNPs occur in clusters, accurate mapping becomes challenging, especially when utilizing a large k-mer size setting and a minimizer k-mer sampling strategy. To mitigate this issue, VACmap adopts a similar approach used in a previous study22 by constructing a local index with a smaller k-mer size, typically a 9-mer. The local index is built by collecting all possible k-mers from the previously computed highest and second highest score chains’ covered reference genome region. Subsequently, all the k-mers in the long read are employed to query the index, obtaining matching information known as anchors. The highest-scoring chain is then recomputed based on these anchors.

Note that the runtime of the VAC algorithm depends on the number of anchors. Using a smaller k-mer size setting often results in an increased number of anchors, consequently escalating computational time. To address this issue, two strategies are employed to reduce the number of anchors and improve runtime efficiency. Firstly, VACmap iterates through the anchors and removes an anchor if the distance between the anchor and the previously computed highest-scoring chain exceeds a certain threshold (2000bp). The distance between an anchor and the previous chain is defined as the distance between the anchor and its closest anchor in the chain under the reference genome coordinate. Determining the closest anchor in the chain can be accomplished through a binary search, with a time complexity of O(logN), where N is the number of anchors. Secondly, VACmap merges two anchors (x i, y i, s i, k i) and (x j, y j, s j, k j) into a new anchor (x i, min(y i, y j), s i, x j+k j –x i) if they satisfy certain conditions (Eq. (16)).

Two anchors can be merged only if they are overlapping or close adjacent and share identical sequences in their overlapping regions. To avoid producing large anchors—which could expand the search space during best predecessor computation—we limit the maximum output size (default: 19) for merged anchors, thereby reducing computational demand in the chaining process without compromising alignment quality.

The sorting of N anchors in ascending order based on their position in the long read allows the merging process to be computed in O(N) time using a hash table. The hash table utilizes an integer key and maintains a list of merged anchors as its value. For anchor i, if s i equals 1, then the key is set to y i −x i; otherwise, it is set to −(y i +x i ). VACmap tests the last anchor j in the list corresponding to a key. If the two anchors are overlapped or closely adjacent and the covered sequence in reference and long read is identical, VACmap updates the last anchor in the list; otherwise, anchor i is appended to the list. If a key does not exist in the hash table, a new key-value pair is inserted. Finally, the merged anchors set can be obtained by traversing the hash table.

Quality control of linear subalignment. The presence of sequencing error and closely situated SNPs complicate the process of gathering sufficient anchors for the chaining process. To accommodate the occurrence of sequencing errors and clustered SNPs, VACmap employs a sizeable parameter, maxgap, with a default setting of 100, within the VAC algorithm. However, a large maxgap value has the potential drawback of mistakenly aligning sequences that are not actually related. To address this problem, we assess the error ratio within each linear subalignment and exclude those that demonstrate a high propensity for errors. This error ratio is ascertained by employing edlib42 to calculate the edit distance between the matched sequences in the reference genome and the long read. The edit distance obtained is then normalized by dividing the length of the shorter sequence involved in the comparison. A linear subalignment is deemed unreliable and is therefore rejected if its error ratio surpasses the predetermined threshold: 0.2 for sequences generated using PacBio CLR or ONT, and 0.1 for sequences generated using PacBio HiFi.

Mode selection. VACmap provides four mapping modes (L, H, S, R) tailored for specific alignment purposes. Mode L is recommended for high-accuracy long-read data (e.g., PacBio HiFi), using a large variation penalty (β=59) and a low sequence divergence threshold (0.1) to retain only high-confidence alignments. Mode H is preferred for noisy long-read data (e.g., ONT), applying a slightly smaller variation penalty (β=40) and a higher sequence divergence threshold (0.3) to accommodate higher sequencing errors. Mode S is designed to improve sensitivity for small-scale rearrangements (β=30) and a high sequence divergence threshold (0.5) to capture subtle variations. Mode R is an experimental mode intended solely for testing purposes, employing a fixed variation penalty (VP =30) and a divergence threshold of 0.5. Due to its developmental nature, Mode R is not recommended for normal usage, and users should rely on Modes L, H, or S for standard applications.

Synthetic data simulation and evaluation
To assess the performance of different mapping methods in downstream SV detection tasks, we employed the CHM13 T2T human reference genome43 to generate two sets of synthetic datasets. For both simple and complex SVs simulation, we utilized VACsim to randomly introduce 30,000 SVs into the CHM13 T2T reference and generated the altered genome sequence containing both simple and complex SVs. VACsim determined the SV complexity (measured by the number of simple structural variation events within) of each simulated SV by randomly sampling from the range 1–20. The types of simple structural variation events were sampled from five simple SV events, namely deletion, insertion, duplication, inversion, and translocation, with respective probabilities of 0.24, 0.24, 0.24, 0.24, and 0.04. The size of each simple SV event ranged from 100 to 1000 bps. Additionally, SVIM (v1.4.2) employed a pattern match strategy to detect SVs, requiring each simple SV event to conform to the “Normal+Variation+Normal” pattern. Therefore, we spaced adjacent simple SV events by a normal sequence of 200 base pairs in size. For the copy number estimation task, VACsim was used to simulate and implant 10000 tandem duplications with repeat unit sizes ranging from 100 to 500 bps and repeat numbers ranging from 1 to 20 into chromosome 1 of the CHM13 T2T reference genome sequence. The SNPs were introduced into the altered sequences, and SURVIVOR generated simulated PacBio and ONT long-read datasets.

Five long-read mapping methods, namely VACmap, NGMLR, Winnowmap2, minimap2, and LRA, were utilized to align the simulated long reads to the T2T reference genome sequences. For parameter setting of VACmap, we use mode L to align synthetic PacBio HiFi data and mode H to align synthetic PacBio CLR and ONT data. The SVIM (v1.4.2) SV detection tool was employed to detect SVs from the alignments generated by the different mapping methods. For 40×, 20×, 10×, and 5× read coverage, the read support of SVIM called SVs is set to 10, 5, 3, and 2, respectively. Following previous studies6, we evaluated the performance of CSV calls by decomposing CSVs into individual simple SV events and evaluating each event separately.

Furthermore, following the procedure outlined in a previous study13, we employed Mashmap44 to identify repetitive regions in the reference genome sequence with sequence similarity exceeding 95% and a sequence length greater than 10,000bps. Subsequently, we analyzed the performance of downstream SV detection within these repetitive regions using alignments produced by different mapping algorithms.

Evaluation with GIAB call set
We evaluated five mapping methods, namely VACmap, NGMLR, Winnowmap2, minimap2, and LRA using the GIAB Tier1 (v0.6) and CMRG (1.00) SV benchmark set for the HG002 human sample. These benchmark sets are available at the following links: https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/NIST_STRUCTURALVARIANT_v0.6/ and https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/AshkenazimTrio/HG002_NA24385_son/CMRG_v1.00/, respectively. PacBio CLR, PacBio HiFi, and ONT sequencing data were accessible at s3://giab/data/AshkenazimTrio/HG002_NA24385_son/. For parameter setting of VACmap, we use mode L to align synthetic PacBio HiFi data and mode H to align synthetic PacBio CLR and ONT data. The SVIM (v1.4.2) and cuteSV (v2.1.2) were used to detect SVs. For PacBio CLR (69×), PacBio HiFi (30×), and ONT (50×) data, the read support of SVIM and cuteSV called SVs is set to 5. Except for the cuteSV called SVs using PacBio CLR data, the read support is set to 10 due to cuteSV needs a higher read support setting to filter out false-positive results (Supplementary Fig. 20). Subsequently, we evaluated the SVs against the GIAB benchmark set using Truvari (v2.0.0). Command line parameters provided to these tools are listed in Supplementary Table 4.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors present an algorithm for aligning long reads in order to improve calling complex structural variants, such as duplications and inversions. They have implemented their algorithm in software called VACmap, which they have also tested with both simulated and biological data, as well as compared it to other commonly used software.\nBoth improving third-generation sequencing read alignment and refining structural variant calling are active topics in bioinformatics research. Due to my background in computer science, I concentrated on assessing the presented algorithm.\nRegrettably the manuscript has several shortcomings in describing the algorithm. In particular, I find Variant-aware Chaining Algorithm section quite difficult to follow and hence comparison to existing methods is hard. To improve the section, the authors could present one equation at a time with an explanation and possibly start by providing a high-level description of the algorithm (but more detailed than those in Figure 1 and Supplementary Figure 20).\nAdditionally, when reviewing previous research, the authors point out why representing complex structural variants using a linear reference genome is problematic but do not consider graph-based approaches. These are topical in bioinformatics research. Specifically, there have been efforts to make use of cyclic [3] and bidirected genome graphs [4]. With these, structural variants such as duplications and inversions can be represented [5]. The authors should compare their approach to those at least on some level.\nSpecific issues have been outlined below.\nMain\n- 61–64: The class of algorithms the authors refer to is commonly known in previous research as seed-and-extend algorithms; stating this would make the paragraph more straightforward.\n- 70–11: It is unclear to me to what the authors refer when stating that minimap2 selects “the highest-scoring sub-alignment at each step until all alignable regions of the read are covered”. Minimap2 uses a co-linear chaining algorithm [1, § 2.1.1] to find the best match of a set of anchors (which I believe the authors call sub-alignments) in the read and in the reference. While a heuristic is applied to only consider a fixed number of matches, this is not a greedy approach. However, a greedy approach is applied when identifying primary chains [1, § 2.1.3].\n- 71–74: The authors state that NGMLR improves the approach of YAHA by utilising dynamic programming. However, YAHA also uses dynamic programming [2, § 2.4].\n- 99–100: While complex structural variants may be difficult to detect using co-linear chaining with a linear reference genome, the authors should compare their approach to graph-based aligners at least to a certain extent.\nMethods\n- 381–382: The authors explain that the acyclicity of the graph generated from the anchors follows from the topological order of its nodes, but the authors do not define an order relation for the nodes. In particular, is the order determined by e.g. x or y only, by all of the elements of an anchor from left to right, or is the subsequence length considered when comparing the co-ordinates (i.e. is (x1, y1, s1, k1) < (x2, y2, s2, k2) only if e.g. x1 + k1 < x2)?\n- 385–393: The equations have been presented as a long list with explanations in the next paragraph. The meaning of some variables (including a, readgap) has been omitted. All variables should be defined clearly.\n- 390: It is unclear to me, is the intention to define a conditionally (for which notation similar to that in e.g. Equation 8 should be used) or is it something else. Additionally I did not see a being used in the other equations.\n- 397: I believe the term hyperparameter is from machine learning and I have not seen it in use elsewhere. The authors could call these just parameters.\n- 423: I could not find a definition for the set M.\n- 448: Why are specifically these conditions considered when merging anchors? What is the significance of the value 19?\n- 473: I was not able to find a URL where the synthetic data is available.\nMinor issues\n- The authors publish their implementation under the GNU General Public Licence version 3.0. I would like to point out that all files published under the GPL should contain a licence notice, which currently is not the case. If needed, the authors can refer to the instructions on how to apply the GNU licences for further clarification: https://www.gnu.org/licenses/gpl-howto.html\n- Ref. 14 only has the URL of the paper in question\n- Ref. 29 seems incomplete in the sense that at least the journal is missing\n- 98: collinearity → co-linearity (or colinearity)\n- 384: Formulae are more commonly called equations in this context.\n- 394: collinearity → co-linearity (or colinearity)\n- 446: O(log^N) → O(log N)\n- 528–529: All the three links are the same.\nReferences\n1. Heng Li. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics, 34(18):3094–3100, May 2018.\n2. Gregory G. Faust and Ira M. Hall. Yaha: fast and flexible long-read alignment with optimal breakpoint detection. Bioinformatics, 28(19):2417–2424, July 2012.\n3. Jyotshna Rajput, Ghanshyam Chandra, and Chirag Jain. Co-linear chaining on pangenome graphs. Algorithms for Molecular Biology, 19(1):4, 2024.\n4. Mikko Rautiainen and Tobias Marschall. Graphaligner: rapid and versatile sequence-to-graph alignment. Genome Biology, 21(1):253, 2020.\n5. Benedict Paten, Adam M Novak, Jordan M Eizenga, and Erik Garrison. Genome graphs and the evolution of genome inference. Genome research, 27(5):665–676, 2017.\n(Remarks on code availability)\nThe Conda environment could not be installed on my macOS laptop. I currently do not have access to a suitable Linux workstation for running the complete experiments. Having a very simple set of test data would have been useful.", "reviewer_2": "The authors are introducing VACmap, a new long-read mapping tool that significantly enhances the representation of all types of SVs. The tool's innovative approach employs a non-linear chaining algorithm to connect subalignments, resulting in markedly improved detection of complex structural variants and more accurate quantification of copy number variations. The extensive supplementary materials, including tables and figures, provide compelling evidence that substantiates the mapping enhancements that lead to improved structural variant calling. The work addresses a highly significant area in clinical genetics, as structural variants are increasingly recognized as important diagnostic targets in medical practice. I particularly appreciate the performance in the case of the gene conversion example.\nVACmap is evaluated against 3 other long-reads aligners using both simulated data and GIAB samples. Their methodology involves mapping long reads to the T2T human reference genome followed by structural variant calling with SVIM. I suggest expanding the comparison to include more contemporary aligners such as LRA or LAST, while mentioning in the manuscript that NGMLR is no longer actively maintained. The manuscript's impact would be substantially enhanced by incorporating an additional structural variant caller like CuteSV or Sniffles2, both of which have demonstrated superior performance in recent benchmarks (like https://www.nature.com/articles/s41598-024-56604-2 Benchmarking long-read aligners and SV callers for structural variation detection in Oxford nanopore sequencing data)\nFurthermore, given that pseudogenes can originate from specific SV events, it would be valuable to examine VACmap's capacity to identify pseudogenes, particularly in comparison to duplomap, aligner dealing well with segmental duplication alignment. A compelling case study would be the STRC gene (implicated in hearing loss) and its adjacent pseudogene, which shares over 99.6% sequence identity with STRC.\nConcerning the VACmap software implementation: How do the mapping quality scores vary across the different mapping modes (H, L, S, R parameter settings)? I recommend including more specific guidance on selecting appropriate modes for various applications. Additionally, please specify in the manuscript which mapping mode was employed to generate the reported F1 scores.\nMinor textual errors:\nline 121 \"As illustrated in Figure 2a, VACmap's alignment data markedly enhanced SVIM’s[21] ability to detect complex SVs in simulated data from PacBio CLR, PacBio HiFi, and ONT, with F1 score improvements ranging between 31.9 to 68.9 percent,\" but on the Fig 2a there is no SVIM shown not mentioned in the figure caption.\nline 210 typo \"arcuate\" accurate\nline 212 \"a gene is involved\", remove is\nline 215, 216 - rephrase sentence \"It’s evidence in the size standard deviation of the inversion in SVIM calls is 254.2 and 2066.4 for VACmap and NGMLR alignment, respectively”\nline 218 add and\" between VACmap and NGMLR? (the SVIM quality score for this inversion is 13, 0 for VACmap, NGMLR alignment) and add \"will\" be discarded\nfigure 4 caption 4f is 4e (reference in text is correct)\nI therefore support publication of this manuscript pending the above minor revisions."}, "author_response": "Dear Editor and Reviewers,\nWe sincerely thank you for your time and effort in reviewing our manuscript, titled \" VACmap: An\nAccurate Long-Read Aligner for Unraveling Complex Genomic Rearrangements\". We are deeply\ngrateful for the insightful and constructive comments provided by the reviewers, which have\nsignificantly improved the clarity, rigor, and impact of our work. In this revision, we have carefully\naddressed each reviewer’s suggestions and incorporated their feedback to strengthen the manuscript.\nBelow, we summarize the major changes made in the revised version of our manuscript:\n1. Enhanced Description of the Variant-aware Chaining Algorithm: Addressing Reviewer\n#1’s concern about the clarity of the Variant-aware Chaining (VAC) algorithm section, we\nhave significantly revised this section to improve readability. We added a high-level\noverview, presented equations individually with clear explanations, defined all variables,\nand restructured the section into clear subsections (Overview, Anchor Identification, Graph\nConstruction, Edge Classification, Penalty Calculations, Optimization). This revision\nfacilitates understanding and comparison with existing linear alignment methods.\n2. Incorporation of Graph-Based Method Comparison: Following Reviewer #1’s\nsuggestion, we included a detailed comparison of VACmap with graph-based aligner\nminigraph in the Discussion section and added Supplementary Note 1: Comparison of\nVACmap with Graph-Based Method. This note includes case studies demonstrating\nVACmap’s superior performance in resolving duplications, inversions, and precise\nbreakpoints, addressing limitations of co-linear matching in graph-based approaches.\n3. Expanded Comparative Analysis with Additional Aligners and SV Callers: In response\nto Reviewer #2’s suggestion, we expanded the evaluation in the Results section to include\nthe contemporary aligner LRA alongside minimap2, NGMLR, and Winnowmap2, using\nGIAB samples. We also incorporated SV calling with cuteSV in addition to SVIM,\nreferencing a recent benchmark study (Helal et al., 2024, Sci Rep). We noted that NGMLR\nis no longer actively maintained, enhancing the manuscript’s relevance.\n4. Case Study on STRC and Pseudogene STRCP1: Per Reviewer #2’s recommendation, we\nadded a case study on the STRC gene and its pseudogene STRCP1 in the Results section.\nUsing GIAB sample NA19240, we demonstrated VACmap’s ability to detect a\nheterozygous deletion in the CKMTIB-STRC-CATSPER2 gene cluster, which was missed\nby NGMLR, LRA, Winnowmap2, minimap2, and Duplomap (Figure 7b). This highlights\nVACmap’s strength in resolving pseudogene-associated SVs.\n5. Improved Guidance on Mapping Modes and Software Accessibility: Addressing\nReviewer #2’s query on mapping quality scores and mode selection, we added a “Mode\nSelection” subsection in the Methods section, providing detailed guidance on choosing\nmodes (L, H, S, R) for different data types (e.g., PacBio HiFi, ONT). We clarified that\nMode L was used for PacBio HiFi and Mode H for PacBio CLR and ONT data to generate\nF1 scores. Additionally, we added a test dataset to the VACmap GitHub repository\n(https://github.com/micahvista/VACmap) to verify installation on Linux systems,\naddressing Reviewer #1’s concern about code accessibility.\n6. Correction of Minor Textual Errors and References: We addressed all minor textual\nerrors, including updating figure captions (e.g., Figure 2a, Figure 4), correcting typos, and\nimproving sentence clarity. We also updated incomplete references and corrected data links.\nIn the following, we present one-by-one responses to the reviewers’ comments.\nReviewer #1 (Remarks to the Author):\nThe authors present an algorithm for aligning long reads in order to improve calling complex\nstructural variants, such as duplications and inversions. They have implemented their\nalgorithm in software called VACmap, which they have also tested with both simulated and\nbiological data, as well as compared it to other commonly used software.\nBoth improving third-generation sequencing read alignment and refining structural variant\ncalling are active topics in bioinformatics research. Due to my background in computer science,\nI concentrated on assessing the presented algorithm.\nResponse: Thank you for recognizing the significance of our research topic and for your time and\neffort in reviewing our manuscript. We sincerely appreciate your thoughtful and constructive\ncomments. In the revised manuscript, we have carefully considered and addressed each of your\ncomments and suggestions. Please find our detailed, point-by-point responses below.\nComments:\n-1）Regrettably the manuscript has several shortcomings in describing the algorithm. In\nparticular, I find Variant-aware Chaining Algorithm section quite difficult to follow and hence\ncomparison to existing methods is hard. To improve the section, the authors could present one\nequation at a time with an explanation and possibly start by providing a high-level description\nof the algorithm (but more detailed than those in Figure 1 and Supplementary Figure 20).\nResponse:\nWe acknowledge the reviewer’s concern that the Variant-aware Chaining (VAC) algorithm section\nwas challenging to follow, particularly due to the dense presentation of equations and technical\ndetails. To address this, we have substantially revised this section to improve readability, enhance\nclarity, and better facilitate comparison with existing methods. The key revisions include:\n(1) High-level Overview: We have revised the introductory subsection to provide a clearer and\nmore accessible overview of the VACmap algorithm. This overview outlines the algorithm’s\ngoals, its distinction from conventional linear alignment methods, and its graph-based approach\nto handling complex genomic rearrangements. This overview serves as a roadmap for the\ndetailed explanation that follows.\n(2) Stepwise Equation Presentation: To enhance clarity, each equation is now presented\nindividually, accompanied by a clear explanation of its purpose and role within the algorithm.\n(3) Enhanced Comparison to Existing Methods: We have expanded the discussion of how\nVACmap differs from traditional linear alignment algorithms. Specifically, we highlight the\nlimitations of linear methods in addressing complex SVs such as duplications, inversions, and\ntranslocations. We then explain how VACmap’s integration of linear and non-linear edits\novercomes these challenges to achieve globally optimal alignments.\n(4) Improved Structure and Flow: The section has been reorganized into clear defined subsections:\nalgorithm overview, anchor identification, graph construction, edge classification, penalty\ncalculations, and optimization. This makes it easier for readers to understand its novelty and\nadvantages over existing methods. The updated section is provided below for reference.\nVariant-aware Chaining Algorithm\nAlgorithm Overview\nTraditional sequence alignment algorithms [13, 18-20, 22] rely on linear edits—\ninsertions, deletions, and substitutions—that preserve sequence order and orientation.\nWhile effective for point mutations and small indels, these methods struggle with complex\ngenomic rearrangements such as duplications, inversions, and translocations. For example,\nduplications may be misidentified as insertions, and inversions are often indistinguishable\nfrom block substitutions. In cases of clustered rearrangements, linear methods attempt to\nreconstruct the structure by selecting from a pool of local subalignments, but this approach\nrarely yields a globally optimal alignment, especially in repetitive or highly rearranged\nregions [41].\nVACmap introduces a novel approach to address these limitations. Unlike linear\nalignment methods, VACmap employs the Variant-aware Chaining (VAC) algorithm,\nwhich integrates linear and non-linear edits within a weighted directed acyclic graph\n(DAG) framework. The algorithm identifies exact k-mer matches (anchors) between a\nreference and a long-read sequence, constructs a DAG where anchors are nodes, and\nconnects them with edges representing possible alignments. Edges are classified as\nnormal (for collinear alignments) or variation (for rearrangements), with weights adjusted\nby penalties to account for indels and rearrangement complexity. By finding the longest\npath in the DAG, VACmap produces a unified, globally optimal alignment that accurately\ncaptures genomic rearrangements without relying on post-alignment selection.\nAnchor Identification\nVACmap begins by identifying exact k-mer matches between the reference and long-\nread sequences. These matches, called anchors, are represented as quadruples (x, y, s, k),\nwhere:\nx: Start position in the long-read sequence.\ny: Start position in the reference sequence.\ns: Orientation (1 for forward strand, -1 for reverse strand).\nk: Length of the matched k-mer.\nGraph Construction\nAnchors are modeled as nodes in a weighted DAG. A directed edge from node[i] to\nnode[j] is created if the nodes satisfy:\n𝑥 +𝑘 ≤ 𝑥 +𝑘 (1)\n𝑖 𝑖 𝑗 𝑗\nThis condition ensures that the alignment respects the read sequence order, allowing\noverlaps between nodes. The overlap size is defined as:\no𝑣𝑒𝑟𝑙𝑎𝑝 𝑠𝑖𝑧𝑒 = 𝑚𝑎𝑥(𝑥 +𝑘 −𝑥 ,0) (2)\n𝑖 𝑖 𝑗\nThis accounts for cases where anchors partially overlap in the long-read sequence.\nEdge Classification and Weighting\nEdges are classified as normal or variation based on the relative positions and\norientations of the connected nodes. The initial edge weight, or bonus, is calculated as:\nb𝑜𝑛𝑢𝑠 = 𝑚𝑖𝑛(𝑥 +𝑘 −𝑥 −𝑘 ,𝑘 ) (3)\n𝑖,𝑗 𝑗 𝑗 𝑖 𝑖 𝑗\nThis bonus reflects the length of the aligned region between nodes. To determine\nedge type, we compute the sequence gain or loss (diff) as:\nd𝑖𝑓𝑓 = |𝑟𝑒𝑎𝑑𝑔𝑎𝑝−𝑟𝑒𝑓𝑔𝑎𝑝| (4)\nwhere readgap and refgap are the distances between node[i]’s first not overlapped\nbase pair to node[i] and node[i]’s last base pair in the query and reference, respectively:\nr𝑒𝑎𝑑𝑔𝑎𝑝 = 𝑚𝑎𝑥(𝑥 −𝑥 −𝑘 ,0) (5)\n𝑗 𝑖 𝑖\ny +𝑜𝑣𝑒𝑟𝑙𝑎𝑝 𝑠𝑖𝑧𝑒−𝑦 −𝑘 ,if 𝑠 = 𝑠 = 1\n𝑗 𝑖 𝑖 j i\ny −𝑦 −𝑏𝑜𝑛𝑢𝑠,if 𝑠 = 𝑠 = −1\n𝑖 𝑗 j i\nr𝑒𝑓𝑔𝑎𝑝 = (6)\ny +𝑜𝑣𝑒𝑟𝑙𝑎𝑝 𝑠𝑖𝑧𝑒−𝑦 +1,if 𝑠 = 1,𝑠 = −1\n𝑗 𝑖 j i\ny +𝑏𝑜𝑛𝑢𝑠−1−𝑦 −𝑘 ,if 𝑠 ≠ 𝑠,𝑠 = −1\n{ 𝑗 𝑖 𝑖 j i i\nAn edge is classified as normal if it satisfies:\nd𝑖𝑓𝑓 < 𝑚𝑎𝑥𝑑𝑖𝑓𝑓,0 ≤ 𝑟𝑒𝑓𝑔𝑎𝑝 < 𝑚𝑎𝑥𝑔𝑎𝑝,𝑠 = 𝑠 (7)\n𝑖 𝑗\nHere, maxdiff (default: 50 for first round chaining, 30 for second) regulates small\nindels, maxgap (default: 1000 for first round chaining, 100 for second) prevents\nmisalignments. Edges not meeting these criteria are classified as variation edges,\nrepresenting rearrangements.\nPenalty Calculations\nFor normal edges, an additional penalty (NP) account for small indels:\n0.5 if 𝑑𝑖𝑓𝑓 < 10\na = { (8)\n2 if 𝑑𝑖𝑓𝑓 ≥ 10\n𝑁𝑃 = 0.01×𝑘 ∗𝑑𝑖𝑓𝑓+α×𝑙𝑜𝑔 (max(𝑑𝑖𝑓𝑓, 1))+0.1×𝑙𝑜𝑔 (max(𝑟𝑒𝑎𝑑𝑔𝑎𝑝, 1)) (9)\n𝑖,𝑗 𝑗 2 2\nFor large readgap (default: >30), we use a larger penalty:\n𝑁𝑃 = 0.01×𝑘 ∗𝑑𝑖𝑓𝑓+α×𝑙𝑜𝑔 (max(𝑑𝑖𝑓𝑓, 1))+0.5∗readgap (10)\n𝑖,𝑗 𝑗 2\nFor variation edges, the penalty (VP) accounts for rearrangements:\n30+0.5×ln(max (1,𝑑𝑖𝑓𝑓))\n𝑉𝑃\n𝑖,𝑗\n= 𝛽 +𝑚𝑖𝑛{min(10, 𝑑𝑖𝑓𝑓 )+min(30, 𝑑𝑖𝑓𝑓 ) (11)\n100 1000\n𝛾\nWhere the default settings for β is 59, 40, 30 and 30, for mode L, H, S and R,\nrespectively. For parameter γ the default settings are 36. After applying the additional\npenalty, the final weight for normal and variation edge are defined in equation 12 and 13,\nrespectively.\nw(𝑖,𝑗) = 𝑏𝑜𝑛𝑢𝑠 −𝑁𝑃 (12)\n𝑖,𝑗 𝑖,𝑗\nw(𝑖,𝑗) = 𝑏𝑜𝑛𝑢𝑠 −𝑉𝑃 (13)\n𝑖,𝑗 𝑖,𝑗\nOptimization\nHere we consider the optimal alignment of the read and reference sequence is the\nlongest path in the weighted directed graph. To find the longest path among N nodes can\nbe computed in O(N2) time using dynamic programming. In detail, for a given node[j],\nwe can determine its maximum score S(j) and its best predecessor node[i] by utilizing the\nequation 14:\n𝑆(𝑗) = 𝑚𝑎𝑥(𝑚𝑎𝑥 (𝑆(𝑖)+𝑤(𝑖,𝑗)),𝑘 ) (14)\n𝑖<j 𝑗\nHowever, we don’t need to test all node[i] because node[j] best predecessor score\n𝑆(𝑖) always higher to the current highest score S[k] minus (𝛽+𝛾+20). So for node[j]\nwe define a set M which include all nodes who score higher than S[k] minus (𝛽+𝛾+\n20). Then we only need to test nodes in set M to identify the best predecessor of node[j].\nj\nThis further reduce the time complexity to O(hN), where h equal to the average size of set\nM.\nHowever, we observed that the best predecessor of node[j] tends to be among the\ntop 20% of nodes in set M when sorted in descending order by their scores S[i]. To exploit\nj\nthis property and enhance computational efficiency, we initialize the score of node[j] as\nS[j] = k and then iterate over all nodes i in M (where i < j) in descending order of their\nj j\nscores S[i] to identify the best predecessor. The iteration stops early if S[i] + k < S[j],\nas\nj\nthis condition indicates that further predecessors cannot improve the score.\nOnce the maximum score of all nodes is computed, we can identify the highest\nscores non-colinear chain by backtracking. The optimal set of colinear subalignments is\nrecovered by discarding variation edges in the highest score chain.\n-2）Additionally, when reviewing previous research, the authors point out why representing\ncomplex structural variants using a linear reference genome is problematic but do not\nconsider graph-based approaches. These are topical in bioinformatics research. Specifically,\nthere have been efforts to make use of cyclic [3] and bidirected genome graphs [4]. With these,\nstructural variants such as duplications and inversions can be represented [5]. The authors\nshould compare their approach to those at least on some level.\nResponse:\nWe sincerely appreciate the reviewer’s insightful comment regarding the relevance of graph-based\napproaches in bioinformatics for representing structural variants (SVs). In response, we have\nincorporated a detailed comparison between VACmap and a representative graph-based method in\nthe supplementary materials and added a brief discussion in the main manuscript to contextualize\nthese findings.\nWe acknowledge the reviewer’s point that graph-based approaches, such as cyclic [1] and bidirected\ngenome graphs [2], have been developed to better represent SVs like duplications and inversions\n[3]. These methods provide a flexible framework for encoding SVs as graph structures, potentially\nenhancing the representation of genetic diversity across populations. However, similar to linear\nreference genome methods, graph-based tools such as minigraph [4] and GraphAligner [5] rely on\nco-linear matching strategies during both read alignment and graph construction. As a result, they\noften fail to accurately capture non-linear SVs, including duplications, inversions, and nested\nvariants.\nTo substantiate this comparison, we conducted an analysis using minigraph with the GRCh38\nreference assembly and the T2T HG002 paternal assembly, evaluating its performance against\nVACmap on three specific SVs. The results are detailed in Supplementary Note 1: Comparison\nof VACmap with Graph-Based Method. This analysis includes a duplication (chr1:1,249,227-\n1,249,623) that misinterpreted as an insertion by minigraph, an inversion (chr1:92,796,716-\n92,796,839) that minigraph failed to detect, and an inversion with precisely defined breakpoints\n(chr17:5,982,441-5,983,595), in which VACmap outperformed minigraph by reducing the number\nof mismatch bases from 23 to 1. These case studies demonstrate that VACmap’s non-linear chaining\nalgorithm, which integrates linear and non-linear edits within a weighted directed acyclic graph,\neffectively overcomes the co-linear limitations of graph-based methods. This enables more accurate\nsequence alignment and graph construction, particularly in the presence of complex SVs.\nTo reflect these findings, we have added a paragraph to the Discussion section that highlights the\nadvantages of graph-based representations while acknowledging the challenges faced by existing\nalgorithms. This discussion refers to Supplementary Note 1 to support our argument that VACmap\naddresses these limitations, positioning it as a valuable tool for advancing graph-based genomics.\n-3）We thank the reviewer for their detailed and constructive feedback. Below, we address each\ncomment and outline the revisions made to improve the manuscript.\nMain Section Comments\n Lines 61–64: The class of algorithms the authors refer to is commonly known in\nprevious research as seed-and-extend algorithms; stating this would make the\nparagraph more straightforward.\nWe agree that specifying the \"seed-and-extend\" classification would enhance clarity. We\nhave revised the text on lines 73–77 as follow: These subalignments are identified through\na process commonly known as the seed-chain-extend algorithm, tailored for long-read\nmapping. In this approach, seeds (exact matches, such as k-mers) are identified, and a co-\nlinear subset of these seeds is selected to form chains, which are then extended\nbidirectionally until significant differences between the read and reference sequences are\nencountered.\n Lines 70–11: It is unclear to me to what the authors refer when stating that minimap2\nselects ‘the highest-scoring sub-alignment at each step until all alignable regions of\nthe read are covered’. Minimap2 uses a co-linear chaining algorithm [1, § 2.1.1] to\nfind the best match of a set of anchors (which I believe the authors call sub-alignments)\nin the read and in the reference. While a heuristic is applied to only consider a fixed\nnumber of matches, this is not a greedy approach. However, a greedy approach is\napplied when identifying primary chains [1, § 2.1.3].\nWe appreciate the reviewer’s clarification. The statement was intended to describe\nminimap2’s alignment selection process but was oversimplified. We have revising lines\n83–90 as following: For instance, minimap2 employs a co-linear chaining algorithm to\nidentify the set of all possible co-linear chains. It then uses a greedy strategy during primary\nchain selection [18], which determines the optimal subset of subalignments to represent the\nread. This process initializes an empty set Q and iteratively processes subalignments from\nthe highest to the lowest chaining scores: if a subalignment overlaps with a chain in Q by\n50% or more of the shorter subalignment’s length, it is marked as secondary; otherwise, it\nis added to Q, ultimately representing the read with the subalignments in Q [18].\n Lines 71–74: The authors state that NGMLR improves the approach of YAHA by\nutilising dynamic programming. However, YAHA also uses dynamic programming [2,\n§ 2.4].\nWe acknowledge the oversight. Both NGMLR and YAHA utilize dynamic programming,\nand our statement was inaccurate. We have revised lines 90–94 as following: In contrast,\nYAHA adopts a graph-based approach, leveraging Optimal Query Coverage (OQC)\nalgorithm to finds the optimal set of subalignments that cover the length of the query [19].\nNGMLR builds on a similar strategy, enhancing it with a refined scoring function to\nidentify the optimal combination of subalignments with the with the highest joint score [20].\n Lines 99–100: While complex structural variants may be difficult to detect using co-\nlinear chaining with a linear reference genome, the authors should compare their\napproach to graph-based aligners at least to a certain extent.\nWe thank the reviewer for this suggestion. We have added a comparative analysis with\ngraph-based aligners, such as minigraph, in the Discussion section and Supplementary\nMaterials. This comparison highlights VACmap’s advantages in handling complex SVs, as\ndetailed in our response to comment -2 in a prior interaction.\nMethods Section Comments\n Lines 381–382: The authors explain that the acyclicity of the graph generated from\nthe anchors follows from the topological order of its nodes, but the authors do not\ndefine an order relation for the nodes. In particular, is the order determined by e.g. x\nor y only, by all of the elements of an anchor from left to right, or is the subsequence\nlength considered when comparing the co-ordinates (i.e. is (x1, y1, s1, k1) < (x2, y2, s2,\nk2) only if e.g. x1 + k1 < x2)?\nThe order relation is based on the end position in the read sequence, ensuring a topological\norder. We have revised the text to specify that nodes are ordered by their end position in\nread sequence, with edges allowed only if 𝑥 +𝑘 ≤𝑥 +𝑘 .\n𝑖 𝑖 𝑗 𝑗\n Lines 385–393: The equations have been presented as a long list with explanations in\nthe next paragraph. The meaning of some variables (including a, readgap) has been\nomitted. All variables should be defined clearly.\nWe have restructured lines 439–519 to present each equation individually with clear\ndefinitions for all variables (e.g., a, readgap, refgap) in the accompanying text, improving\nreadability and comprehension.\n Line 390: It is unclear to me, is the intention to define a conditionally (for which\nnotation similar to that in e.g. Equation 8 should be used) or is it something else.\nAdditionally I did not see a being used in the other equations.\nThe intention was to define a conditionally, and we thank the reviewer for pointing out the\ninconsistency. We have corrected this equation which assign value to a to use the\nconditional notation and ensured a is applied in the penalty calculation for normal edge\n(Equation 9).\n Line 397: I believe the term hyperparameter is from machine learning and I have not\nseen it in use elsewhere. The authors could call these just parameters.\nWe acknowledge the reviewer’s concern. While “hyperparameter” is common in machine\nlearning, we have replaced it with “parameter” to align with broader bioinformatics usage.\n Line 423: I could not find a definition for the set M.\nWe apologize for the omission. We revised the text as follow: However, we don’t need to\ntest all node[i] because node[j] best predecessor score 𝑆(𝑖) always higher to the current\nhighest score S[k] minus (𝛽+𝛾+20). So for node[j] we define a set M which include\nj\nall nodes who score higher than S[k] minus (𝛽+𝛾+20). Then we only need to test nodes\nin set M to identify the best predecessor of node[j]. This further reduce the time complexity\nj\nto O(hN), where h equal to the average size of set M.\n Line 448: Why are specifically these conditions considered when merging anchors?\nWhat is the significance of the value 19?\nWe thank the reviewer for their insightful question regarding the conditions considered\nwhen merging anchors and the significance of the value 19. The anchor merging process,\nnow detailed in the revised Lines 554-558, anchors are merged if they are overlapping or\nclosely adjacent and share identical sequences in their overlapping regions.\nThe significance of the value 19 bp lies in its role as a maximum limit on the k (k-mer\nj\nlength at node[j] of merged anchors. The merging process can produce larger k values,\nj\nwhich increase the search space and iteration count in the best predecessor computation\nduring the dynamic programming step. To address this computational burden, we\nempirically determined a maximum k of 19 bp. The revised manuscript reflects this\nj\nexplanation, clarifying how the 19 bp cap mitigates the impact of large k values on the\nj\nchaining process.\n Line 473: I was not able to find a URL where the synthetic data is available.\nThe synthetic data generated by VACsim is available at\nhttps://github.com/micahvista/VACmap/vacsim.\nMinor Issues\n The authors publish their implementation under the GNU General Public Licence\nversion 3.0. I would like to point out that all files published under the GPL should\ncontain a licence notice, which currently is not the case. If needed, the authors can\nrefer to the instructions on how to apply the GNU licences for further clarification:\nhttps://www.gnu.org/licenses/gpl-howto.html\nWe thank the reviewer for this observation. We have added the required GPL v3.0 license\nnotice to all source files in the VACmap GitHub repository\n(https://github.com/micahvista/VACmap) and updated the documentation accordingly.\n Ref. 14 only has the URL of the paper in question.\nWe have corrected Reference 14 by including full citation details, now listed as: [14]\nSchloissnig, S., Pani, S., Rodriguez-Martin, B., Ebler, J., Hain, C., Tsapalou, V., Söylev, A.,\nHüther, P., Ashraf, H., Prodanov, T., Asparuhova, M., Hunt, S., Rausch, T., Marschall, T.,\n& Korbel, J. O. (2024). Long-read sequencing and structural variant characterization in\n1,019 samples from the 1000 Genomes Project. bioRxiv (Cold Spring Harbor Laboratory).\nhttps://doi.org/10.1101/2024.04.18.590093.\n Ref. 29 seems incomplete in the sense that at least the journal is missing.\nWe apologize for the error. Reference 29 has been updated to: [29] O'Regan G, deSouza\nRM, Balestrino R, Schapira AH. Glucocerebrosidase Mutations in Parkinson Disease. J\nParkinsons Dis. 2017;7(3):411-422. doi: 10.3233/JPD-171092. PMID: 28598856.\n Line 98, Line 394:: collinearity → co-linearity (or colinearity) We have corrected\nrevised line 118 to read “co-linearity” for consistency with standard usage.\n Line 384: Formulae are more commonly called equations in this context.\nWe have replaced “Formulae” with “equations” to align with common scientific\nterminology.\n Line 446: O(log^N) → O(log N)\nWe thank the reviewer for catching this typo. We have corrected revised line 549 to read\nO(log N).\n Lines 528–529: All the three links are the same.\nWe apologize for the duplication. The correct links have been updated:\ns3://giab/data/AshkenazimTrio/HG002_NA24385_son/, s3://giab/data/AshkenazimTrio/\nHG003_NA24149_father/, and s3://giab/data/AshkenazimTrio/ HG004_NA24143_mother\n/.\nReviewer #1 (Remarks on code availability):\nThe Conda environment could not be installed on my macOS laptop. I currently do not have access\nto a suitable Linux workstation for running the complete experiments. Having a very simple set of\ntest data would have been useful.\nResponse:\nWe thank Reviewer #1 for their feedback regarding the code availability and installation process.\nWe acknowledge the challenge with installing the Conda environment on macOS and clarify that\nVACmap is currently designed to run exclusively on Linux systems due to its dependency on\nspecific Linux-based libraries and optimization requirements. We apologize for any inconvenience\nthis may cause and are exploring future enhancements to broaden platform compatibility.\nTo assist with verifying the installation, we have uploaded a small test dataset, including a\nreference.fasta and a read.fasta, to the repository. This dataset is intended to allow users to confirm\nthat the VACmap installation is successful by running a basic alignment test. We hope this addition\naddresses the reviewer’s suggestion for a simple test set and facilitates accessibility for those with\nLinux environments. We appreciate the reviewer’s input and welcome further suggestions to\nimprove usability.\nReviewer #2 (Remarks to the Author):\nThe authors are introducing VACmap, a new long-read mapping tool that significantly\nenhances the representation of all types of SVs. The tool's innovative approach employs a non-\nlinear chaining algorithm to connect subalignments, resulting in markedly improved detection\nof complex structural variants and more accurate quantification of copy number variations.\nThe extensive supplementary materials, including tables and figures, provide compelling\nevidence that substantiates the mapping enhancements that lead to improved structural\nvariant calling. The work addresses a highly significant area in clinical genetics, as structural\nvariants are increasingly recognized as important diagnostic targets in medical practice. I\nparticularly appreciate the performance in the case of the gene conversion example.\nResponse: We sincerely appreciate Reviewer’s thorough and constructive feedback on our\nmanuscript introducing VACmap, a novel long-read mapping tool designed to enhance the detection\nof structural variants (SVs). We are grateful for the reviewer’s recognition of the tool’s innovative\nnon-linear chaining algorithm, its improved performance in detecting complex SVs and quantifying\ncopy number variations, and its relevance to clinical genetics. Below, we address each comment and\noutline the revisions made to strengthen the manuscript.\nComments:\n-1）VACmap is evaluated against 3 other long-reads aligners using both simulated data and\nGIAB samples. Their methodology involves mapping long reads to the T2T human reference\ngenome followed by structural variant calling with SVIM. I suggest expanding the comparison\nto include more contemporary aligners such as LRA or LAST, while mentioning in the\nmanuscript that NGMLR is no longer actively maintained. The manuscript's impact would be\nsubstantially enhanced by incorporating an additional structural variant caller like CuteSV\nor Sniffles2, both of which have demonstrated superior performance in recent benchmarks\n(like https://www.nature.com/articles/s41598-024-56604-2 Benchmarking long-read aligners\nand SV callers for structural variation detection in Oxford nanopore sequencing data).\nResponse:\nWe thank the reviewer for suggesting the inclusion of more contemporary aligners such as LRA, as\nwell as noting that NGMLR is no longer actively maintained. We have expanded the comparison in\nthe Results section to include LRA alongside minimap2, NGMLR and Winnowmap2, using GIAB\nsample, followed by SV calling with SVIM and cuteSV, referencing the recent study\n(https://www.nature.com/articles/s41598-024-56604-2). The updated analysis demonstrates\nVACmap’s superior performance across these aligners. We have also added a note in the Results\nsection stating that NGMLR is no longer actively maintained.\n-2）Furthermore, given that pseudogenes can originate from specific SV events, it would be\nvaluable to examine VACmap's capacity to identify pseudogenes, particularly in comparison\nto duplomap, aligner dealing well with segmental duplication alignment. A compelling case\nstudy would be the STRC gene (implicated in hearing loss) and its adjacent pseudogene, which\nshares over 99.6% sequence identity with STRC.\nResponse:\nWe thank the reviewer for the insightful suggestion to evaluate VACmap’s capacity to identify\npseudogenes, particularly in comparison to Duplomap, which excels in aligning reads within\nsegmental duplication regions, and for recommending a case study on the STRC gene. The STRC\ngene, associated with mild-to-moderate autosomal recessive hearing loss, is notoriously difficult to\ndetect due to its location within a tandem duplication region and the presence of its highly\nhomologous pseudogene, STRCP1, which shares >99.6% sequence identity. This high similarity\noften leads to misalignments in standard aligners, complicating SV detection.\nTo address this, we analyzed GIAB sample (NA19240) using VACmap, with alignment results\nvisualized in Figure 7b. VACmap successfully identified a heterozygous deletion involving the loss\nof one copy of the CKMT1B-STRC-CATSPER2 gene cluster. The alignment tracks in Figure 7b\ndemonstrate VACmap’s ability to accurately map reads across this complex region, with clear\ndelineation of the deletion boundary (highlighted in red bar). In contrast, other aligners—NGMLR,\nLRA, Winnowmap2, and minimap2—failed to detect this deletion. Duplomap, despite its strength\nin segmental duplications, also missed the deletion, likely due to its reliance on minimap2’s internal\nmapping, which tends to avoid splitting reads in highly identical regions, resulting in incomplete\nrepresentation (see Duplomap track in Figure 7b).\nThis case study underscores VACmap’s superior performance in resolving pseudogene-associated\nSVs, attributed to its non-linear chaining algorithm, which effectively handles high sequence\nidentity and tandem duplications. We have added a detailed discussion of this analysis in the Results\nsection (Lines 344–353).\n-3) Concerning the VACmap software implementation: How do the mapping quality scores vary\nacross the different mapping modes (H, L, S, R parameter settings)? I recommend including more\nspecific guidance on selecting appropriate modes for various applications. Additionally, please\nspecify in the manuscript which mapping mode was employed to generate the reported F1 scores.\nResponse:\nWe sincerely thank the reviewer for their concern regarding VACmap’s software implementation,\nspecifically how mapping quality (mapq) scores vary across the different mapping modes (H, L, S,\nR parameter settings), and for their recommendation to provide specific guidance on selecting\nappropriate modes for various applications. We have addressed these points in the revised\nmanuscript, ensuring clarity and practical utility for users.\nThe map quality scores in VACmap are computed using a formula adapted from minimap2, as\ndetailed in the “Map Quality Calculations” subsection in the Results section:\n:\n𝑚\n𝑚𝑎𝑝 𝑞𝑢𝑎𝑙𝑖𝑡𝑦 =40×(1−𝑓2/𝑓1)×min (1, )× ln ( 𝑓1). (1)\nand m is the number of anchors on the highest-scoring chain. Since mapq primarily depends on the\n𝑓2\nratio where a larger difference between f1 and f2 yields a higher mapq—and the number of\n𝑓1\nanchors m, the scores remain consistent across modes for comparable alignments. The mapping\nmodes (H, L, S, R) adjust alignment strategies (e.g., variation penalties and divergence thresholds)\nto optimize for specific data types or genomic regions, but these adjustments do not significantly\nalter the relative difference between f1 and f2.\nTo provide specific guidance on selecting appropriate modes, we have added a new subsection,\n“Mode Selection,” in the Methods section. Mode L, with a high variation penalty (β=59) and a low\ndivergence threshold (0.1), is recommended for high-accuracy long-read data (e.g., PacBio HiFi).\nMode H, using a smaller penalty (β=40) and a higher divergence threshold (0.2), is preferred for\nnoisy long-read data (e.g., ONT). Mode S, with a smaller penalty (β=30) and a high divergence\nthreshold (0.5), is designed for recover more small-scale rearrangements or SVs located in regions\nwith high homology. Mode R is an experimental mode intended solely for testing purposes, with a\nfixed variation penalty (VP =30) and a divergence threshold of 0.5. Due to its developmental nature,\ni,j\nMode R is not recommended for normal usage, and users should instead rely on Modes L, H, or S\nfor standard applications.\nFor parameter setting of VACmap to generate the reported F1 scores, we use Mode L to align PacBio\nHiFi sequence data and Mode H to align PacBio CLR and ONT sequence data. Which mentioned\nin Method section (Lines 618-619, 639-641).\nResponse to Reviewer #2 Comment on Minor Textual Errors\nWe sincerely thank Reviewer #2 for their meticulous review and for identifying several textual\nerrors in our manuscript. We have addressed each concern, ensuring accuracy and clarity in the\nrevised manuscript. Below, we provide a point-to-point response to each issue raised:\n Line 121: The text states, \"As illustrated in Figure 2a, VACmap's alignment data markedly\nenhanced SVIM’s[21] ability to detect complex SVs in simulated data from PacBio CLR,\nPacBio HiFi, and ONT, with F1 score improvements ranging between 31.9 to 68.9 percent,\"\nbut SVIM is not shown or mentioned in the Figure 2a caption.\nResponse: We have revised the Figure 2a caption to explicitly mention SVIM and its role\nin the analysis, ensuring consistency with the text. The caption now reads: \"Figure 2a: The\nprecision, recall, and F1 scores of SVIM’s complex SV detection performance in all\nchromosomes and the repetitive region across different read depth and sequence\ntechnologies.\"\n Line 210: The typo \"arcuate\" should be corrected to \"accurate.\"\nResponse: We have corrected the typo on Line 240, replacing \"arcuate\" with \"accurate\" to\nensure proper terminology.\n Line 212: The phrase \"a gene is involved\" should have \"is\" removed.\nResponse: We have revised the phrase on Line 241 to \"a gene involved\" by removing \"is,\"\nimproving grammatical accuracy.\n Lines 215–216: The sentence \"It’s evidence in the size standard deviation of the inversion\nin SVIM calls is 254.2 and 2066.4 for VACmap and NGMLR alignment, respectively\"\nrequires rephrasing for clarity.\nResponse: We have rephrased the sentence on Lines 244–246 to: \"The size standard\ndeviation of inversion calls by SVIM is 291.4 for VACmap alignments and 2066.4 for\nNGMLR alignments, respectively,\" enhancing clarity and grammatical correctness.\n Line 218: The text \"the SVIM quality score for this inversion is 13, 0 for VACmap,\nNGMLR alignment\" should include \"and\" between \"VACmap\" and \"NGMLR\" and add\n\"will be discarded.\"\nResponse: We have updated Lines 247-248 to: \"The SVIM quality score for this inversion\nis 14 and 0 for VACmap and NGMLR alignments, respectively, and will be discarded,\"\nadding \"and\" and \"will be discarded\" as requested.\n Figure 4 Caption: The caption incorrectly states \"4f is 4e,\" while the reference in the text\nis correct.\nResponse: We have corrected the Figure 4 caption by updating the panel label from \"4f\" to\n\"4e,\" ensuring consistency with the text references and resolving the discrepancy.\nThese revisions address all minor textual errors identified by the reviewer, improving the\nmanuscript’s accuracy and readability. We appreciate the reviewer’s attention to detail, which has\nhelped enhance the quality of our submission"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I find that the authors have resolved the issues that I pointed out. In particular, the Variant-aware Chaining Algorithm section is now substantially easier to follow. I only have some minor comments but I recommend the manuscript for publication. I successfully ran VACmap with the provided test data.\n- In the Variant-aware Chaining Algorithm section, I would use “node i” and “node j” to refer to the nodes in question, instead of ”node[i]” and “node[j]”.\n- Similarly, since in Equation 14 the score of node i is defined as S(i), I would use the notation in question instead of S[i].\n- Considering Equation 1, the authors could clarify whether the condition is tested for all j or only the smallest j that satifies it.\n- Additionally I am interested in the overall time complexity of the algorithm (average case or otherwise), taking the sorting steps into account.\n- Figure 4e: Alignment score growth slowly → Alignment score grows slowly\n- Line 299: gene → genes\n- Line 380: we show existing → we show that existing\n- Line 423: collinear → co-linear\n- Lines 453–454: To determine edge type → To determine the edge type\n- Line 490: equation 12 and 13 → Equations 12 and 13\n- Lines 501–502: the equation 14 → Equation 14\n- Line 554: can merge → can be merged\n- Lines 555–558: please consider rewriting the last sentence of the paragraph.\n- Line 571: hyperparameter → parameter\n- Lines 572 and 573: maxgap → maxgap"}, "author_response": "Dear Editor and Reviewer,\nWe sincerely thank you for your time and effort in reviewing our manuscript, titled \" VACmap: An\nAccurate Long-Read Aligner for Unraveling Complex Genomic Rearrangements\". We are deeply\ngrateful for the insightful and constructive comments provided by the reviewer, which have\nsignificantly improved the clarity, rigor, and impact of our work.\nIn the following, we present one-by-one responses to the reviewer’s comments.\nReviewer #1 (Remarks to the Author):\nI find that the authors have resolved the issues that I pointed out. In particular, the Variant-aware\nChaining Algorithm section is now substantially easier to follow. I only have some minor comments\nbut I recommend the manuscript for publication. I successfully ran VACmap with the provided test\ndata.\nResponse: We’re glad the previous revision improved the Variant-aware Chaining Algorithm section.\nIn the new revised manuscript, we have carefully addressed each of your comments and suggestions.\nPlease find our detailed, point-by-point responses below.\nComments:\n-1）In the Variant-aware Chaining Algorithm section, I would use “node i” and “node j” to refer to\nthe nodes in question, instead of ”node[i]” and “node[j]”.\n- Similarly, since in Equation 14 the score of node i is defined as S(i), I would use the notation in\nquestion instead of S[i].\nResponse:\nWe have updated the notation throughout the Variant-aware Chaining Algorithm section to use\n“node i” and “node j” and use S(i) and S(j) instead of S[i] and S[j] in the relevant text and equations.\n-2）Considering Equation 1, the authors could clarify whether the condition is tested for all j or only\nthe smallest j that satifies it.\nResponse:\nWe have revised the text and clarify this condition tests for all node j that satisfy it.\n-3）Additionally I am interested in the overall time complexity of the algorithm (average case or\notherwise), taking the sorting steps into account.\nResponse:\nWe have revised the manuscript to add the text about the overall time complexity as follows (added\nto the Optimization subsection): In practice, we maintain an array of nodes in descending order by\ntheir scores using binary search. Then the overall time complexity of the algorithm is O(kN + N log\nN). The average value of k for HiFi, ONT, and CLR data in the first round of chaining is 21, 33, and\n44, respectively.\n-4） Figure 4e: Alignment score growth slowly → Alignment score grows slowly\nResponse:\nWe have corrected the caption for Figure 4e to \"Alignment score grows slowly.\"\n-5）Line 299: gene → genes\nResponse:\nWe have corrected \"gene\" to \"genes\".\n-6) Line 380: we show existing → we show that existing\nResponse:\nWe have inserted \"that\": \"we show that existing...\"\n-7) Line 423: collinear → co-linear\nResponse:\nWe have changed \"collinear\" to \"co-linear\" for consistency with earlier usage.\n-8) Lines 453–454: To determine edge type → To determine the edge type\nResponse:\nWe have added \"the\": \"To determine the edge type...\"\n-9) Line 490: equation 12 and 13 → Equations 12 and 13\nResponse:\nWe have capitalized and pluralized to \"Equations 12 and 13\".\n-10) Lines 501–502: the equation 14 → Equation 14\nResponse:\nWe have corrected to \"Equation 14\".\n-11) Line 554: can merge → can be merged\nResponse: We have changed \"can merge\" to \"can be merged\".\n-12) Lines 555–558: please consider rewriting the last sentence of the paragraph.\nResponse:\nWe have rewritten the last sentence of the anchor merging paragraph for improved clarity: \"To avoid\nproducing large anchors—which could expand the search space during best predecessor\ncomputation—we limit the maximum output size (default: 19) for merged anchors, thereby reducing\ncomputational demand in the chaining process without compromising alignment quality.\"\n-13) Line 571: hyperparameter → parameter\nResponse:\nWe have replaced \"hyperparameter\" with \"parameter\", consistent with previous revisions.\n-14) Lines 572 and 573: maxgap → maxgap\nResponse: We have not made changes here, as the formatting adheres to the editor’s requirement.\nWe believe these changes fully address the reviewer's minor comments and further polish the\nmanuscript. Thank you again for your valuable input."}]

Decision letter