Use of 3D chaos game representation to quantify DNA sequence similarity with applications for hierarchical clustering

The Chaos Game Representation (CGR), especially in 2D, has been frequently utilized as an alignment-free method to capture information about DNA sequences (Vinga et al., 2012; Almeida et al., 2001; Hoang et al., 2016; Huang and Shi, 2010). In 2D CGR, each vertex of the unit square is assigned to one of the four possible nucleotides, e.g., A = (0, 0), C = (0, 1), G = (1, 1), T = (1, 0). The representation starts at the center of the square, at $\left(\frac{1}{2},\frac{1}{2}\right)$. For each nucleotide in the DNA sequence of interest, the subsequent CGR value is positioned halfway between the previous CGR value and the corresponding corner of the nucleotide. An example of how a 2D CGR is iteratively built for a small subsequence is presented in Figure 1.

CGR provides a computationally efficient, lossless representation, wherein all possible sequences have a unique representation. Theoretically, an entire sequence can be recovered from its final CGR coordinate, with the resolution of the final coordinate being the only limiting factor as to how much of the sequence can be recovered. CGR has numerous useful properties: the 2D CGR plane is a scale-independent generalization of Markov chain probability tables (Almeida et al., 2001); it can be used to identify local alignment points in DNA sequences (Almeida et al., 2001; Joseph and Sasikumar, 2006); and the entropic profiles of sequences can be identified through estimations of local densities (Vinga and Almeida, 2007).

However, a thorough exploration of the entire trajectory and structure of a chaos game representation to quantify global similarity between entire sequences has been mostly unexamined. A point-by-point comparison of the trajectory is typically restricted to local comparisons because it requires sequences be of equal length and/or that the sequences be aligned (similar to traditional, alignment-based sequence comparison methods) prior to comparison. Even in the circumstance where all DNA sequences to be compared are identical in length, the 2D CGR may force some nucleotides to be arbitrarily farther apart in space in a way that has potential downstream consequences. For example, consider using a Euclidean distance to compare sequences of CA, GA, and TA. The trajectory of CA would be $(0.25,0.75)\rightarrow(0.125,0.375)$; the trajectory of GA would be $(0.75,0.75)\rightarrow(0.375,0.375)$; the trajectory of TA would be $(0.75,0.25)\rightarrow(0.375,0.125)$. At every step of the trajectory, CA is closer to GA than it is to TA because the C vertex is positioned closer the G vertex in the 2D-plane; this is despite the fact that both GA and TA are a single nucleotide different from CA. So the distance between any two sequences is not only a function of sequence similarity but also depends on the spatial representation of the 2D CGR. As such, the inflated dissimilarity between the CA and TA (or any general pair of sequences) is not grounded in biological reality.

Many applications of 2D CGR handle unequal sequence lengths using frequency counts (FCGR) at various grid resolutions or using a function of the CGR points to calculate a distance between sequences and find local homology between them (Almeida et al., 2001; Joseph and Sasikumar, 2006). From these studies, the use of the 2D CGR to form local-level comparisons of DNA sequences is well established. But its use to determine a global similarity based on point-level comparisons or the entirety of its trajectory has not been fully realized. The method proposed by Hoang et al. (2016) is one attempt to perform a global similarity comparison by examining the CGR in the frequency domain, but this approach also suffers from the same issues noted in Yin et al. (2014a) where deletion mutations strongly influence the final distance metrics.

In this paper, we propose an alignment-free sequence analysis method that utilizes a 3D chaos game encoding to represent DNA sequences and quantifies similarity and/or distance between sequences based on the resulting 3D trajectory by using shape-comparison techniques. Our method requires no sequence padding or stretching that distorts the information embedded. Instead, we focus on a shape-similarity comparison of the overall trajectory of sequences to evaluate their distances. Our method is unique in that distances generated by deletion mutations vs substitution mutations are often comparable in a way that we hope links the estimated distance metric to the number of evolutionary events.

