Structural Analysis of Biodiversity

Large, recently-available genomic databases cover a wide range of life forms, suggesting opportunity for insights into genetic structure of biodiversity. In this study we refine our recently-described technique using indicator vectors to analyze and visualize nucleotide sequences. The indicator vector approach generates correlation matrices, dubbed Klee diagrams, which represent a novel way of assembling and viewing large genomic datasets. To explore its potential utility, here we apply the improved algorithm to a collection of almost 17000 DNA barcode sequences covering 12 widely-separated animal taxa, demonstrating that indicator vectors for classification gave correct assignment in all 11000 test cases. Indicator vector analysis revealed discontinuities corresponding to species- and higher-level taxonomic divisions, suggesting an efficient approach to classification of organisms from poorly-studied groups. As compared to standard distance metrics, indicator vectors preserve diagnostic character probabilities, enable automated classification of test sequences, and generate high-information density single-page displays. These results support application of indicator vectors for comparative analysis of large nucleotide data sets and raise prospect of gaining insight into broad-scale patterns in the genetic structure of biodiversity.


Introduction
Genetic study of biodiversity has been hampered by the relatively small number of species represented in databases. For example, the largest set of alignable sequences in GenBank (small subunit ribosomal RNA) represents fewer than 21,000 species and the second largest (cytochrome b) includes fewer than 14,000 [1]. This is modest coverage compared to the approximately 1.9 million named species of plants and animals and likely much larger numbers of protozoa, fungi, bacteria, and archaea [2]. Usually, a primary goal of comparative genetic study is assembling a Tree of Life that represents the temporal sequence of evolutionary divergences. As it is computationally difficult to construct a phylogenetic tree for more than a few thousand taxa, most analyses focus on a taxonomicallyrestricted subset and select a few exemplars from each group (e.g., [3,4]). Beyond computational challenges, potential limitations to tree representations include difficulty in representing discontinuities among species or groups of species, as all taxa are linked in a continuous structure; visualizing horizontal affinities across groups, as taxa within each group are joined in a single branch; and comparing data sets such as from ecological surveys, as branching diagrams challenge visual comparison.
Large, newly-available data sets [5] offer the possibility of studying genetic diversity on a wide scale. In an earlier paper, we described a method for creating ''indicator vectors'' representative of sets of nucleotide sequences [6]. Our aim is to develop an approach to genetic biodiversity that is computationally efficient and enables quantitative display of affinities at various taxonomic scales. Here we extend and refine this method and first apply it to large-scale differences, using sequences drawn from 12 diverse sets of animal species. On a finer scale we apply this mathematical apparatus to delineate affinities within one of the groups, North American birds, and examine biological implications of discontinuities that appear in structural representations of nucleotide sequence correlations.

Data Preparation
We considered the 648-nucleotide region of COI employed as a standard for distinguishing animal species [5]. Inspection of terminal regions of barcode sequence alignments deposited in (BOLD) http://www.barcodinglife.org showed a high degree of ambiguous or missing nucleotides, presumably reflecting incomplete sequencing runs. To reduce this noise we restricted attention to base pair (bp) positions 100 through 600 in the downloaded alignments, a 501nucleotide span representing 167 complete codons.
For the correlation analysis of the present framework nucleotide positions that are conserved lead to an uninformative increase in correlation, i.e., these carry no differential information. Among the 16,876 sequences of the 12 groupings considered below, we found that 161 of the 501 positions were conserved (Table 1); for the purposes of this analysis, these were dropped from analysis.
The stretch of 501 nucleotide characters can each be uniquely translated into a digital vector under the nucleotide convention as follows.
In schematic form a sequence s transforms to a vector s as follows s~:::GATTC:::?½:::G,A,T,T,C,:::? ½:::0,0,1,0,1,0,0,0,0,0,0,1,0,0,0,1,0,1,0,0,:::~s: ð2Þ There are various metrics for calculating sequence distances based on models of nucleotide substitution. Among these the Hamming distance, d H , i.e., the number of substitutions required to bring two sequences of like length into agreement, is the freest of additional assumptions. More complex distances distinguish between transitions and transversions, codon positions, and equilibrium based frequencies, as for example [7][8][9][10]. These forms are based on evolutionary considerations, while for our approach, which is based on the present state of correlations, the Hamming distance is the metric of choice. Each COI sequence thus becomes a vector s of 2004 entries; after removal of the 161 conserved nucleotides, 1,360 entries remain. The transformation of eq. (2) is not unique. An alternate transformation is which doubles, instead of quadrupling the sequence length as in eq.
(2). This does not lead to the desired form of eq. (5) given below. Other alternatives that have been tried also lead to problems.

