Not Seeing the Forest for the Trees: Size of the Minimum Spanning Trees (MSTs) Forest and Branch Significance in MST-Based Phylogenetic Analysis

Trees, including minimum spanning trees (MSTs), are commonly used in phylogenetic studies. But, for the research community, it may be unclear that the presented tree is just a hypothesis, chosen from among many possible alternatives. In this scenario, it is important to quantify our confidence in both the trees and the branches/edges included in such trees. In this paper, we address this problem for MSTs by introducing a new edge betweenness metric for undirected and weighted graphs. This spanning edge betweenness metric is defined as the fraction of equivalent MSTs where a given edge is present. The metric provides a per edge statistic that is similar to that of the bootstrap approach frequently used in phylogenetics to support the grouping of taxa. We provide methods for the exact computation of this metric based on the well known Kirchhoff’s matrix tree theorem. Moreover, we implement and make available a module for the PHYLOViZ software and evaluate the proposed metric concerning both effectiveness and computational performance. Analysis of trees generated using multilocus sequence typing data (MLST) and the goeBURST algorithm revealed that the space of possible MSTs in real data sets is extremely large. Selection of the edge to be represented using bootstrap could lead to unreliable results since alternative edges are present in the same fraction of equivalent MSTs. The choice of the MST to be presented, results from criteria implemented in the algorithm that must be based in biologically plausible models.


Introduction
The use of trees for phylogenetic representations started in the middle of the 19th century. One of their most popular uses is Charles Darwin's sole illustration in "The Origin of Species" [1]. The simplicity of the tree representation makes it still the method of choice today to easily convey the diversification and relationships between species. Many different methods have been proposed to reconstruct phylogenies, mostly concerned with recovering evolutionary relationships over long periods of time [2]. Each algorithm or method used to infer and draw a tree, makes a series of implicit or explicit assumptions that limit the types of trees generated. This variability in the trees generated by different algorithms using the same data, has important repercussions that frequently go unappreciated by those who use them. At shorter timescales and with limited diversity, conditions that are encountered in population genetics and microevolutionary studies of a single species, the assumptions made by these methods may not be equally valid [3] and a number of other methods have been used when analyzing this data.
Minimum Spanning Trees (MSTs) are becoming increasingly used for representing relationships between strains in epidemiological and population studies of bacterial pathogens. Although MST computation is a classical mathematical problem and its application to evolutionary studies was suggested more than a decade ago [3], it was not until recently, with the advent of multilocus sequence typing (MLST) [4] and particularly whole genome sequencing [5,6], that they gained popularity as an alternative to eBURST [7]. One appeal of MSTs is the simplicity of their assumptions that reflect the concept of minimal evolution. MSTs simply link together the more closely related individuals in the population, generating a single tree representing all individuals. The Steiner trees [8], generated by the more classical methods for phylogenetic inference, place individuals exclusively in branch tips. By allowing individuals to be placed in interior nodes, spanning trees and MSTs in particular, may better convey the peculiarities of short-term intraspecific evolution [3].
It was also recently pointed out that the optimal implementation of the BURST rules in goe-BURST, results in a set of disjoint MSTs [9]. These trees group sequence types (STs) that differ by a maximum threshold number of alleles from at least one other ST in the group. These groups or connected components are frequently referred to as clonal complexes (CCs). In fact, goeBURST is a maximum weight problem that together with MSTs are particular cases of graphic matroids [9]. But, as it is well known, MSTs are in general not unique for a given network and this was also recognized in the context of phylogenetic trees [3,10]. The fact that a single tree is reported from a multitude of possible and equally optimal solutions and that no statistical metrics exist to evaluate them, justified a recent heuristic approach to address these issues [10]. The authors suggested a method based on a mark-recapture approach to estimate the number of possible trees and a bootstrap procedure to evaluate tree credibility.
The problem of counting MSTs has been a challenge for the last decades, namely the development of efficient approaches for counting MSTs in weighted graphs, and different approaches have been described. In 1987, Gavril [11] addressed the problem of counting the number of MSTs by constructing a tree-like recursive structure, the root of which is the subgraph G 0 formed by removing all non-maximum-weight edges from G, and each sub-tree of which is constructed recursively from the components of G\G 0 . The minimum spanning trees of G can then be counted by multiplying together the numbers of spanning trees at each node of this structure. This method runs in O(nM(n)) time, where M(n) is the time required to multiply two n × n matrices. Later, in 1997, Broder and Mayr [12] improved this bound by proposing a method based on a generating function that can be expressed as a simple determinant, where the weights of the edges appear as exponents of polynomials. This method proceeds by factoring the determinant and it works for nonnegative integral edge weights. It runs in O(M(n)) time.
Eppstein [13] took still a different approach and created the concept of equivalent graph. Specifically, one constructs from any given edge-weighted graph G an equivalent graph EG without weights, with a sliding transformation, such that the minimum spanning trees of G correspond one-for-one with the spanning trees of EG. Having translated the weighted graph to an equivalent unweighted graph, one can compute the number of MSTs by just applying the Kirchhoff's matrix tree theorem to the new graph.
Note that most of these approaches aim at generating and sampling MSTs, a harder problem than just counting the number of MSTs. Moreover, our approach may be applied to the general case of graphic matroids. As discussed above, the problem of finding an MST is a particular case of graphic matroids [14] and, thus, finding a solution for a given graph consists of solving an instance of graphic matroids [14][15][16], which can be optimally solved with a greedy approach [17]. One of those greedy approaches is precisely Kruskal's algorithm [18]. In the general case of graphic matroids, edges may be unweighted, which is usually the case. We just need to define a total order for the edges based on specific criteria, which is precisely what we have in phylogenetic approaches using MSTs [9].
Here, we present an improvement to the mark-recapture/bootstrap approach by introducing a new edge centrality metric taking advantage of determining exactly the number of possible trees and the proportion of the tree space that includes each of the possible edges through an expansion of the Kirchhoff's matrix tree theorem [19,20]. Contrary to other methods that depend on edges being weighted, our method just depends on sorting edges in increasing order and, thus, we just require a total order to be defined. This simple approach allowed better performance and can be applied to the general case of graphic matroids.