In the Methods section, we: 1) describe the construction of the 3D CGR, 2) discuss similarity/distance measures used to compare the 3D CGR, 3) discuss the phylogenetic algorithm used to test downstream analysis of CGR-based distances, and 4) discuss the simulations used to evaluate the effect of different mutation types on distance under the proposed methodology. In the Data section, we describe the data we used to test our proposed method. In the Results section, 1) we show how the method performs for phylogenetic analysis, 2) demonstrate the relative stability of estimated distances across a wide range of parameter settings, and 3) through simulations, show how the method performs under different types of mutations. In the Application section, we apply our method to a set of SARS-CoV-2 sequences collected between 2020 and early 2023 to demonstrate its effect on a substantially larger dataset as well as to engage in a criticism of the literature on alignment-based methods and the use of phylogenetic trees as the standard by which we judge their performance. Finally, in the Discussion section, we summarize our findings and discuss the performance of our method relative to other existing alignment-free sequence comparison tools in the AFproject (Zielezinski et al., 2019) which provides objective benchmarking against known references.

The proposed 3D CGR forms a Sierpinksi tetrahedron where each of the four possible nucleotides are vertices on the tetrahedron. The CGR vertices are coded as coordinates in three-dimensional space in the following manner: A = $\alpha\left(-1,-1,-1\right)$; T = $\alpha\left(1,1,-1\right)$; C = $\alpha\left(1,-1,1\right)$; G = $\alpha\left(-1,1,1\right)$, where multiplier $\alpha$ is some constant. In this paper, we assume a default of $\alpha=\frac{1}{\sqrt{3}}$ to set the Euclidean distance between each vertex to the origin equal to 1.

This specification ensures that each nucleotide is equidistant from every other possible nucleotide and the angle between any two nucleotides through the origin is identical.¹¹1Chang et al. (2003) offers an alternative 3D encoding of DNA sequences, not specific to chaos games, with the same characteristics: A = $(0,0,1)$; T = $\left(0,\frac{2\sqrt{2}}{3},-\frac{1}{3}\right)$; C = $\left(-\frac{\sqrt{2}}{3},\frac{\sqrt{6}}{3},-\frac{1}{3}\right)$; G = $\left(-\frac{\sqrt{2}}{3},-\frac{\sqrt{6}}{3},-\frac{1}{3}\right)$ Starting with the center at $X_{0}=(0,0,0)$, the CGR representation $(X_{n})=(x_{n},y_{n},z_{n})$ for a DNA sequence $s_{1}s_{2}...s_{n}$ where $s_{i}\in\{A,C,G,T\}$ is built iteratively in the following manner:

where $W$ maps nucleotides to vertices as specified above. This algorithm is identical to that of the 2D CGR but with a new 3D system of coordinates. A 3D CGR of the mitochondrially encoded NADH:ubiquinone oxidoreductase core subunit 4 (MT-NADH) gene for the Macaca fascicular (a species of old-world monkey) is shown in Figure 2.

In the CGR of a DNA sequence, each nucleotide’s corresponding coordinate is determined by its position in the sequence as well as by all preceding nucleotides. This encoding method preserves uniqueness and each DNA sequence will always produce a distinct CGR coordinate. However, an interesting feature of CGR is that matching substrings within two different sequences will have coordinates that begin to converge and occupy the same position as the length of the match increases; the longer the match, the closer in proximity their 3D CGR regardless of all preceding elements in the sequences.

Take, for example, the small sequences ATCAGGCAG AND TGTAGGCAG that have mismatching nucleotides in the first 3 positions, but share the same subsequence in their remaining 6 positions. The 3D CGR of these sequences, shown in Table 1, illustrate how although the coordinates are unique to each sequence, the matching subsequences have coordinates that are very similar and whose differences begin to disappear the longer the sequence match. This convergence in coordinates happens even when the mismatching prefixes are of unequal length. Because the trajectory of identical subsequences will converge, subsequences in the 3D CGR are, in a way, self-aligning. This makes this representation of DNA sequences amenable to alignment-based shape comparison techniques.

In order to compare DNA sequences encoded in 3D CGR, we look to distance/similarity measures that have been developed to compare shapes and volumes for inspiration. Several approaches for comparing 3D conformations, in particular with application to molecular structure, have been developed for performing shape similarity quantification.