Distances
If two sequences of length N disagree at K positions the Hamming distance is On the other hand from eq. (2) the square of its Euclidean distance d E is and therefore d H~d In normalized form this can be written as which places the sequences vectors measured from a zero origin on the unit sphere and also uniquely associates the correlation coefficient cos h, and the angle h, as a consequence of the law of cosines, i.e., the right hand side of eq. (7). d H is the ratio of substitutions to site number, a customary representation. Equations (6) & (7) are special cases of a more general recipe for associating a correlation coefficient with a metric. If d(,) denotes a metric (distance function), then we recall that for elements A, B & C by definition the triangle inequality is satisfied where c~d(A,B),a~d(B,C),b~(A,C): One may then show from eq. (8) that which fulfills the requirement of a correlation. And if the ratio in eq. (10) is written as cos h we obtain the law of cosines.
In a vector space this is exactly the case. In the construction eq. (7) C is taken as the origin.

Indicator Vectors
For purposes of exposition consider the particular grouping of ''Canadian freshwater fish'' see Table 2. After the above preparation of sequence data we denote a typical fish sequence by the row vector s(F ). The Canadian fish dataset has 1,324 members. Next we chose M distinct sequences at random from this set and form the fish set.
In general if there are N groupings we consider N sets G a~f s j (a)g, where a ranges over the N groupings.
An indicator unit vector v j for each set is then determined on the basis that it have a maximal correlation with the selected j th taxon, and minimal correlation with all other taxa [6]. As a simple but useful illustration consider N sequence vectors s j , j~1,::,N, say one representative from each of N groups, or each an average of each group. We then seek v l , the l th indicator vector such that where vw signifies the average. It is straightforward to show that under the reasonable assumption that if fs j g are linearly independent then the criterion function C l has a positive maximum and that it is determined as the eigenvector with the largest (positive) eigenvalue of One consequence of the particular criterion for choosing the v j is that it provides a natural structural representation expressed as auto-and cross-correlations, given by and referred to as the structure matrix. We also define the diversity matrix as given by This notation denotes the mean over all M(M{1)=2 inner products pairs of the members of G i with those of G j , which thus gives a depiction of within and among group correlations.
A fixed number of members, M, in the sets G confers equal weights on each of the taxa. These may be considered as the ''training set,'' for the indicator vector and the remaining sequences are used as a ''test set.'' There is reason to make M relatively small in initial calculations. Once past the testing stage there may be reason to take M as large as possible within the restriction of equal weightings.

Probabilities
Another consequence of embedding a character sequence into a vector space, eq. (2), is that the average of an ensemble of sequences fs j g can be defined as Which through the inverse operation of eq. (2) furnishes the probability of occurrence of (A,T,C,G) at each nucleotide position and that which is a consequence of eq. (2). Conservation of Probability. Eq. (4) allows us to regard the 4-vectors as specifying the probabilities of the associated symbols. We now demonstrate that this property is inherited by the indicator vectors, i.e., its 4-vectors sum to unity. To see this define where p, the number of rows, is also the number of bps. Multiplication of (15) by U yields but Us k~u~1 1 . . .    17). Numerical forms of matrices given in Table 3. Differing color bar scales in (A) and (B) are used to emphasize off diagonal resemblance between matrices. doi:10.1371/journal.pone.0009266.g001 Table 3. Numerical representations of Figure 1A and 1B, respectively. for any k and from this it follows that which proves the assertion. (This proof depends specifically on regarding an unknown bp as ½1=4,1=4,1=4,1=4, which we deem to be reasonable.) Therefore each indicator vector can be regarded as p quartets of probability in the four possible symbols.

Tree Construction
A customary practice is to express sequence separations as distances, which play a role in the construction of trees. It is straightforward to show the connection of distances to the correlations contained in eq. (16) and of eq. (17). In fact it directly follows from eq. (7) that is the matrix of average Hamming distances between taxons i and j. By the same token is the distance matrix between the i & j indicator vectors. It is important to note that evolutionary considerations do not figure in the calculation of the above distances.