Methods
In this section, we describe an exact method to compute the significance of a branch/edge in a given MST representation, and we present a module for the PHYLOViZ software [21] implementing the proposed metric. We start with the formalization of the problem under study and the proposed metric. Then, we show how the metric can be effectively computed.

Spanning edge betweenness
Let G = (V, E) be a connected, undirected and weighted graph, with weight function w: E ! IR, where V is the set of vertices and E & V × V is the set of edges. A Minimum Spanning Tree (MST) T = (V, E 0 ) is a subgraph of G that is a tree and contains all the vertices of G, i.e., that spans over all vertices in V, with jE 0 j = jVj − 1, and such that ∑ e2E 0 w(e) is minimum among all possible spanning trees. It is clear that we can have more than one MST for a given graph G and we would like to count how many MSTs exist in G. The solution to this problem is provided by the Kirchhoff's matrix tree theorem [20] for unweighted graphs and by Eppstein [13] for weighted graphs, where the Kirchhoff's matrix tree theorem is still used but only after some graph transformations.
However, in this paper we are interested in a slightly different question. Given an edge e 2 E, we want to know the fraction δ G (e) of MSTs where e occurs. The value δ G (e) is what we call the spanning edge betweenness for e and it is formally defined as where τ G is the number of different MSTs for G and τ G (e) is the number of different MSTs for G where e occurs. Note that δ G (e) may be zero whenever an edge e is not present in any MST, causing τ G (e) to be zero. In what follows we write δ(e), τ(e) and τ whenever G is clear from the context. It remains to see how to compute, as efficiently as possible, the spanning edge betweenness τ G (e) for a given e 2 E. In next sections, we show how to compute τ G (e) and δ G (e) when G = (V, E) is a connected, undirected and unweighted graph, with n = jVj vertices and m = jEj edges. Note that in this case the number τ of MSTs in G is equal to the number of spanning trees in G and it can be computed directly from the Kirchhoff's matrix tree theorem [19]. Then we extend our result to weighted graphs and we discuss implementation details.