Phylogenetic trees are built using distances derived from the methods detailed in Section 2.2 to test their performance. For comparison, two baseline methods were also implemented. For the first baseline method, sequences were aligned using Clustal Omega alignment using the msa package in R (Bodenhofer et al., 2015) and distances for the aligned sequences were calculated using pairwise polymorphism p-distances as implemented by the dist.gene() in the ape package (Paradis and Schliep, 2019).²²2Evolutionary models such as the Jukes-Cantor model were not used because they are generally not amenable to deletion mutations. For the second baseline method, distances between the sequences were calculated using k-mer distances using the kmer package. We optimized k-mer parameters by selecting those that best reproduced the phylogenetic trees built under Clustal Omega and p-distances. We use $k=6$ and Euclidean k-mer distances. We avoid the default distance metric, proposed by Edgar (2004), as it produced non-monotone distances and reversals in the dendrograms, and generally resulted in poor phylogenetic reconstruction.³³3Other baseline approaches that were tested and excluded due to poor performance include: Nucleotide Amino Acid K-mer Vector (Bao et al., 2022), Extended Natural Vector method (Pei et al., 2019), and Positional Correlation Natural Vector method (He et al., 2020). Both the Extended Natural Vector method and the Positional Correlation Natural Vector method performed considerably worse than the presented baseline methods; the Nucleotide Amino Acid K-mer method was comparable to the k-mers method and not sufficiently different to warrant deviating away from the more general and well-known k-mers approach..

Phylogenetic trees are built using the neighbor-joining method as implemented by the bionj()⁴⁴4bionj() uses an improved implementation of the NJ algorithm by Gascuel (1997). function in the ape package (Paradis and Schliep, 2019). This agglomeration method was selected because it is a popular and efficient method. Because it does not assume a constant rate of evolution across all lineages, it is applicable to a wide range of biological datasets and can be applied to the different types of DNA sequences we examine. Other agglomeration methods were tested as well but, in some instances, the reference methods performed worse than the 3D CGR methods under those alternative agglomeration approaches; to ensure that our method is robust, we always opt for modeling choices that benefit the reference methods. Ultimately, however, the choice of agglomeration method does not impact the substantive conclusions made regarding the 3D CGR method as the CGR method performed equally well across different agglomeration methods.

The effects of deletion mutations and substitution mutations on the 3D CGR are explored using the coding sequences (CDS) of the first exon of the $\beta$-globin gene for the chimp which consists of 105 base pairs. Point mutations are introduced to this parent sequence, increasing from 2 mutations to 76 total mutations in intervals of 2. One child lineage receives only substitution mutations and the other receives only deletion mutations with each subsequent generation receiving two more mutations than the previous; both substitution and deletion mutation locations correspond to the exact same locations on the original parent sequence and once a location has been mutated it cannot be mutated again. Bandwidths for the volume intersection method ranged from 0.0001 to 0.0009 in intervals of 0.0001, 0.001 to 0.009 in intervals of 0.001, and 0.01 to 0.1 in intervals of 0.01. The number of bins used for histograms of the shape signature methods, at the low end, is set to equal 5% of the average sequence length across all lineages and, at the high end, is equal the average sequence length.

Code to produce these analyses (as well as under alternative agglomeration methods) is available on https://github.com/tomato-stats/3D_CGR/.

For the CDS of the $\beta$-globin gene, a heatmap of the distances generated using multiple sequence alignment (MSA) method Clustal Omega and p-distances, k-mer distances, intersection of volumes for Gaussian spheres in 3D CGR representation, and shape signature are provided in Figure 5. For all heatmaps represented, each organism’s distance to itself is excluded from the plot and all remaining distances are scaled to be between 0 and 1; color scales are chosen to improve color diversity within each method and to better highlight the similarities and differences between methods. The resulting phylogenetic trees built from those distances are shown in Figure 5.

The heatmaps show that the alignment-based method results in a smaller distance between the rabbit and the human than between the chimp and the human. A closer look at the rabbit, human, and chimp sequences in Figure 6 reveals why this is the case: there are nine substitutions between the human and the rabbit; in contrast, the chimp sequence is identical to the human sequence with the exception of an insertion of 13 additional nucleotides. The total number of positions with mismatched nucleotides is smaller between the human and rabbit than between the human and chimp.

