OutbreakFinder: a visualization tool for rapid detection of bacterial strain clusters based on optimized multidimensional scaling

Multidimensional scaling method

Multidimensional scaling is a widely used approach to projecting high-dimensional data onto a low-dimensional space, where the relationships among individual samples are preserved. This characteristic of MDS is often used to present high-dimensional data on a two-dimensional plane. A simple derivation of the MDS method is presented as follows:

Given n points in a p-dimensional space, the Euclidean distance between any two points, x_r and x_s, can be defined as $d_{rs}^2=\sum _{i=1}^p{\left(x_{ri}-x_{si}\right)}^2$. The MDS method aims to minimize the objective function as follows: $$\text{tress=}\sqrt{\frac{{\displaystyle \sum {}^{{\left(d_{rs}-{\widehat{d}}_{rs}\right)}^2}}}{{\displaystyle \sum {}^{{}^{d_{rs}}}}2}},$$

where ${\hat{d}}_{rs}$ is the predicted Euclidean distance. Let B = XX^T be an n × n inner product matrix. In metric MDS, matrix B can be decomposed into QMQ^T, where Q is an eigenvector and M is a diagonal matrix of eigenvalues. The two forms of matrix B are combined as follows: $$QM{Q}^T=Q\sqrt{M}\sqrt{M}{Q}^T={\left(\sqrt{M}{Q}^T\right)}^T\left(\sqrt{M}{Q}^T\right)=X^{T}X$$

Then, we obtain $X=\sqrt{M}{Q}^T$. Therefore, the aforementioned derivation demonstrates that the traditional MDS method is designed to find an approximate solution rather than the optimal solution. For detailed derivation of the traditional MDS method, please refer to Supplemental Appendix A.

MDS by using Newton’s method

The 100 points in Fig. S1 are generated by the function genSimulationData() as shown in Supplemental Appendix B, every 10 points is a group, 0–9, 10–19, 20–29…, and so on. Figure S1A is an example of the distortion produced by the traditional MDS method, as shown in Supplemental Appendix C. Points 0–9 belong to the same group; however, points 9 and 3 evidently do not belong to the same group in traditional MDS (Fig. S1A). Therefore, we propose to perform MDS by using Newton’s method (NMDS) to compensate for the deficiencies of the traditional MDS method. The proposed method has two parts; one is NMDS and the other is determination of the largest error edge to be amended. The detailed steps of the algorithm are as follows:

1. ${\hat{d}}_{jk}=\sqrt{\sum _{i=1}^p{\left(x_{ji}-x_{ki}\right)}^2}$

2. $r_{jk}=d_{jk}/{\hat{d}}_{jk}$, where d_jk is the observed distance

3. ${\hat{x}}_{ji}=\frac{1}{n}\sum _{k=1}^n\left(x_{ji}-x_{ki}\right)\left(r_{jk}-1\right)$

4. Find the maximum error edge and force it to move to an ideal position (find the maximum ratio of ${\hat{d}}_{jk}/d_{jk}$ as the maximum error edge)

Here, ${\hat{d}}_{jk}$ represents the predicted distance between nodes j and k; r_jk represents the ratio of the observed and predicted distance between nodes j and k, namely $d_{jk}/{\hat{d}}_{jk}$; ${\hat{x}}_{ji}$ represents the new value of dimension i; and x_ji represents the current value of dimension i.

Steps 1–3 in NMDS minimize the objective function. Step 4 finds the maximum error edge and forces it to move to an ideal position. As shown in Fig. S1B, the proposed method can reduce grouping errors.

Multidimensional scaling by using Newton’s method is implemented in OutbreakFinder, which provides parameters through which a user can change the distance scale between data points, such as using a logarithmic scale or power 2 scale for distances. Users can specify the color of each data point or use the AP algorithm to cluster and specify the color of each group. In addition, OutbreakFinder generates a text file of data point coordinates and color, as well as an MDS plot graph. Users can reproduce the MDS plot figure by using other tools. OutbreakFinder can directly read Lyve-SET, Parsnp, and multiple sequence alignment files and plot MDS graphs. If a user does not use the aforementioned tools, he or she can generate a distance matrix to plot an MDS graph.

Implementation

OutbreakFinder is written in Java and compiled into a standalone executable jar file that can be executed in Java Runtime Environment 1.8 or a later version. In addition, users can download the source code and compile it into a preferred version.

OutbreakFinder can parse the results of the Lyve-SET, Parsnp and Multi Alignment tools and generate MDS coordinates and graphs, as shown in Fig. 1. The user can define the color of each data point, and OutbreakFinder will export the corresponding image file according to the defined color map. Users can also use the AP module in OutbreakFinder for automatic clustering. The advantage of using an AP module for clustering is that there is no need to provide the number of groups, which makes cluster analysis easier.

Example analysis

To demonstrate how to use OutbreakFinder, we employed 59 Salmonella Heidelberg isolate genomes from three outbreaks with the same PFGE type (Bekal et al., 2016). Whole-genome sequencing raw reads (Table S1) were downloaded from the National Center for Biotechnology Information’s Sequence Read Archive database (https://www.ncbi.nlm.nih.gov/sra). The downloaded sra files were converted into fastq files; then, Lyve-SET was used to extract SNPs from these fastq files and produce a distance matrix. OutbreakFinder uses the distance matrix to generate an MDS plot, as shown in Fig. 2A. Subsequently, we used the same dataset to test Parsnp. OutbreakFinder can directly parse the results of Parsnp output and generate an MDS plot, as shown in Fig. 2B. We use published benchmark datasets (Timme et al., 2017) to verify the availability of OutbreakFinder. This dataset contains four major foodborne bacterial pathogens, such as Listeria monocytogenes, Campylobacter jejuni, Escherichia coli, and Salmonella enterica. These pathogens are listed in Tables S2–S5. The procedure is the same as described above, Lyve-SET was used to extract SNPs and produce the distance matrixes, and then OutbreakFinder uses the distance matrixes to export the MDS plot. As shown in Fig. 3, outbreaks and outgroups can be clearly distinguished from the MDS plot. Although there are two outbreak clusters in Fig. 3A, the outbreak and outgroup can still be clearly distinguished.

Performance

Finally, we used a large dataset to stress OutbreakFinder, which included 38,360 16S rRNA sequences downloaded from rrnDB (Stoddard et al., 2015). Although 38,360 sequences constitute a very large dataset, the results revealed that OutbreakFinder could still output results normally, as shown in Fig. 4. By contrast, if the traditional MDS method is used for dimension reduction analysis, more computational resources are required because the time complexity of traditional MDS is O(n³). In the case of NMDS, the time complexity is only O(cn²), where c is the number of iterations and usually does not exceed 1,000. Therefore, OutbreakFinder can calculate the MDS of 38,360 16S rRNA in 1 day using a single core and 30 G memory computing resources. In contrast, the manifold.mds cannot handle the same dataset, even if the computing resources are increased to 20 core and 120 G memory, the work cannot be completed. Manifold.mds is a method in the scikit-learn suite in Python. To confirm that the proposed method is robust enough, we performed 1,000 simulations to compare the performance of NMDS and manifold.mds. As shown in Fig. S2, in most of the simulated fn is smaller than fc, fn, and fc are the values of the error function of NMDS and manifold.mds, respectively, which means that the performance of NMDS is better than that of manifest.mds.