Results
We first considered 12 animal groups, using COI sequences deposited in BOLD taxon-specific projects ( Table 2). In all cases analysis was restricted to sequences of sufficient length, and excluded those containing excessive blank positions.
The structure matrix for the 12 groups displays correlations among their respective indicator vectors ( Figure 1A). These are arranged in large-scale taxonomic divisions [Chordata, Arthropoda (Insecta, Malacostraca), Mollusca, Cnidaria], and sub-ordered based on correlations, e.g., within the upper 4|4 matrix (Chordata), groups are ordered by vector correlation as quantified by Thus amphibians have the highest relationships with the others in this set. The next 5|5 block representing Class Insecta, are ordered by relationship as above. The diversity matrix eq. (17) quantifies the degree of diversity within and among data sets ( Figure 1B), and has an impressionistic similarity to the structure matrix of unitary indicator vectors ( Figure 1A). The diagonal elements of Figure 1B illustrate the high internal diversity of amphibians, ants, crayfish, and jellyfish, and relative lack of internal diversity for flies, butterflies, wasps, and aphids. Numerical equivalents of Figure 1 are given in Table 3. Lack of diversity might be consistent with these data being drawn from single families or subfamilies. Diversity as defined by (24) introduces an objective measure of diversity based on variance.
We applied the 12 indicator vectors to the remaining set of 10,876 test sequences, generating a structure matrix of correlations, (Figure 2). With one interesting set of exceptions, there were no assignment errors, i.e., each individual test sequence was most highly correlated with its respective group-level vector. The exceptions were 33 sequences, .09% of all sequences, in the fish dataset which, according to the metric, were more closely correlated with the amphibian than the fish indicator vector. Inspection revealed that each error was caused by a lamprey (Class Cephalospidomorphi) sequence and all lamprey sequences produced this erroneous assignment. The remaining sequences in the Canadian fish dataset represented ray-finned fishes (Class Actinopterygii). Viewed taxonomically, the lampreys appear to be inadvertently included in fish dataset; when removed there was 100% accuracy of assignment of test sequences plus training sequences.
We applied the indicator vector approach at a finer scale, analyzing differences within the dataset of North American birds, which contained 1,693 sequences representing 558 species. As a compromise between a large M and a large test set, we chose M~3, giving 262 admissible species and 471 test sequences. With the input ordered alphabetically by taxonomic genus, the resulting structure matrix appears to be disordered with small regions of high correlation ( Figure 3A). When arranged in a taxonomic order representing phylogenetic relationships [11] (Table 4), these correlations coalesced into a coherent picture ( Figure 3B), which could be viewed as taxonomy organizing the structure matrix according to closeness of correlations. Discontinuities in the correlation among North American birds, evident as ''boxes'' or ''blocks'' in the color matrix, corresponded to avian taxonomic divisions (Figure 4). Most of the blocks represented families, with some blocks corresponding to lower (genera) or higher (suborder) groupings ( Figure 4).
Among the 471 test bird sequences, there were 16 apparently incorrect assignments distributed among 4 species pairs (Junco phaneotus/J. hyemalis; Anas platyrhynchos/A. rubripes; Larus smithsonianus/L. glaucescens; Sphyrapicus ruber/S. nuchalis). In the first instance each sequence set of M~3 were identical so that the indicator vectors were also identical. In the remaining cases the indicator vectors were close but not equal reflecting the fact that the defining sequence sets shared some identical members. While such singular behavior is revealed by the present algorithm, these sets of species were previously noted to be indistinguishable by COI barcode [12].
As indicated in eq. (24) the structure matrix S ij can be directly associated with a matrix of inter-species distances D ij . Since such a matrix can be made the basis of tree constructions we can apply the neighbor joining (NJ) algorithm of Saitou and Nei [13] to D ij . Using consistency arguments [14,15], Bryant [16] has demonstrated that the NJ construction is a unique clustering algorithm of the distance matrix [17]. Since the distance matrix is based on genomic distances, and not on evolutionary hypotheses, we can view the resulting NJ tree as intrinsic to the data. The species ordering according to this tree produces the structure matrix shown in Figure 5. This demonstrated the same set of clusters as seen in Figure 4; only the order of clusters differed. Thus at this level of resolution the indicator vector approach to classification coupled with NJ provides a self-generating ranking that is in general agreement with established taxonomy. Figure 6 compares the NJ tree that emerges from the structure matrix with the tree that derives from the averaged Hamming distance matrix between species, is equivalent to the diversity matrix (17).