In contrast, all of the alignment-free methods place the chimp closest to the human. Since the 13-nucleotide insertion at the end of the chimp sequence is likely to be a single mutation event, the total number of mutation events between the human and chimp is likely fewer than the number between the human and rabbit. In this scenario, it is not unjustified for the alignment-free methods to quantify the distance between the human and chimp as less than the distance between the human and rabbit.

Both the Clustal Omega and k-mer reference methods produce a similar tree topology under the neighbor-joining agglomeration method despite quite disparate distance estimates. Although none of the shape-based approaches perfectly reproduce the trees built by the reference methods, they do produce near-identical leaf-node clusterings with the ungulates (bovine and goat) clustered together, the primates (human, chimp, and gorilla) clustered together, and the rodents (mouse and rat) clustered together. The lemur’s total distance from all other organisms is the highest under both the volume intersection and shape signature methods, which is not concordant with the reference methods. Additionally, placement varies under the shape-based approaches; under the chosen parameter settings, the lemur clusters with no other sequence under the volume intersection method and clusters with the rabbit in the shape signature method.

While it is technically possible to select parameter settings such that the both 3D CGR methods reproduce the general topology of the reference methods, what remains constant is that: 1) the lemur’s clustering position tends to be less stable relative to other organisms and 2) the lemur is consistently amongst the set of most dissimilar organisms, regardless of parameterization, rivaled in total distance only by the oppossum (a marsupial) and the gallus (the only non-mammal organism in the set).

The mutations observed in the lemur are uniquely different from those observed in the other organisms. Firstly, the types of mutations observed in lemurs tends to be different from those of other organisms. Existing examinations of the lemur’s $\beta-$globin gene have found that the lemur lineage of the $\beta$-globin gene contains a significant excess of non-synonymous mutations over the neutral mutation rate, especially in the first 13 codons of the first exon (Harris et al., 1986), which we confirm. In contrast, mutations observed in other organisms are primarily synonymous mutations and/or transition mutations.

Secondly, a large proportion (9 out of 22) of the lemur’s mutations are private mutations; the locations of its mutations are frequently disparate from those of other organisms and occur where all other sequences match the consensus. Excluding from consideration the chimp’s insertion of 13 base-pairs at the end of its sequence, the only other organisms with private mutations are the gallus with 3 out of 22 substitution mutations being private and the opossum with only 1 out of its 23 substitution mutations being private. So when the lemur has a mutation, it is often in a location that is unusual compared to all other sequences. These singleton mutations result in the lemur being considered more distant under the shape-based methods due to the fact that its 3D CGR trajectory is distinct because it diverges at different points than the other organisms.

The 3D CGR methods are picking up on more than just a mismatch in nucleotides; they are also implicitly capturing the location of the mismatch and identifying it as being rare. Since the lemur’s mismatches tend to occur in locations that are only ever mutated in the lemur, this further separates the lemur from the rest of the organisms. A number of potential explanations have been suggested for why the lemur’s mutations are so disparate from sequences such as that of the human and rabbit which have mutation levels in line with a neutral evolution model: adaptive evolution, gene conversion, and relaxed selective constraints (Harris et al., 1986). The 3D CGR captures the contrasting evolutionary pressures that the lemur’s sequence experienced. Whether this is a desirable feature in a distance metric depends on whether a researcher is interested in assigning greater distances to rarer mutations.

A final important note is that that both the Clustal Omega p-distances as well as the k-mer distances required an agglomeration method based on reducing total tree length such as the neighbor-joining algorithm in order to perform well. Under both the unweighted pair group method with arithmetic mean (UPGMA) and the weighted pair group method with arithmetic mean (WPGMA), the sequence-alignment approach placed the rabbit between the chimp and the human/gorilla cluster. Under both UPGMA and WPGMA, the k-mer distance approach consistently resulted in the goat/bovine cluster being placed between the rat and the mouse; increasing the k-mer length up to 10 did not resolve the issue. None of these issues were observed with either of the shape signature methods whose leaf-node clusterings remained robust and consistent with expectations. This is noteworthy because the k-mer distances look very similar to the volume intersection distances. Yet the subtle differences between the two are sufficient to ensure that the volume intersection tends to perform better. This is indicative that the volume intersection approach is capturing additional information that is not captured by the k-mers approach.