Unweighted graphs
Let F 2 {−1, 0, 1} n × m be the incidence matrix for G such that F i, e = 1 and F j, e = −1, for e = (i, j) 2 E where i < j without loss of generality. Let us also consider the reduced incidence matrix F (i) obtained from F by deleting row i. Note that rank(F) = n − 1, rank(F (i) ) = n − 1, and the determinant for any square submatrix of F (i) , for any i, is either 0, −1, or 1. A more interesting observation due to Kirchhoff is that a submatrix (n − 1) × (n − 1) of F (i) , for any i, is non-singular if and only if its columns correspond to the edges of a spanning tree.
Theorem 1 (Kirchhoff [19]). The spanning trees of a connected and undirected graph G with n vertices are the non-singular (n − 1) × (n − 1) submatrices of the reduced incidence matrix F (i) , for any i, and the determinants of the submatrices are all ±1.
Hence, by using Cauchy-Binet theorem on determinants, the number of spanning trees τ is given by the Kirchhoff's well known formula where S ranges over the subsets of E with size n − 1, L = FF > is the Laplacian matrix for G, and L (i) denotes the matrix obtained from L by deleting row and column i. We extend this result to compute τ(e), for e 2 E, as follows. Theorem 2. Given G = (V, E) an undirected and connected graph, let e = (i, j) 2 E and L (ij) denote the matrix obtained from L by deleting rows i and j and columns i and j. Then, det(L (ij) ) is the number of spanning trees τ(e) that contain e.
Proof. As discussed above, the total number of spanning trees is given by det(L (i) ), for any i. Let G 0 be the graph where we remove the edge (i, j) and L 0 be the Laplacian for G 0 . Hence, the total number of spanning trees for G 0 is given by det(L 0 (i) ), for any i, and the number of MSTs that contain (i, j) is simply given by det( ). Let us show that det(L (ij) ) = det(L (i) ) − det(L 0 (i) ) or, equivalently, that det(L (i) ) = det(L 0 (i) ) + det(L (ij) ). We have that is obtained from F by removing rows i and j, and, using Cauchy-Binet's formula, we can show instead that where F 0 is the incidence matrix for G 0 , S ranges over the subsets of E with size n − 1, S 0 ranges over the subsets of E\{(i, j)} with size n − 1, and S Ã ranges over the subsets of E with size n − 2. Since S 0 ranges over the subsets of E\{(i, j)}, we can replace F 0 by F in previous equation. Note also that because adding edge (i, j) to S Ã and considering F (i) instead of F (i, j) just adds a term ±1 to each matrix determinant. Therefore, which is an equality as the first term on the right side ranges over all subsets of E with size n − 1 that do not contain (i, j) and the second term ranges over all subsets of E with size n − 1 that do contain (i, j). Hence, using both results, we can easily compute δ(e) for any e 2 E. Note also that the same is true for multigraphs, graphs that allow multiple edges between the same pair of vertices, as both results above hold with the following changes in the Laplacian matrix L [22]: if vertex i is adjacent to vertex j in G, then L ij is equal to the number of edges between i and j; when counting the degree of a vertex, all loops are excluded.

