Mapping the Space of Genomic Signatures
Lila Kari, Kathleen A. Hill, Abu S. Sayem, Rallis Karamichalis, Nathaniel Bryans, Katelyn Davis, Nikesh S. Dattani
1 Department of Computer Science, University of Western Ontario, London, ON, N6A 5B7, Canada
2 Department of Biology, University of Western Ontario, London, ON, N6A 5B7, Canada
3 Physical and Theoretical Chemistry Laboratory, Department of Chemistry, Oxford University, Oxford, OX1 3QZ, United Kingdom
We propose a computational method to measure and visualize interrelationships among any number of DNA sequences allowing, for example, the examination of hundreds or thousands of complete mitochondrial genomes. An “image distance” is computed for each pair of graphical representations of DNA sequences, and the distances are visualized as a Molecular Distance Map: Each point on the map represents a DNA sequence, and the spatial proximity between any two points reflects the degree of structural similarity between the corresponding sequences. The graphical representation of DNA sequences utilized, Chaos Game Representation (CGR), is genome- and species-specific and can thus act as a genomic signature. Consequently, Molecular Distance Maps could inform species identification, taxonomic classifications and, to a certain extent, evolutionary history. The image distance employed, Structural Dissimilarity Index (DSSIM), implicitly compares the occurrences of oligomers of length up to (herein ) in DNA sequences. We computed DSSIM distances for more than 5 million pairs of complete mitochondrial genomes, and used Multi-Dimensional Scaling (MDS) to obtain Molecular Distance Maps that visually display the sequence relatedness in various subsets, at different taxonomic levels.
This general-purpose method does not require DNA sequence homology and can thus be used to compare similar or vastly different DNA sequences, genomic or computer-generated, of the same or different lengths. We illustrate potential uses of this approach by applying it to several taxonomic subsets: phylum Vertebrata, (super)kingdom Protista, classes Amphibia-Insecta-Mammalia, class Amphibia, and order Primates. This analysis of an extensive dataset confirms that the oligomer composition of full mtDNA sequences can be a source of taxonomic information. This method also correctly finds the mtDNA sequences most closely related to that of the anatomically modern human (the Neanderthal, the Denisovan, and the chimp), and that the sequence most different from it belongs to a cucumber.
Every year biologists discover and classify thousands of new species, with an average of 18,000 new species named each year since 1940 , . Other findings, , suggest that as many as 86% of existing species on Earth and 91% of species in the oceans have not yet been classified and catalogued. In the absence of a universal quantitative method to identify species’ relationships, information for species classification has to be gleaned and combined from several sources, morphological, sequence-alignment-based phylogenetic anaylsis, and non-alignment-based molecular information.
We propose a computational process that outputs, for any given dataset of DNA sequences, a concurrent display of the structural similarities among all sequences in the dataset. This is obtained by first computing an “image distance” for each pair of graphical representations of DNA sequences, and then visualizing the resulting interrelationships in a two-dimensional plane. The result of applying this method to a collection of DNA sequences is an easily interpretable Molecular Distance Map wherein sequences are represented by points in a common Euclidean plane, and the spatial distance between any two points reflects the differences in their subsequence composition.
The graphical representation we use is Chaos Game Representation (CGR) of DNA sequences, [4, 5], that simultaneously displays all subsequence frequencies of a given DNA sequence as a visual pattern. CGR has a remarkable ability to differentiate between genetic sequences belonging to different species, and has thus been proposed as a genomic signature. Due to this characteristic, a Molecular Distance Map of a collection of genetic sequences may allow inferrences of relationships between the corresponding species.
Concretely, to compute and visually display relationships in a given set of DNA sequences, we propose the following computational process:
(i) Chaos Game Representation (CGR), to graphically represent all subsequences of a DNA sequence , , as pixels of one image, denoted by ;
(ii) Structural Dissimilarity Index (DSSIM), an “image-distance” measure, to compute the pairwise distances , , for each pair of CGR images , and to produce a distance matrix;
(iii) Multi-Dimensional Scaling (MDS), an information visualization technique that takes as input the distance matrix and outputs a Molecular Distance Map in 2D, wherein each plotted point with coordinates represents the DNA sequence whose CGR image is . The position of the point in the map, relative to all the other points , reflects the distances between the DNA sequence and the other DNA sequences in the dataset.
We apply this method to analyze and visualize several different taxonomic subsets of a dataset of 3,176 complete mtDNA sequences: phylum Vertebrata, (super)kingdom Protista, classes Amphibia-Insecta-Mammalia, class Amphibia only, and order Primates. We illustrate the usability of this approach by discussing, e.g., the placement of the genus Polypterus within phylum Vertebrata, of the unclassified organism Haemoproteus sp. jb1.JA27 (#1466) within the (super)kingdom Protista, and the placement of the family Tarsiidae within the order Primates. We also provide an interactive web tool, MoD-Map, that allows an in-depth exploration of all Molecular Distance Maps in the paper, complete with zoom-in features, search options, and easily accessible additional information for each species-point sequence.
Overall, this method groups mtDNA sequences in correct taxonomic groups, from the kingdom level down to the order and family level. These results are of interest both because of the size of the dataset analyzed and because this information was extracted from DNA sequences that normally would not be considered homologous. Our analysis confirms that sequence composition (presence or absence of oligomers) contains taxonomic information that could be relevant to species identification, taxonomic classification, and identification of large evolutionary lineages. Last but not least, the appeal of this method lies in its simplicity, robustness, and generality, whereby the exact same measuring tape is able to automatically yield meaningful measurements between non-specific DNA sequences of species as distant as those of the anatomically modern human and a cucumber, and as close as those of the anatomically modern human and the Neanderthal.
A CGR [4, 5] associates an image to each DNA sequence as follows. Begin with a unit square with corners labelled A, C, G, and T, clockwise starting from the bottom-left corner. The first point of any CGR plot is the center of the square. To plot the CGR corresponding to a given DNA sequence, start reading the letters of the sequence from left to right, one by one. The point corresponding to the first letter is the point plotted in the middle of the segment determined by the center of the square and the corner labelled by the first letter. For example, if the center of the square is labelled “O” and the first letter of the sequence is “A”, then the point of the plot corresponding to the first “A” is the point situated halfway between O and the corner A. Subsequent letters are plotted iteratively as the middle point between the previously-drawn-point and the corner labelled by the letter currently being read.
CGR images of genetic DNA sequences originating from various species show rich fractal patterns containing various motifs such as squares, parallel lines, rectangles, triangles and diagonal crosses, see, e.g., Figure 1. CGRs of genomic DNA sequences have been shown to be genome and species specific, [4, 5, 6, 7, 8, 9, 10]. Thus, sequences chosen from each genome as a basis for computing “distances” between genomes do not need to have any relation with one another from the point of view of their position or information content. In addition, this graphical representation facilitates easy visual recognition of global string-usage characteristics: Prominent diagonals indicate purine or pyrimidine runs, sparseness in the upper half indicates low G+C content, etc., see, e.g., .
If the generated CGR image has a resolution of pixels, then every pixel represents a distinct DNA subsequence of length : A pixel is black if the subsequence it represents occurs in the DNA sequence, otherwise it is white. In this paper, for the CGR images of all 3,176 complete mtDNA sequences in our dataset, we used the value , that is, occurrences of subsequences of lengths up to 9 were being taken into consideration. In general, a length of the DNA sequence of 4,000 bp is necessary to obtain a sharply defined CGR, but in many cases 2,000 bp give a reasonably good approximation, . In our case, we used the full length of all analyzed mtDNA sequences, which ranged from 288 bp to 1,555,935 bp, with an average of 28,000 bp.
Other visualizations of genetic data include the 2D rectangular walk  and methods similar to it in ,, vector walk , cell , vertical vector , Huffman coding , and colorsquare  methods. Three-dimensional representations of DNA sequences include the tetrahedron , 3D-vector , and trinucleotide curve  methods. Among these visualization methods, CGR images arguably provide the most immediately comprehensible “signature” of a DNA sequence and a desirable genome-specificity, [4, 9]. In addition, the images produced using CGR are easy to compare, visually and computationally. Coloured versions of CGR, wherein the colour of a point corresponds to the frequency of the corresponding oligomer in the given DNA sequence (from red for high frequency, to blue for no occurrences) have also been proposed [22, 23].
Note that other alignment-free methods have been used for phylogenetic analysis of DNA strings, such as computing the Euclidean distance between frequencies of -mers () for the analysis of 125 GenBank DNA sequences from 20 bird species and the American alligator, . Another study, , analyzed 459 dsDNA bacteriophage genomes and compared them with their host genomes to infer host-phage relationships, by computing Euclidean distances between frequencies of -mers for . In , 75 complete HIV genome sequences were compared using the Euclidean distance between frequencies of 6-mers (), in order to group them in subtypes. In , 27 microbial genomes were analyzed to find implications of 4-mer frequencies () on their evolutionary relationships. In , 20 mammalian complete mtDNA sequences were analyzed using the “similarity metric”. Our method uses a larger dataset (3,176 complete mtDNA sequences), an “image distance” measure that was designed to capture structural similarities between images, as well as a value of .
Structural Similarity (SSIM) index is an image similarity index used in the context of image processing and computer vision to compare two images from the point of view of their structural similarities . SSIM combines three parameters - luminance distortion, contrast distortion, and linear correlation - and was designed to perform similarly to the human visual system, which is highly adapted to extract structural information. Originally, SSIM was defined as a similarity measure whose theoretical range between two images and is where a high value amounts to close relatedness. We use a related DSSIM distance , with the distance being 0 between two identical images, 1 between e.g. a black image and a white image, and 2 if the two images are negatively correlated; that is, if and only if every pixel of image has the inverted value of the corresponding pixel in image while both images have the same luminance (brightness). For our particular dataset of genetic CGR images, almost all (over 5 million) distances are between 0 and 1, with only half a dozen exceptions of distances between 1 and 1.0033.
MDS has been used for the visualization of data relatedness based on distance matrices in various fields such as cognitive science, information science, psychometrics, marketing, ecology, social science, and other areas of study . MDS takes as input a distance matrix containing the pairwise distances between given items and outputs a two-dimensional map wherein each item is represented by a point, and the spatial distances between points reflect the distances between the corresponding items in the distance matrix. Notable examples of molecular biology studies that used MDS are  (where it was used for the analysis of geographic genetic distributions of some natural populations),  (where it was used to provide a graphical summary of the distances among CO1 genes from various species), and  (where it was used to analyze and visualize relationships within collections of phylogenetic trees).
Classical MDS, which we use in this paper, receives as input an distance matrix of the pairwise distances between any two items in the set. The output of classical MDS consists of points in a -dimensional space whose pairwise spatial (Euclidean) distances are a linear function of the distances between the corresponding items in the input distance matrix. More precisely, MDS will return points such that for all where is the spatial distance between the points and , and is a function linear in . Here, can be at most and the points are recovered from the eigenvalues and eigenvectors of the input distance matrix. If we choose (respectively ), the result of classic MDS is an approximation of the original -dimensional space as a two- (respectively three-) dimensional map.
In this paper all Molecular Distance Maps consist of coloured points, wherein each point represents an mtDNA sequence from the dataset. Each mtDNA sequence is assigned a unique numerical identifier retained in all analyses, e.g., #1321 is the identifier for the Homo sapiens sapiens mitochondrial genome. The colour assigned to a sequence-point may however vary from map to map, and it depends on the taxon assigned to the point in a particular Molecular Distance Map. For example, in Figure 2 all mammalian mtDNA sequence-points are coloured red, while in Figure 6 the red points represent mtDNA sequences from the primate suborder Haplorhini and the green points represent mtDNA sequences from the primate suborder Strepshirrini. For consistency, all maps are scaled so that the - and the -coordinates always span the interval . The formula used for scaling is , , where and are the minimum and maximum of the -coordinates of all the points in the original map, and similarly for and .
Each Molecular Distance Map has some error, that is, the spatial distances are not exactly the same as . When using the same dataset, the error is in general lower for an MDS map in a higher-dimensional space. The Stress-1 (Kruskal stress, ), is defined in our case as
where the summations extend over all the sequences considered for a given map, and is a linear function whose parameters are determined by linear regression for each dataset and corresponding Molecular Distance Map. A benchmark that is often used to assess MDS results is that Stress-1 should be in the range , see .
The dataset consists of the entire collection of complete mitochondrial DNA sequences from NCBI as of 12 July, 2012. This dataset consists of 3,176 complete mtDNA sequences, namely 79 protists, 111 fungi, 283 plants, and 2,703 animals. This collection of mitochondrial genomes has a great breadth of species across taxonomic categories and great depth of species coverage in certain taxonomic categories. For example, we compare sequences at every rank of taxonomy, with some pairs being different at as high as the (super)kingdom level, and some pairs of sequences being from the exact same species, as in the case of Silene conica for which our dataset contains the sequences of 140 different mitochondrial chromosomes . The prokaryotic origins and evolutionary history of mitochondrial genomes have long been extensively studied, which will allow comparison of our results with known relatedness of species. Lastly, this genome dataset permits testing of both recent and deep rooted species relationships, providing fine resolution of species differences.
The creation of the datasets, acquisition of data from NCBI’s GenBank, generation of the CGR images, calculation of the distance matrix, and calculation of the Molecular Distance Maps using MDS, were all done (and can be tested with) the free open-source MATLAB program OpenMPM . This program makes use of an open source MATLAB program for SSIM, , and MATLAB’s built-in MDS function. The interactive web tool MoD-Map (Molecular Distance Map), , allows an in-depth exploration and navigation of the Molecular Distance Maps in this paper.
Results and Discussion
The Molecular Distance Maps we analyzed, of several different taxonomic subsets (phylum Vertebrata, (super)kingdom Protista, classes Amphibia-Insecta-Mammalia, class Amphibia only, and order Primates), confirm that the presence or absence of oligomers in mtDNA sequences may contain information that is relevant to taxonomic classifications. These results are of interest both because of the large dataset considered and because this information has been extracted from DNA sequences that, by normal criteria, would be considered nonhomologous. The main contributions of the paper are the following:
The use of an “image distance” (designed to detect structural similarities between images) to compare the graphic signatures of two DNA sequences. For any given , this distance simultaneously compares the occurrences of all subsequences of length up to of the two sequences. In all computations of this paper we use . This image distance (with parameter set to ) is highly sensitive and succeeds to successfully group hundreds of CGRs that are visually similar, such as the ones in Figure 1(left) and Figure 1(right), into correct taxonomic categories.
The use of an information visualization technique to display the results as easily interpretable Molecular Distance Maps, wherein the spatial position of each sequence-point in relation to all other sequence-points is quantitatively significant. This is augmented by an interactive web tool which allows an in-depth exploration of the Molecular Distance Maps in this paper, with features such as zoom-in, search by scientific name or NCBI accession number, and quick access to complete information for each of the full mtDNA sequences in the map.
A method that is general-purpose, simple, computationally efficient and scalable. Since the compared sequences need not be homologous or of the same length, this method can be used to provide comparisons among any number of completely different DNA sequences: within the genome of an individual, across genomes within a single species, between genomes within a taxonomic category, and across taxa.
The use of a large dataset of 3,176 complete mitochondrial DNA sequences.
An illustration of potential uses of this approach by the discussion of several case studies such as the placement of the genus Polypterus within phylum Vertebrata, of the unclassified organism Haemoproteus sp. jb1.JA27 (#1466) within the (super)kingdom Protista, and the placement of the family Tarsiidae within the order Primates.
This method could complement information obtained by using DNA barcodes  and Klee diagrams , since it is applicable to cases where barcodes may have limited effectiveness: plants and fungi for which different barcoding regions have to be used , , ; protists where multiple loci are generally needed to distinguish between species ; prokaryotes ; and artificial, computer-generated, DNA sequences. This method may also complement other taxonomic analyses by bringing in additional information gleaned from comparisons of non-homologous and non-coding sequences.
An example of the CGR/DSSIM/MDS approach is the Molecular Distance Map in Figure 2 which depicts the complete mitochondrial DNA sequences of all 1,791 jawed vertebrates in our dataset. (In the legends of Figures 2-6, the number of represented mtDNA sequences in each category is listed in paranthesis after the category name.) Note that the position of each point in a map is determined by all the distances between the sequence it represents and the other sequences in the dataset. In the case of Figure 2, the position of each sequence-point is determined by the 1,790 numerical distances between its sequence and all the other mtDNA vertebrate sequences in that dataset.
Observe that all five different subphyla of jawed vertebrates are separated in non-overlapping clusters, with very few exceptions. Examples of fish species bordering or slightly mixed with the amphibian cluster include Polypterus ornatipinnis (#3125, ornate bichir), Polypterus senegalus (#2868, Senegal bichir), both with primitive pairs of lungs; Erpetoichthys calabaricus (#2745, reedfish) who can breathe atmospheric air using a pair of lungs; and Porichtys myriaster (#2483, specklefish midshipman) a toadfish of the order Batrachoidiformes. It is noteworthy that the question of whether species of the Polypterus genus are fish or amphibians has been discussed extensively for hundreds of years . Interestingly, all four represented lungfish (a.k.a. salamanderfish), are also bordering the amphibian cluster: Protopterus aethiopicus (#873, marbled lungfish), Lepidosiren paradoxa (#2910, South American lungfish), Neoceratodus forsteri (#2957, Australian lungfish), Protopterus doloi (#3119, spotted African lungfish). Note that, in answer to the hypothesis in  regarding the diversity of signatures across vertebrates, in Figure 2, the avian mtDNA signatures cluster neither with the mammals nor with the reptiles, and form a completely separate cluster of their own (albeit closer to reptiles than to mammals).
We applied our method to visualize the relationships among all represented species from the (super) kingdom Protista whose taxon, as defined in the legend of Figure 3, had more than one representative. As expected, the maximum distance between pairs of sequences in this map was higher than the maximum distances for the other maps in this paper, all at lower taxonomic levels.
The most obvious outlier in Figure 3 is Haemoproteus sp. jb1.JA27 (#1466), sequenced in  (see also ), and listed as an unclassified organism in the NCBI taxonomy. Note first that this species-point belongs to the same kingdom (Chromalveolata), superphylum (Alveolata), phylum (Apicomplexa), and class (Aconoidasida), as the other two species-points that appear grouped with it, Babesia bovis T2Bo (#1935), and Theileria parva (#3173). This indicates that its position is not fully anomalous. Moreover, as indicated by the high value of Stress-1 for this figure, an inspection of DSSIM distances shows that this species-point may not be a true outlier, and its position may not be as striking in a higher dimensional version of the Molecular Distance Map. Overall, this map shows that our method allows an exploration of diversity at the level of super kingdom, obtains good clustering of known subtaxonomic groups, while at the same time indicating a lack of genome sequence information and paucity of representation that complicates analyses for this fascinating taxonomic group.
We then applied our method to visualize the relationships between all available complete mtDNA sequences from three classes, Amphibia, Insecta and Mammalia (Figure 4), as well as observe relationships within class Amphibia and three of its orders (Figure 5). Note that a feature of MDS is that the points are not unique. Indeed, one can translate or rotate a map without affecting the pairwise spatial distances . In addition, the obtained points in an MDS map may change coordinates when more data items are added to or removed from the dataset. This is because the output of the MDS aims to preserve only the pairwise spatial distances between points, and this can be achieved even when some of the points change their coordinates. In particular, the -coordinates of a point representing an amphibian species in the amphibians-insects-mammals map (Figure 4) will not necessarily be the same as the -coordinates of the same point when only amphibians are mapped (Figure 5).
In general, Molecular Distance Maps are in good agreement with classical phylogenetic trees at all scales of taxonomic comparisons, see Figure 5 with , and Figure 6 with. In addition, our approach may be able to weigh in on conflicts between taxonomic classifications based on morphological traits and those based on more recent molecular data, as in the case of tarsiers, as seen below.
Zooming in, we observed the relationships within an order, Primates, with its suborders (Figure 6). Notably, two extinct species of the genus Homo are represented: Homo sapiens neanderthalensis and Homo sapiens ssp. Denisova. Primates can be classified into two groups, Haplorhini (dry-nosed primates comprising anthropoids and tarsiers) and Strepsirrhini (wet-nosed primates including lemurs and lorises). The map shows a clear separation of these suborders, with the top-left arm of the map in Figure 6, comprising the Strepsirrhini. However, there are two Haplorhini placed in the Strepsirrhini cluster, namely Tarsius bancanus (#2978, Horsfield’s tarsier) and Tarsius syrichta (#1381, Philippine tarsier). The phylogenetic placement of tarsiers within the order Primates has been controversial for over a century, . According to , mitochondrial DNA evidence places tarsiiformes as a sister group to Strepsirrhini, while in contrast,  places tarsiers within Happlorhini. In Figure 6 the tarsiers are located within the Strepsrrhini cluster, thus agreeing with . This may be partly because both this study and  used mitochondrial DNA, whose signature may be different from that of chromosomal DNA as seen in Figure 1(left) and Figure 1(center).
The DSSIM distances computed between all pairs of complete mtDNA sequences varied in range. The minimum distance was 0, between two pairs of identical mtDNA sequences. The first pair comprised the mtDNA of Rhinomugil nasutus (#98, shark mullet, length 16,974 bp) and Moolgarda cunnesius (#103, longarm mullet, length 16,974 bp). A base-to-base sequence comparison between these sequences (#98, NC_017897.1; #103, NC_017902.1) showed that the sequences were indeed identical. However, after completion of this work, the sequence for species #103 was updated to a new version (NC_017902.2), on 7 March, 2013, and is now different from the sequence for species #98 (NC_017897.1). The second pair comprises the mtDNA sequences #1033 and #1034 (length 16,623 bp), generated by crossing female Megalobrama amblycephala with male Xenocypris davidi leading to the creation of both diploid (#1033) and triploid (#1034) nuclear genomes, , but identical mitochondrial genomes.
The maximum distance was found to be between Pseudendoclonium akinetum (# 2656, a green alga, length 95,880) and Candida subhashii (#954, a yeast, length 29,795). Interestingly, the pair with the maximum distance featured neither the longest mitochondrial sequence, with the darkest CGR (Cucumis sativus, #533, cucumber, length 1,555,935 bp), nor the shortest mitochondrial sequence, with the lightest CGR (Silene conica, #440, sand catchfly, a plant, length 288 bp).
An inspection of the distances between Homo sapiens sapiens and all the other primate mitochondrial genomes in the dataset showed that the minimum distance to Homo sapiens sapiens was , the distance to Homo sapiens neanderthalensis (#1720, Neanderthal), with the second smallest distance to it being , the distance to Homo sapiens ssp. Denisova (#1052, Denisovan). The third smallest distance was to Pan troglodytes (#3084, chimp). Figure 7 shows the graph of the distances between the Homo sapiens sapiens mtDNA and each of the primate mitochondrial genomes. With no exceptions, this graph is in full agreement with established phylogenetic trees, . The largest distance between the Homo sapiens sapiens mtDNA and another mtDNA sequence in the dataset was 0.9957, the distance between Homo sapiens sapiens and Cucumis sativus (#533, cucumber, length 1,555,935 bp).
In addition to comparing real DNA sequences, this method can compare real DNA sequences to computer-generated sequences. As an example, we compared the mtDNA genome of Homo sapiens sapiens with one hundred artificial, computer-generated, DNA sequences of the same length and the same trinucleotide frequencies as the original. The average distance between these artificial sequences and the original human mitochondrial DNA is 0.8991. This indicates that all “human” artificial DNA sequences are more distant from the Homo sapiens sapiens mitochondrial genome than Drosophila melanogaster (#3120, fruit fly) mtDNA, with . This further implies that trinucleotide frequencies may not contain sufficient information to classify a genetic sequence, suggesting that Goldman’s claim  that “CGR gives no futher insight into the structure of the DNA sequence than is given by the dinucleotide and trinucleotide frequencies” may not hold in general.
The Stress-1 values for all but one of the Molecular Distance Maps in this paper were in the “acceptable” range . The exception was Figure 3 with Stress-1 equal to 0.26. Note that Stress-1 generally decreases with an increase in dimensionality, from to . Note also that, as suggested in , the Stress-1 guidelines are not absolute: It is not always the case that only MDS representations with Stress-1 under are acceptable, nor that all MDS representations with Stress-1 under are good.
In all the calculations in this paper, we used the full mitochondrial sequences. However, since the length of a sequence can influence the brightness of its CGR and thus its Molecular Distance Map coordinates, further analysis is needed to elucidate the effect of sequence length on the positions of sequence-points in a Molecular Distance Map. The choice of length of DNA sequences used may ultimately depend on the particular dataset and particular application.
We now discuss some limitations of the proposed methods. Firstly, DSSIM is very effective at picking up subtle differences between images. For example, all vertebrate CGRs present the triangular fractal structure seen in the human mtDNA, and are visually very similar, as seen in Figure 1(left) and Figure 1(right). In spite of this, DSSIM is able to detect a range of differences that is sufficient for a good positioning of all 1,791 mtDNA sequences relative to each other. This being said, DSSIM may give too much weight to subtle differences, so that small and big differences in images produce distances that are numerically very close. This may be a useful feature for the analysis of datasets of closely related sequences. For large-scale taxonomic comparisons however, refinements of DSSIM or the use of other distances needs to be explored, that would space further apart the values of distances arising from small differences versus those arising from big-pattern differences between images.
Secondly, MDS always has some errors, in the sense that the spatial distance between two points does not always reflect the original distance in the distance matrix. For fine analyses, the placement of a sequence-point in a map has to be confirmed by checking the original distance matrix. Possible solutions include increasing the dimensionality of the maps to three-dimensional maps, which are still easily interpretable visually and have been shown in some cases to separate clusters which seemed incorrectly intermeshed in the two-dimensional version of the map. Other possibilities include a colour-scheme that would colour points with low stress-per-point differently from the ones with high stress-per-point, and thus alert the user to the regions where discrepancies between the spatial distance and the original distance exist.
Thirdly, we note that the use of the particular distance measure (DSSIM) or particular scaling technique (classical MDS) does not mean that these are the optimal choices in all cases.
Our analysis suggests that the oligomer composition of mitochondrial DNA sequences can be a source of taxonomic information. These results are of interest both because of the large dataset considered (see, e.g., the correct grouping in taxonomic categories of 1,791 mitochondrial genomes in Figure 2), and because this information is extracted from DNA sequences that, by normal criteria, would not be considered homologous.
Potential applications of Molecular Distance Maps - when used on a dataset of genomic sequences, whether coding or non-coding, homologous or not homologous, of the same length or vastly different lengths – include identification of large evolutionary lineages, taxonomic classifications, species identification, as well as possible quantitative definitions of the notion of species and other taxa.
Possible extensions include generalizations of MDS, such as 3-dimensional MDS, for improved visualization, and the use of increased oligomer length (higher values of ) for comparisons of longer subsequences in case of whole chromosome or whole genome analyses. We note also that this method can be applied to analyzing sequences over other alphabets. For example binary sequences could be imaged using a square with vertices labelled 00, 01, 10, 11, and then DSSIM and MDS could be employed to compare and map them.
We thank Ronghai Tu for an early version of our MATLAB code to generate CGR images, Tao Tao for assistance with NCBI’s GenBank, Steffen Kopecki for generating artificial sequences and discussions. We also thank Andre Lachance, Jeremy McNeill, and Greg Thorn for resources and discussions on taxonomy. We thank the Oxford University Mathematical Institute for the use of their Windows compute server Pootle/WTS.
- 1. IISE (2012) Retro SOS 2000 - 2009: A decade of species discovery in review. International Institute for Species Exploration : Retrieved Sept. 26, 2014 from http://www.esf.edu/species/SOS.htm.
- 2. Frazer J (2014) Top 10 New Species of 2014. National Geographic News, May 22, 2014 : http://news.nationalgeographic.com/news/2014/05/140522-top-ten-new-species-2014-biodiversity/.
- 3. Mora C, Tittensor D, Adl S, Simpson A, Worm B (2011) How many species are there on earth and in the ocean? PLoS Biology 9: 1-8.
- 4. Jeffrey H (1990) Chaos Game Representation of gene structure. Nucleic Acids Research 18: 2163–2170.
- 5. Jeffrey H (1992) Chaos game visualization of sequences. Comput Graphics 16: 25-33.
- 6. Hill K, Schisler N, Singh S (1992) Chaos Game Representation of coding regions of human globin genes and alcohol dehydrogenase genes of phylogenetically divergent species. J Mol Evol 35: 261-9.
- 7. Hill K, Singh S (1997) Evolution of species-type specificity in the global DNA sequence organization of mitochondrial genomes. Genome 40: 342-356.
- 8. Deschavanne P, Giron A, Vilain J, Fagot G, Fertil B (1999) Genomic signature: characterization and classification of species assessed by Chaos Game Representation of sequences. Molecular Biology and Evolution 16: 1391–1399.
- 9. Deschavanne P, Giron A, Vilain J, Dufraigne C, Fertil B (2000) Genomic signature is preserved in short DNA fragments. In: IEEE Intl. Symposium on Bio-Informatics and Biomedical Engineering. pp. 161–167.
- 10. Wang Y, Hill K, Singh S, Kari L (2005) The spectrum of genomic signatures: From dinucleotides to Chaos Game Representation. Gene 346: 173–185.
- 11. Gates M (1986) A simple way to look at DNA. J Theor Biology 119: 319–328.
- 12. Nandy A (1994) A new graphical representation and analysis of DNA sequence structure: Methodology and application to globin genes. Current Science 66: 309 - 314.
- 13. Leong P, Morgenthaler S (1995) Random walk and gap plots of DNA sequences. Computer applications in the biosciences : CABIOS 11: 503-507.
- 14. Liao B (2005) A 2D graphical representation of DNA sequence. Chemical Physics Letters 401: 196–199.
- 15. Yao Y, Wang T (2004) A class of new 2D graphical representation of DNA sequences and their application. Chemical Physics Letters 398: 318–323.
- 16. Yu C, Liang Q, Yin C, He R, Yau S (2010) A novel construction of genome space with biological geometry. DNA Research 17: 155-168.
- 17. Qi Z, Li L, Qi X (2011) Using Huffman coding method to visualize and analyze DNA sequences. Journal of Computational Chemistry 32: 3233-3240.
- 18. Zhang Z, et al. (2012) Colorsquare: A colorful square visualization of DNA sequences. Comm in Math and in Comp Chemistry 68: 621-637.
- 19. Randic M, Vracko M, Nandy A, Basak S (2000) On 3D graphical representation of DNA primary sequences and their numerical characterization. J Chem Inf and Comp Sci 40: 1235-1244.
- 20. Yuan C, Liao B, Wang T (2003) New 3D graphical representation of DNA sequences and their numerical characterization. Chemical Physics Letters 379: 412 - 417.
- 21. Yu J, Sun X, Wang J (2009) TN curve: A novel 3D graphical representation of DNA sequence based on trinucleotides and its applications. Journal of Theoretical Biology 261: 459 - 468.
- 22. Makula M, Benuskova L (2009) Interactive visualization of oligomer frequency in DNA. Computing and Informatics 28: 695-710.
- 23. Hao B, Lee H, Zhang S (2000) Fractals related to long DNA sequences and complete genomes. Chaos, Solitons and Fractals 11: 825-836.
- 24. Edwards S, Fertil B, Girron A, Deschavanne P (2002) A genomic schism in birds revealed by phylogenetic analysis of DNA strings. Systematic Biology 51: 599-613.
- 25. Deschavanne P, DuBow M, Regeard C (2010) The use of genomic signature distance between bacteriophages and their hosts diplays evolutionary relationships and phage growth cycle determination. Virology Journal 7: 163.
- 26. Pandit A, Sinha S (2010) Using genomic signatures for HIV-1 subtyping. BMC Bioinformatics 11: S26.
- 27. Pride D, Meinersmann R, Wassenaar T, Blaser M (2003) Evolutionary implications of microbial genome tetranucleotide frequency biases. Genome Research 13: 145-158.
- 28. Li M, Chen X, Li X, Ma B, Vitany P (2004) The similarity metric. IEEE Transactions on Information Theory 50: 3250-3264.
- 29. Wang Z, Bovik A, Sheikh H, Simoncelli E (2004) Image quality assessment: From error visibility to structural similarity. IEEE Transactions on Image Processing 13: 600-612.
- 30. Borg I, Groenen P (2010) Modern Multidimensional Scaling: Theory and Applications. Springer, 2nd edition.
- 31. Lessa E (1990) Multidimensional analysis of geographic genetic structure. Systematic Zoology 39(3): 242–252.
- 32. Hebert P, Cywinska A, Ball S, Dewaard J (2003) Biological identifications through DNA barcodes. Proc Biol Sci 270: 313–321.
- 33. Hillis D, Heath T, StJohn K (2005) Analysis and visualization of tree space. Systematic Biology 54: 471-482.
- 34. Kruskal J (1964) Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika 29: 1–27.
- 35. Sloan D, et al. (2012) Rapid evolution of enormous, multichromosomal genomes in flowering plant mitochondria with exceptionally high mutation rates. PLoS Biology 10: e1001241.
- 36. Dattani N, Sayem A, Tu R, Bryans N (2013) OpenMPM. Computer Program : http://git.io/Ypa_jA.
- 37. Wang Z (2003) SSIM index. Computer Program : https://ece.uwaterloo.ca/z̃70wang/research/ssim/.
- 38. Karamichalis R (2014) MoD-Map. Web Tool : http://www.csd.uwo.ca/MoDMap/.
- 39. Sirovich L, Stoeckle M, Zhang Y (2010) Structural analysis of biodiversity. PLoS ONE 5: e9266.
- 40. Kress W, Wurdack K, Zimmer E, Weigt L, Janzen D (2005) Use of DNA barcodes to identify flowering plants. PNAS 102: 8369–8374.
- 41. Hollingsworth P, et al. (2009) A DNA barcode for land plants. PNAS 106: 12794-2797.
- 42. Schoch C, et al. (2012) Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi. PNAS 109: 6241-6246.
- 43. Hoef-Emden K (2012) Pitfalls of establishing DNA barcoding systems in protists: the Cryptophyceae as a test case. PLoS One 7: e43652.
- 44. Unwin R, Maiden M (2003) Multi-locus sequence typing: a tool for global epidemiology. Trends Microbiol 11: 479–487.
- 45. Hall B (2001) John Samuel Budgett (1872-1904): In pursuit of Polypterus. BioScience 51: 399-407.
- 46. Beadell J, Fleischer R (2005) A restriction enzyme-based assay to distinguish between avian hemosporidians. Journal of Parasitology 91: 683-685.
- 47. Valkiunas G, et al. (2010) A new Haemoproteus species (Haemosporida: Haemoproteidae) from the endemic Galapagos dove Zenaida galapagoensis, with remarks on the parasite distribution, vectors, and molecular diagnostics. Journal of Parasitology 96: 783-792.
- 48. Pyron R, Wiens J (2011) A large-scale phylogeny of amphibia including over 2800 species, and a revised classification of extant frogs, salamanders, and caecilians. Molecular Phylogenetics and Evolution 61: 543-583.
- 49. Shoshani J, et al. (1996) Primate phylogeny: morphological vs molecular results. Molecular Phylogenetics and Evolution 5: 102-154.
- 50. Jameson N, et al. (2011) Genomic data reject the hypothesis of a prosimian primate clade. Journal of Human Evolution 61: 295-305.
- 51. Chatterjee H, Ho S, Barnes I, Groves C (2009) Estimating the phylogeny and divergence times of primates using a supermatrix approach. BMC Evolutionary Biology 9: 259.
- 52. Perelman P, Johnson W, Roos C, Seuánez H, Horvath J, et al. (2011) A molecular phylogeny of living primates. PLoS Genetics 7: e1001342.
- 53. Hu J, et al. (2012) Characteristics of diploid and triploid hybris derived from female Megalobrama amblycephala Yih male Xenocypris davidi Bleeker. Aquaculture 364-365: 157-164.
- 54. Goldman N (1993) Nucleotide, dinucleotide and trinucleotide frequencies explain patterns observed in Chaos Game Representations of DNA sequences. Nucleic Acids Research 21: 2487-2491.
For a complete version, including figures, please see