For the MT-NADH gene, a heatmap of the distances generated using MSA with Clustal Omega and p-distances, k-mer distances, intersection of volumes for Gaussian spheres in 3D CGR representation, and shape signature are provided in Figure 7. The resulting phylogenetic trees built from those distances are shown in Figure 8.

The distributions of distances are quite varied amongst the different approaches. But although the magnitude of the distances between the methods can be quite different, the patterns of distances are somewhat similar across the methods. For example, the darker red blocks for most pairwise comparisons with Tarsisus syrchita, Lemur catta, and Saimiri sciureus indicate that these primates tend to be the most different from all others regardless of method chosen. And the blue block of old world primates also appears fairly consistently across the different methods. In general, the alignment-free methods tend to consider the old world primates to be much more distant from the hominoids than the alignment-based method. For the k-mers approach, this is likely due to the fact that multiple single-nucleotide substitutions have much larger downstream consequences for k-mer counts. For the volume intersection approach, the disparity is similar as multiple single-nucleotide substitutions in close, but not adjacent positions, can repeatedly alter the 3D CGR trajectory that, when paired with a very small bandwidth, can minimize the calculated volume overlap/similarity.

All methods depicted in Figure 8 correctly cluster the old world primates together and the five hominoid species together. Both the k-mer distance and the distances derived from the volume intersection method produce a tree topology that closely resembles the one built under alignment and p-distances. The shape signature distances also produce very similar trees with the exception of the placement of the Lemur Catta as, unlike the other methods, it is not clustered with the Tarsisus syrichta. However, as there remains some differences of opinion regarding whether the tarsiers truly belong under the group of strepsirrhines/prosimians (tarsiers do not have moist noses while all other strepsirrhines/prosimians do), the difference in the trees appear to reflect ongoing discussions regarding the correct taxonomy of primates. The taxonomy of the primates is likely to be revised in subsequent years as more DNA data is collected and sequenced and so the topology under the shape signature may yet still turn out to accurately reflect biological evolution.

For the synthetic DNA sequences, we provide some ground-truth information in Figure 9. In the first panel of Figure 9, we provide the exact number of separating mutation events based on how the sequences were constructed. For instance, the number of mutation events separating $\Delta$50NT from $+$100NT/50:149, a sequence with a 100-contiguous nucleotide insertions, is exactly 51; one mutation to return $+$100NT/50:149 to the original parent sequences and 50 mutations to $\Delta$50NT. The number of mutation events separating the $\Delta$50NT from the sequences with transpositions (e.g., $\pm$ 50NT/50:99 $\rightarrow$ 150:199) is 51 if the transposition is counted as a single mutation and 52 if the deletion and insertion events are each counted separately; for simplicity, we counted a transposition as 1.5 mutation events. The sequences with the greatest number of separating events are $\Delta$50NT and $-$50NT; both descend from the parent sequences with 50 mutations each and, thus, there are 100 total mutation events separating the two from one another.

In the second panel of Figure 9, we provide the theoretical number of mutation events observable between any pair of sequences under ground-truth alignment. When examining any pair of sequences without contextual knowledge of the parent sequences or information about how those sequences came about, the number of observable differences must be less than or equal to the number of mutation events. For instance, by construction, there are exactly 100 mutation events separating $\Delta$50NT and $-$50NT (50 mutation events for each sequences to return to the original parent sequence). However, in the absence of any context or information about the parent sequence and, given that the location of deletions/substitutions are identical, it would appear that there are only 50 mutation events separating those sequences when those two sequences are correctly aligned. For concreteness, consider the toy example where we compare the sequences ACG and A-G (where the hyphen denotes a deletion mutation), both mutated from their shared ancestor ATG. These two sequences are separated by exactly two mutation events–one substitution event and one deletion event. However, without information about the shared ancestor, we only observe a single difference between ACG and A-G at the second nucleotide position. For comparisons between sequences such as $\Delta$50NT and one of the sequences with 100-contiguous nucleotide deletions, the number of observable differences is 51 if none of the substitution mutations occur within the deleted region and fewer than 51 otherwise. This second heatmap essentially serves as an upper-bound on the number of observable mutations under correct pairwise alignment.