Weighted graphs
Let G = (V, E) be a connected, undirected and weighted graph, with weight function w: E ! IR. We can compute a MST for G by using the Kruskal's algorithm [18]: 1. sort E with respect to w in increasing order; 2. create a forest M where each u 2 V is a tree; 3. iterate over E in increasing order and, for each (u, v) 2 E, if u and v are in different trees, add (u, v) to M combining both trees as single tree; 4. return M.
Note that we may get different MSTs by changing the order obtained in step 1, where we can exchange positions of edges with the same weight. In particular, since it is well known that the sorted list of edge weights is the same for any MST, changing the order allow us to obtain all different MSTs.
We can take this a step further. Consider the algorithm SEB for computing the number of MSTs and the spanning edge betweenness for each edge: The algorithm SEB works similarly to the Kruskal's algorithm by iterating over edges in increasing order with respect to w and, at each main iteration (steps 3 to 6), it considers sets of edges with the same weight. Let e 2 E and let M 0 be the forest obtained in Kruskal's algorithm after processing all edges e 0 2 E such that w(e 0 ) < w(e). Let also H be a graph where each tree in M 0 is a vertex, i.e., where each tree was contracted, and where we add all edges in E with weight equal to w(e). Since for some main iteration of SEB algorithm we stop after adding edges e 0 2 E such that w(e 0 ) w(e), H does not contain edges e@ 2 E such that w(e@) > w(e). Moreover, since we contracted all edges e 0 2 E such that w(e 0 ) < w(e), all edges in H have the same weight w(e), and we can treat it as an unweighted graph (or, since H may be a multigraph, as an unweighted multigraph). Hence, if we consider the connected component C of H that contains edge e, and by using results in previous section, we are able compute the number τ C of spanning trees for that component and also the number τ C (e) of spanning trees for that component where e occurs. The key observations clarified in the following lemmas are that we can use this approach to compute the number of spanning trees in G and that δ G (e) = δ C (e).
Lemma 1. Given G = (V, E) a connected, undirected and weighted graph, with weight function w: E ! IR, the algorithm SEB computes the number of spanning trees in G.
Proof. It is clear that an edge e 2 E can only permute with another edge e 0 2 E to form a different MST iff w(e) = w(e 0 ) and, if a MST M contains e, adding e 0 to M leads to a cycle. Moreover, that cycle can only contain edges with weight equal to or lower than w(e), otherwise M would not be an MST. If we add all edges with weight w(e) to M and contract all edges with weight lower than w(e), we obtain the graph H and the product of the number of trees in each connected component of H is the number of ways we can select edges with weight w(e) for each MST of G. By doing this for each different weight in G and then multiplying all values, we obtain the number of MSTs τ for G.
Lemma 2. Given G = (V, E) a connected, undirected and weighted graph, with weight function w: E ! IR, an edge e 2 E, H the graph obtained in algorithm SEB while processing edges with weight equal to w(e), and C the connected component of H that contains e, we have that δ G (e) = δ C (e).
Proof. Since a given edge e only has influence on the number of trees for the component C of H where it occurs, the number of trees for all other components C 0 in H, and in any other graph H in remaining algorithm iterations, remains the same. In particular, by inspecting algorithm SEB, and, by a similar construction, where the last equality holds because edge e 2 E does not occur in any C 0 . Therefore, by Equa-

Implementation in PHYLOViZ
We have implemented our metric as a module for PHYLOViZ [21], available at http://www. phyloviz.net/. Our implementation uses the Colt library (http://acs.lbl.gov/software/colt/) for linear algebra operations, including in particular the computation of matrix determinants. Since we are dealing with relatively large sparse graphs, we use the class SparseDoubleMatrix2D in Colt. We also use a disjoint-set data structure to track connected components similarly to what is common in Kruskal's algorithm implementations [23].
The time complexity of the proposed approach is dominated by the time required to compute the determinants, since the Kruskal's algorithm runs in O(m log n) time, for a graph with n vertices and m edges. Computing the determinant for a n × n matrix can be done in O(n 1.5 ) time [24]. Hence, for sparse graphs with m = O(n), this method runs in O(n 2.5 ) time since we have to compute a determinant for each edge. In practice, it runs faster as connected components are usually much smaller than the original graph.