Discussion
This paper describes a mathematical approach to comparative analysis of nucleotide sequences using digital transformation in vector space. We term the resulting structure matrices ''Klee diagrams'', in acknowledgement of the geometric paintings of artist Paul Klee (see Figure 7). This approach is of general utility and could be applied to any set of aligned sequences. In this study we explore its potential by analyzing a large, diverse set of DNA barcodes, the short segment of mitochondrial COI gene employed as a standard for identification of animal species (6). The resulting Klee diagrams display the structure of present-day mitochondrial genetic diversity, a ''macroscopic'' view of the products of evolution [18,19]. This approach is akin to a distance metric (see Methods), and in fact the matrix of indicator vector correlations can be used to generate an NJ tree ( Figure 6).
As compared to standard distance metrics with neighborjoining, indicator vectors preserve character probabilities that distinguish sequence sets, enable automated classification of test sequences, and generate high-information density displays without constraints of tree diagrams. Regarding the latter point, as one example, the 12-group Klee diagram displays affinity among flies and crayfish, a finding which might be of interest for further exploration, and yet this sort of horizontal similarity is not represented in the NJ tree diagram, shown in Figure 8. Discontinuities in indicator vector correlations, evident as blocks in Klee diagrams, corresponded to branches in the tree; for example, in North American bird matrix, these blocks represent families, genera, and sub-orders (Figures 4, 5). These results, generated with a small sample of world birds, suggest that this approach might be usefully applied to generate a classification for poorly-studied groups by combining DNA barcodes with indicator vector analysis. Such a classification could be refined when additional morphologic, ecological, and genetic study was available.
The results so far suggest natural discontinuities, or fractures, in the genetic structure of biodiversity, at least as reflected in animal mitochondrial genomes. In quantitative terms, blocks represent higher correlation within than among sets of sequences. Further study will help determine the nature of underlying mitochondrial differences, for instance whether species-and family-level blocks, for example, reflect differences in coding or non-coding positions. The present-day discontinuities seen in Klee diagrams may not be evident from a historical perspective, such as in a phylogenetic tree which links all forms in a continuous structure. It is of interest to reconcile these two perspectives, namely the continuous nature of evolution with the fractures in present-day genetic biodiversity;  these might be viewed, respectively, as ''time-like'' and ''spacelike''. One may speculate on the relation of such jump phenomena to adaptive radiations and the punctuated equilibrium model of evolution [20]. It may be possible to make useful observations for time-like behavior from space-like behavior as was done through the ergodic theory of statistical physics, [21]. As currently developed, our approach is limited to complete sets of homologous sequences, rather than overlapping sets of incomplete data as are often used in phylogenetic inference. In addition, the present analysis employing COI shares problems inherent to mitochondrial biology, including maternal inheritance, introgression, hybridization, male-biased dispersal patterns, and recent speciation among others [22]; most of these are likely to apply only at the fine-scale level of distinguishing closely-related species. As noted, the indicator method is of general utility and could readily be applied to longer sequences or concatenated multi-gene alignments without substantially increasing computation time, which might address some of these limitations. In this regard, it of interest to compare indicator vector affinities using mitochondrial and nuclear genes in puzzling cases that appear to represent convergent evolution [23].
Although the output is different, it may be revealing to compare the efficiency of the indicator vector approach to that of phylogenetic treebuilding programs. Due to computational demands, data sets in analyses beyond 1000 species are exceptional (e.g., [24][25][26]) and calculation times for larger studies are typically several CPU-months. The largest published phylogenetic tree includes 73,060 eukaryote taxa [1] and took 2.5 months with 16 processors, and the next largest analyzed 13,533 plant taxa [27]. The present study ranks with the largest biodiversity analyses in terms of number of organisms, and is at least two orders of magnitude faster. For example, the case of 12 animal groups deals with almost 17,000 sequences and required times of roughly 10-20 minutes on an ordinary desktop computer. This suggests the potential for analyzing the largest datasets available, including, for example, BOLD (w700,000 sequences) http://www.barcodinglife.org, NCBI Influenza Virus Resource (w50,000 complete genomes) http://www.ncbi.nlm.nih.gov/ genomes/FLU/FLU.html, or Los Alamos HIV Sequence Database (http://www.hiv.lanl.gov).
In addition to animals, cytochrome c oxidase is present in plants, protozoa, fungi, and some bacteria, which raises the prospect of insight into broad-scale patterns in the genetic structure of biodiversity. Also, the methodology as present here applies to nucleotide sequences of any sort and so might usefully be applied to a variety of questions.
From the point of view of accuracy, density of information and assimilation it would seem compelling that any properly ordered distance matrix should be viewed as a Klee diagram. It may be that the focus on evolution and therefore trees impeded this direction. In this connection we note that the distance matrix for a species count of N contains N(N{1)=2 distances and for large N a treebuilding algorithm cannot accommodate this number of condi-   tions, and an increasing number of larger and larger errors occur with increasing N. Klee diagrams accurately display distances for any species count.
An important advance in the present treatment derives from the vectorization of nucleotide sequences, (1), which has been accomplished with the exact preservation of Hamming distances. Advantages flow from a vector space framework, an example of which is the optimization procedure leading to the indicator vectors. Another consequence is that bps occupation is rigorously transformed to the probability of occurrence of the four nucleotides, which opens the possibility of introducing information theory into these considerations.