In the third and fourth panel of Figure 9, we perform MSA of all sequences manually, using knowledge of each sequence’s construction. Due to disparate sequence lengths, multiple gap insertions are necessary to align sequences. The original sequences ranged from a minimum of 693 base-pairs to a maximum of 893 base-pairs. The aligned sequences resulted in all sequences being 993 base-pairs long. The third panel heatmap shows the number of mismatching nucleotides after alignment, but without penalizing mismatches due to gaps in the sequences. The fourth panel heatmap shows the same, but with penalties for gaps in the sequences. The fourth panel corresponds to p-distances under a perfect MSA method.

For comparison to the ground-truth data of the synthetic sequences, a heatmap of the distances generated using Clustal Omega with p-distances, k-mers, intersection of volumes for Gaussian spheres in 3D CGR representation, and shape signature are provided in Figure 10. Though not especially meaningful, the resulting phylogenetic trees built from those distances are shown in Figure 11 for completeness.

From comparing the ground-truth plots and the distance plots, we can see that the Clustal Omega method with p-distances produces distances that are similar to the fourth ground-truth plot, p-distance after alignment (not ignoring indels). The shape-based methods produce distances similar to the heatmap for the number of pairwise-observable difference. The k-mers method produces distances that are also similar to the number of pairwise-observable differences, but the resemblance is weaker than for the shape-based methods.

For all distance metrics, as the number of random point substitutions and deletions increases, the distance from the original parent sequence also increases. However, under p-distances with sequence alignment, the distance between the parent sequence and the sequence with 25 random nucleotide substitutions is much smaller than the distance generated from deletion of those same 25 nucleotides even though, in theory, they should be identical. This disparity is due to the erroneous alignment for sequences with deletion mutations. The misalignment is so influential that the sequence with 50-nucleotide substitutions is considered closer to the parent sequence than the sequence with only 25 random deletions. In contrast, all alignment-free methods produce distances where the 25/50-nucleotide substitution is much more commensurate with its deletion counterpart.

Another difference between the alignment-based method and the alignment-free methods is the effect of larger contiguous mutations on the estimated distances. For alignment-free methods, the distances from the parent sequence generated by 100-nucleotide indel mutations are smaller than those generated by the 25- and 50-nucleotide random substitutions and deletions; the opposite is the case under the alignment-based approaches. As in the case of the CDS of $\beta$-globin gene for the chimp, these large contiguous mutations were single mutation events and so the smaller distance reflects the history of mutation events.

Since the differences between the k-mer distances and the volume intersection method were so subtle for the real DNA sequences examined, we take advantage of these synthetic sequences with its larger recombination mutations to more closely examine the commonalities between the two methods as well as where they diverge. For all alignment-free methods, sequences that have transposition mutations remain closest to the original parent sequence as all of the original elements of the parent sequence are still retained in the child after transposition; additionally, the two sequences with transposition mutations produce very similar distances, which can be seen in the first two rows of the heatmaps. In contrast, under the alignment-based method, the distance generated by the transposition sequence $\pm$ 50NT/50:99 $\rightarrow$ 350:399 is inadvertently considered to be more distant than the sequence $\pm$ 50NT/50:99 $\rightarrow$ 150:199; this is counter-intuitive because the two mutations are quite similar and affect the same number of nucleotides.

The differences between the shape-based methods and the k-mer distance are much more subtle, but still noteworthy. For example, under k-mer distances, the sequence that was generated by two separate deletion events ($-$100NT/50:149,$-$100NT/500:599) produces a much greater distance from the parent sequence than the sequence generated by two similar insertion events ($-$100NT/50:149,$-$100NT/500:599). However, the total number of mutation events is identical between the two sequences; the disparity in the k-mer distance is due merely to the changing distribution of the k-mers resulting from the recombination of sequences. In contrast, these two sequences produce comparable distances to each other for the shape-based methods; their distance from the original parent sequence is also very small, better reflecting the low number of mutation events.