A more efficient implementation
Beside the motivation of an implementation of a module to PHYLOViZ application, we also implemented an offline version of the module where we used some extra settings to accelerate its execution and allow to obtain results that are not meant to be shown in PHYLOViZ. In this offline implementation, we used the MTJ library (Matrix Toolkit Java, available at https:// github.com/fommil/matrix-toolkits-java/) that is is a high-performance library for developing linear algebra applications. MTJ is based on BLAS (http://www.netlib.org/blas/) and LAPACK (http://www.netlib.org/lapack/) for its dense and structured sparse computations.
With this library, we use the LU decomposition to calculate the determinant. We create an upper triangle dense matrix and then we go through all the elements of the diagonal. Instead of multiplying all the determinants as in the module developed for PHYLOViZ we sum the logarithm of each absolute diagonal value, obtaining instead the logarithm of the determinant.
To improve our running time, we used the Java concurrent library for computing edges statistics in parallel. Since the statistics for each edge can be computed independently, we could parallelize statistics computation in a straightforward manner. Note that for computing statistics for each edge we must compute the determinant for a given matrix and, since these computations are independent, we can compute edge statistics in parallel. Package available at https://bitbucket.org/phyloviz/popsim-analysis.

Results and discussion
The spanning edge betweenness was applied to nine publicly available MLST databases of important human pathogens: Burkholderia pseudomallei, Campylobacter jejuni, Enterococcus faecium, Haemophilus influenzae, Neisseria spp., Pseudomonas aeruginosa, Streptococcus agalactiae, Staphylococcus aureus, and Streptococcus pneumoniae. These databases were retrieved on June 24th, 2014, from public repositories available in different websites (see Acknowledgments for more details). From all publicly available databases, we considered only those that generated graphs with more than 500 unique STs linked to at least one other ST at the single-locus variant (SLV) level. Analyses were performed both with PHYLOViZ, using a new module publicly available, and with a command line implementation developed to take advantage of high performance numerical libraries and of parallelization in multi-core platforms (see Methods for more details). We determined the goeBURST forest of each species by linking STs at SLV level, double-locus variant (DLV) level and triple-locus variant (TLV) level. Unless otherwise stated, the analyses were performed on the forest generated by creating trees linking STs at the SLV level. Details on how to reproduce this study, including copies of used databases, are also available at https://bitbucket.org/phyloviz/popsim-analysis.
We calculated the number of possible MSTs in the largest CC of each of these species ( Table 1). As expected, even only for the largest CC, the number of possible MSTs is quite large, in fact it exceeds a googol [25] for most of the species considered. When MST results are presented, a single tree is usually shown. This tree is chosen from among the space of possible trees, following a set of rules or simply as a consequence of the algorithm used and the input order of the nodes [10]. The goeBURST algorithm implemented in PHYLOViZ, selects the final tree according to a set of well defined rules that guarantee the uniqueness and consistency of the selected tree, independently of the input order of the nodes [9,21]. The impact of the application of each of the rules on the space of possible trees for the largest CC of each species is presented in Table 2. For most species, a single tree is obtained when applying up to the second tiebreak rule (higher number of DLVs), but in the case of B. pseudomallei, C. jejuni and Neisseria spp. a single tree is only obtained when invoking rules up to the third tiebreak rule (higher number of TLVs). In the case of S. pneumoniae only the last tiebreak rule (higher number of STID) results in a single tree. Large reductions in the available tree space can be seen with the application of each goeBURST rule and this can be used to evaluate the impact of each rule on the final phylogenetic hypothesis proposed by the algorithm.
The magnitude of the reduction of tree space varies between the species considered ( Table 2). It is clear that the number of STs influences the number of possible trees, with the Table 1. Statistics relative to the largest CC linking STs at the SLV level. Columns represent the number of STs, the number of edges, the total number of possible MSTs and the compactness and clustering indexes. The compactness index quantifies how directly connected individuals in the network are. The clustering index quantifies how close the neighbors of a given individual are from a complete graph (clique). Each index is an average after computing the index individually for each ST. number of possible trees increasing with the number of STs (Table 1). But this relationship is complex, with the number of possible edges linking STs at the SLV level having a similar and equally significant influence on tree space. For instance, when comparing the largest CCs of B. pseudomallei and E. faecium, although both have a similar number of STs, the latter has a higher number of possible edges and trees (Table 1). An even more striking example is the comparison between the largest CCs of S. aureus and S. pneumoniae, with the former having a smaller number of STs, but a higher number of possible edges and trees ( Table 1). The measurements of compactness and clustering of the tree of the largest CC capture properties that may be related to intrinsic characteristics of each species. These measures have several formulations in the literature [26], but for the results presented here, these were defined as follows. Given a graph G = (V, E), compactness describes how well a vertex u is connected in the graph, being defined as the quotient between the vertex degree d u and the maximum number of possible neighbors jVj − 1, i.e., d u /(jVj − 1). The clustering coefficient describes how well connected is the neighborhood of a vertex v, being defined as the quotient between the number of edges among neighbors N v of v and the maximum number of possible edges among them, i.e., 2j{(v,

Data sets
. The compactness and clustering coefficient of the graph are defined as the average of the vertex compactness and clustering coefficient over vertices in G, respectively. These definitions allow us to also compute these values for each connected component, which is of particular interest for the data under analysis. We note also that, although these two measures are related, they allow us to discriminate some interesting graph characteristics, which for the data under analysis may be related to mutation and recombination events. For instance, values of compactness < 0.010 are associated with B. pseudomallei, C. jejuni, Neisseria spp. and S. pneumoniae; species that also reach higher tiebreak rules to identify a single tree (Table 2). These species are known to have high rates of recombination [27][28][29][30].
The existence of recombination can generate STs with multiple possible pathways of descent, which in turn would be expected to affect a graph's compactness. The goeBURST algorithm in PHYLOViZ can be run by creating sets of disjoint trees linking STs at DLV or TLV level and the result of this analysis for the largest CC of each of the species considered is presented in Table 3. As expected, as we go from SLV to TLV, the higher number of STs and possible edges in the largest CC results in higher numbers of possible trees. The tree space at any given level, when considering the entire forest, is the product of the number of trees for each CC and is greatly influenced by the largest CC, hence our decision to present the analysis of the largest CC only for simplicity. We have previously proposed that the tiebreak rule reached before deciding if an edge should be drawn, could be used to evaluate the reliability of the represented hypothetical pattern of descent [9]. The spanning edge betweenness can be used for the same purpose, with results that are similar to those of the bootstrap procedure used frequently to support the grouping of taxa on trees [2]. In Fig. 1A we represent all possible edges that could be drawn between STs differing at a single locus (SLVs) in CC1439 of S. pneumoniae. While several STs are only linked by one possible edge to another ST, others are linked by several edges to a number of different STs of which they are SLVs. The goeBURST algorithm will then choose which edges should be represented in the final tree from among the edges found in the 88,833,024 possible MSTs of CC1439. In Fig. 1B is represented the MST identified using the goeBURST rules. On each edge is also indicated the percentage of the equivalent MSTs where that edge is found (the spanning edge betweenness). As expected, all the edges that were unique in Fig. 1A were found in all equivalent MSTs, such as the edge between ST6544 and ST4560. On the other hand, the represented edges of STs that could be linked by multiple possible edges, such as ST369, are present in a lower number of possible trees (in this case 33.3% or 29,611,008 trees). With a given set of STs, a higher proportion of equivalent MSTs including the represented edge, means that fewer alternatives are possible and this can be interpreted as a higher confidence in the represented edge.
Another important conclusion that can be drawn from calculating spanning edge betweenness is that all alternative edges linking a given ST to other STs are present in the same proportion of equivalent MSTs. For instance, ST1931 is linked by two possible edges, either to ST4519 or to ST547 (Fig. 1A), and both edges are represented in 66.7% of all equivalent MSTs (see the chosen edge by goeBURST in Fig. 1B). This means that the mark-recapture/bootstrap approach suggested previously to choose between alternative edges [10], cannot provide an adequate solution to the problem of selecting the most adequate edge. In the space of MSTs, alternative edges are equally represented and so any dominance of a given edge in the mark-recapture procedure will be a consequence of the limitations of the procedure and should not be used as selection criterion. The choice between alternative edges must be based on well defined criteria that should reflect an underlying model of microbial evolution. The goeBURST rules have such an underlying model [7,9,31] and offer a robust method to select a tree from a forest of MSTs. Table 3. Statistics relative to the largest CC linking STs at the SLV, DLV and TLV levels of construction. SLV means that the graph contains only edges linking STs that are SLVs, DLV means that edges were drawn between STs that are SLVs or DLVs of each other, and TLV that edges were drawn between STs that are SLVs, DLVs or TLVs of each other, according to the rules implemented in goeBURST.

Data sets
Statistics for the largest CC  In contrast to the number of MSTs discussed above, there is not a dominant role of recombination in determining the shape of the distribution, since the species identified previously as being recombinogenic are not homogeneous in their distributions. These differences possibly reflect differences in size of the data set considered, as well as a more complex interplay of the intrinsic properties of each species, such as mutation and recombination rates and possibly their ratio.
Centrality measures are important in a large number of graph applications, from search and ranking to social and biological network analysis [26]. Most of these measures are calculated upon the nodes/vertices. With node centrality we can measure the relative importance of nodes within a graph [32] but our interest can be to study the importance of links/edges on a network. A first approach was done by Girvan and Newman [33] where they defined edge betweenness, generalizing Freeman's betweenness centrality [34] to edges, as the number of shortest paths between pairs of vertices that run along an edge, with a direct application on the identification of community structures in networks. There are, however, other problems where alternative definitions of edge centrality are required, as is the case with the statistical evaluation of phylogenetic trees.
Here, we present a new edge centrality metric, the spanning edge betweenness, defined as the fraction of MSTs containing a given edge. We provide the required results and methods to exactly compute this metric. Since we rely on the Kirchhoff's matrix tree theorem, thus needing to compute several determinants for slightly different matrices, we plan to investigate how to accelerate these computations by reusing previous computations and by using more efficient methods for sparse positive semi-definite matrices decomposition, such as those based on Cholesky's decomposition [35]. Currently, our solution allows the algorithm to run on a common laptop in reasonable time (Table 4). For data sets with few STs, such as H. influenzae and P. aeruginosa (Table 1), it takes less than 2 seconds. However, for data sets with a larger number of STs, such as C. jejuni, it can take almost one hour. The running time will depend, mostly, on the number of STs of each data set, that is clearly related to the dimensions of the matrix representing the relationships between STs. Hence, the number of STs is directly related to the number of operations required to calculate determinants.
The comparison between this metric and other well known centrality metrics should also be investigated in the context of complex network analysis, as it provides a rather different approach for evaluating edge relevance or significance. The analyses of MLST data sets available in public databases show the usefulness of spanning edge betweenness in evaluating MSTs as proposals for phylogenetic relationships, and in providing confidence levels for each selected edge in the final tree. These analyses also highlight the impossibility of selecting an MST based on the statistical support of the edges, and reinforce the importance of the biological plausibility of the model underlying the criteria for edge selection in presenting the best possible MSTbased-proposal for the phylogenetic relationship of the entities under analysis. The use of bootstrap values became a key method to assess nodal support in phylogenetic trees [2]. The spanning edge betweenness proposed here offers a similar tool for the evaluation of MSTs in phylogenetic studies.