Another noticeable difference between the shape-based distances and k-mer distances is in, for instance, all the pairwise distances between the sequence with 50 random point substitutions ($\Delta$50NT) and all other sequences. As stated earlier, the number of modifications between $\Delta$50NT and any sequences with large contiguous deletions, insertions, and transpositions is always between 51 and 52; the same is true for the $-$50NT sequence. The shape-based approach accurately identifies the two sequences, $\Delta$50NT and $-$50NT, as the most consistently different in all pairwise comparisons with other sequences (excluding their progenitor sequences, $\Delta$25NT and $-$25NT), as shown by the red rows and columns corresponding to $\Delta$50NT and $-$50NT. But, additionally, the range of observed distances between them and the sequences with large recombinant mutations is also small, reflecting the narrow range for the number of separating mutation events. The shape-based distances for $\Delta$50NT (and $-$50NT) to all other sequences with recombinant mutations are consistently high in magnitude when compared against the other sequences. In contrast, for the k-mer method with $k=10$, the distance for $\Delta$50NT (and $-$50NT) is largest when it is evaluated against the sequence with two large deletion events, $-$100NT/50:149,$-$100NT/500:599, but does not stay similarly high for the other recombinant mutation sequences. Overall, we can see that the color patterns in the shape-based methods closely follow the heatmap for the number of pairwise-observable differences of Figure 9, which is not seen in the k-mer method.

The resulting phylogenetic trees for these synthetic sequences reflect not only the distances shown in the heatmap, but the effectiveness of the neighbor-joining agglomeration method. For all trees, the pair of sequences with overlapping substitution mutations are clustered together and those with overlapping deletion mutations are always clustered together. The sequences with shred insertion mutations are also quite close in the tree. As in the case of the $\beta$-globin gene, the agglomeration method used on the synthetic DNA does affect the sequence alignment tree; under UPGMA, the sequence alignment approach fails to place the 50-nucleotide deletion sequence near the 25-nucleotide deletion sequence, its progenitor/antecedent sequence. The other methods are less affected by agglomeration method.

In Figure 12, we show heatmaps for the distances of the alignment-free approaches under alternative parameter values. Figure 12 shows heatmaps for: 1) k-mers of length 3 and 10 and with distances set to either Euclidean or the default Edgar (2004) distance,⁶⁶6The k-mers with $k>10$ were computationally intractable. 2) volume intersection distances with bandwidths equal to $0.003\left(\frac{1}{16}\right),0.003\left(\frac{1}{4}\right),2\times 0.003$, and $4\times 0.003)$; and 3) shape signature distances for varying histogram bin counts as a fraction, $\frac{1}{60}$th and $\frac{1}{30}$th, or a multiple, 1 and 2, of the average sequence length (equivalent to changing the number of bins in each feature from the original 46 to 12, 23, 685, and 1370). Much more granular parameterizations were explored; the subset presented here are sufficient to illustrate the observed changes in distances as a function of parameter settings.

Most notably, the distances and phylogenies derived from k-mers is quite different under the Edgar method. As expected, there is a vast deterioration in the quality of the distances when reducing k-mer length down to 3. And there is an improvement in the quality of the distances when increasing the k-mer length to 10 as the pairwise distances between recombined sequences and $\Delta$50NT increase and they become relatively more equal; the same applies to the pairwise distances between recombined sequences and $-$50NT. But the noted elevation in distance between the 50-nucleotide substitution and the 100-nucleotide insertion into positions 50:149 and 500:599 remains unchanged. Additionally, there is a limit to how much the distances can be improved under k-mers given that longer k-mer lengths become computationally difficult as the number of unique k-mers is $4^{k}$.

Both the volume intersection method and the shape signature method rely only on one parameter. Both decreasing the size of the bandwidth as well as increasing the number of histogram bins improve the distance metrics by bringing their relative distances more in line with our expectations shown in the ground-truth plots. However, in contrast to the k-mers approach, it is not more computationally demanding to change the bandwidth of the volume intersection method and the computational demands of the shape signature is driven primarily by the sequence length rather than bins (though practical issues like memory usage and cache performance can cause slowdowns for large bin counts). The changes in the heatmaps as a function of parameter values are subtle and gradual, but when we increase the granularity of the shape parameters the heatmaps begin to better reflect the true number of observable mutation events separating the sequences.

Finally, not only are the shape-based methods cheap to reparameterize, they produce relative distances that are quite consistent across different parameterizations. We show extremely wide parameter settings in order to demonstrate the effects of parameter selection; but we can see that the heatmaps for our presented results (Figure 11), with a bandwidth of 0.003 for the volume intersection method and bin count equal to 1/15th the average sequence length (46 bins), are quite similar to heatmaps of parameters with the much higher resolution in Figure 12. When testing different parameter settings, we observed that the heatmaps were relatively stable unless very extreme values were chosen. This also translates into much more stable phylogenies that are largely invariant, especially at leaf node clusterings, to parameter specifications. Both of the baseline methods using sequence alignment and k-mers required particular agglomeration methods in order to perform well across all 3 DNA sequences examined here; this was not the case for the shape-based methods.

The overall summary of results from mutation aggregation simulations for both deletion and substitution mutations are plotted in Figure 13. The figures show that, on average, the distance calculated by shape-based methods for a deletion mutation is roughly mirrored by the distance created by a substitution mutation. In contrast, the method based on sequence alignment heavily relies on correct sequence alignment; since a deletion mutation creates greater difficulty in alignment, the effect is that a deletion mutation often appears to induce a greater distance than an equivalent substitution mutation. A relationship between k-mer distance due to substitution versus deletion mutations is similar to shape-based methods for low k-mer lengths. The bulk of the linear relationship observed between substitution and deletion mutations occur for the lowest values of $k$. As the value of $k$ increases, the shorter the linear portion of the relationship becomes. Increasing the $k$ comes at the cost of the linear relationship as deletion mutations become increasingly more penalized for higher values of $k$. However, keeping $k$ low to preserve the linearity leads to poorer characterization of the DNA sequences. For k-mers, as mutations accumulate, the distance generated by a deletion mutation deviates away from that generated by a substitution mutation regardless of the length of k; it merely takes more mutations for lower $k$ values to do so. The k-mer distances calculated under the Edgar method (not plotted) does produce a linear relationship between substitution and deletion mutation distances, but, as noted previously, the Edgar distances also perform poorly for phylogenetic analyses.

A more detailed breakdown of the shape-based methods are provided in Figure 14 which shows how the shape parameters (bandwidth for the volume intersection method and the number of bins for the shape signature method) affect the induced distances. Results from the volume intersection method are relatively straightforward; the distance created by a deletion mutation is roughly mirrored by the distance created by a substitution mutation, especially when bandwidths are small. The choice of bandwidth affects how quickly mutations approach a distance of 1 between a parent and child sequence; smaller bandwidths, which reduce the size of the Gaussian spheres and reduce the likelihood of intersections, approach a distance of 1 more quickly than larger bandwidths. For the shape signatures, a similar dynamic is observed as the number of histogram bins increase. Although the CGR for each sequence is unique, the shape signature is not necessarily unique; however, increasing the number of bins used in the histogram reduces the chance that two different sequences have the same shape signature. Both smaller bandwidths and greater number of bins correspond to capturing more detail in the data and so smaller bandwidths and greater bin counts both correspond to greater congruence between the distance created by a deletion mutation and that created by a substitution mutation. However, the shape signature method can potentially suffer when the bin counts are increased so much that very similar CGR values fall into separate bins and cannot contribute to the estimated similarity. The volume intersection method with its Gaussian spheres is better at allowing for small misalignments due to the lack of a hard edge cutoff.

The relationship between distances caused by deletions and distances caused by substitutions only start to deviate for the shape signature method when: 1) a large proportion of the original sequence has been mutated and 2) the shape parameter chosen is not very fine-grained (i.e., when bin counts are low). Additionally, larger bandwidths and fewer bin counts also cause greater variability in estimated distances due to deletion and substitution.