Rapid identification of pathogenic bacteria using Raman spectroscopy and deep learning

TL;DR

This study demonstrates a rapid, culture-free method for identifying bacterial pathogens and antibiotic resistance using Raman spectroscopy combined with deep learning, achieving high accuracy in clinical and laboratory settings.

Contribution

The paper introduces the largest bacterial Raman spectra dataset and applies advanced deep learning to accurately identify pathogens and resistance, advancing rapid diagnostics.

Findings

Achieved 99.0% species identification accuracy with clinical samples.

Successfully distinguished MRSA from MSSA and genetically identical strains.

Demonstrated potential for rapid diagnostics in blood, urine, and sputum samples.

Abstract

Rapid identification of bacteria is essential to prevent the spread of infectious disease, help combat antimicrobial resistance, and improve patient outcomes. Raman optical spectroscopy promises to combine bacterial detection, identification, and antibiotic susceptibility testing in a single step. However, achieving clinically relevant speeds and accuracies remains challenging due to the weak Raman signal from bacterial cells and the large number of bacterial species and phenotypes. By amassing the largest known dataset of bacterial Raman spectra, we are able to apply state-of-the-art deep learning approaches to identify 30 of the most common bacterial pathogens from noisy Raman spectra, achieving antibiotic treatment identification accuracies of 99.0$\pm$0.1%. This novel approach distinguishes between methicillin-resistant and -susceptible isolates of Staphylococcus aureus (MRSA and MSSA) as well as a pair of isogenic MRSA and MSSA that are genetically identical apart from deletion of the mecA resistance gene, indicating the potential for culture-free detection of antibiotic resistance. Results from initial clinical validation are promising: using just 10 bacterial spectra from each of 25 isolates, we achieve 99.0$\pm$1.9% species identification accuracy. Our combined Raman-deep learning system represents an important proof-of-concept for rapid, culture-free identification of bacterial isolates and antibiotic resistance and could be readily extended for diagnostics on blood, urine, and sputum.

Equations

$$\sigma =\sqrt{\frac{(n_1-1){\sigma }_1^2+(n_2-1){\sigma }_2^2}{n_1+n_2-2}}.$$

$$\text{se}=\sigma \times \sqrt{\frac{1}{n_1}+\frac{1}{n_2}}.$$

$$t=\frac{{\mu }_1-{\mu }_2}{\text{se}},$$

Full text

Introduction

Bacterial infections are a leading cause of death in both developed and developing nations, taking more than 6.7 million lives each year[1, 2]. These infections are also costly to treat, accounting for 8.7% of annual healthcare spending, or $33 billion, in the United States alone[3]. Current diagnostic methods require sample culturing to detect and identify the bacteria and its antibiotic susceptibility, a slow process that can take days even in state-of-the-art labs[4, 5]. Broad spectrum antibiotics are often prescribed while waiting for culture results[6], and according to the Centers for Disease Control and Prevention, over 30% of patients are treated unnecessarily[7]. New methods for rapid, culture-free diagnosis of bacterial infections are needed to enable earlier prescription of targeted antibiotics and help mitigate antimicrobial resistance.

Raman spectroscopy has the potential to identify the species and antibiotic resistance of bacteria, and when combined with confocal spectroscopy, can interrogate individual bacterial cells (Figure 1a, b). Different bacterial phenotypes are characterized by unique molecular compositions, leading to subtle differences in their corresponding Raman spectra. However, because Raman scattering efficiency is low ($\sim 10^{-8}$ scattering probability[8]), these subtle spectral differences are easily masked by background noise. High signal-to-noise ratios (SNRs) are thus needed to reach high identification accuracies[9], typically requiring long measurement times that prohibit high-throughput single-cell techniques. Additionally, the large number of clinically relevant species, strains, and antibiotic resistance patterns require comprehensive datasets that are not gathered in studies that focus on differentiating between species[10, 11], isolates (typically referred to as strains in the literature)[12, 13], or antibiotic susceptibilities[14, 15, 16, 17, 18, 19]. In this work, we address this challenge by training a convolutional neural network (CNN) to classify noisy bacterial spectra by isolate, empiric treatment, and antibiotic resistance.

Results

Deep learning for bacterial classification from Raman spectra

In order to gather a training dataset, we measure Raman spectra using short measurement times on dried monolayer samples, as illustrated in Figure 1. We ensure that the majority of individual spectra are taken over single cells and preparation conditions are consistent between samples (See Methods). We construct reference datasets of 60,000 spectra from 30 bacterial and yeast isolates for 3 measurement times — these 30 isolate classes cover over 94% of all bacterial infections treated at Stanford Hospital in the years 2016-17 and are representative of the majority of infections in intensive care units worldwide[20]. We further augment our reference dataset with 12,000 spectra from clinical patient isolates, including MRSA and MSSA isolates (see Methods for full dataset information). Previously, the lack of large datasets prohibited the use of CNNs due to the high number of spectra per bacterial class needed for training.

In recent years, CNNs have been applied with tremendous success to a broad range of computer vision problems[21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. However, while classical machine learning techniques have been applied to spectral data[11, 12, 14, 31, 32], relatively little work has been done in adapting deep learning models to spectral data[33, 34, 35, 36]. In particular, state-of-the-art CNN techniques from image classification such as residual connections have previously not been applied to low SNR, 1D spectral data. Our CNN architecture consists of 25 1D convolutional layers and residual connections[37] — instead of two-dimensional images, it takes one-dimensional spectra as input (see Methods for further detail). Unlike previous work, we do not use pooling layers and instead use strided convolutions with the goal of preserving the exact locations of spectral peaks[38]. Empirically, we find that this strategy improves model performance.

We train the neural network on a 30-class isolate identification task, where the CNN outputs a probability distribution across the 30 reference isolates and the maximum is taken as the predicted class. The model is trained on the reference dataset and tested on an independent test dataset gathered from separately cultured samples.

A performance breakdown for individual classes is displayed in the confusion matrix in Figure 2a. Here, we show data for 1 s measurement times, corresponding to a SNR of 4.1 — roughly an order of magnitude lower than typical reported bacterial spectra[11, 12, 10]; classification accuracies increase with SNR, as shown in Supplementary Figure 1. On the 30-class task, the average isolate-level accuracy is $82.2\pm 0.3$% ($\pm$ calculated as standard deviation across 5 train and validation splits). Gram-negative bacteria are primarily misclassified as other Gram-negative bacteria; the same is generally true for Gram-positive bacteria, where additionally, the majority of misclassifications occur within the same genus. In comparison, our implementations of the more common classification techniques of logistic regression and support vector machine (SVM) achieve accuracies of 75.7% and 74.9%, respectively.

Identification of empiric treatments and antibiotic resistance

Species-level classification accuracy is the standard metric for bacterial identification, but in practice, the priority for physicians is choosing the correct antibiotic to treat a patient. Common antibiotics often have activity against multiple species, so the 30 isolates can be arranged into groupings based on the recommended empiric treatment if the bacterial species is known. Classification accuracies can thus be condensed into a new confusion matrix grouped by empiric antibiotic treatment (Figure 2b), where the average accuracy of our method is 97.0$\pm$0.3%. In comparison, logistic regression and SVM achieve accuracies of 93.3% and 92.2%, respectively.

Beyond empiric first choice antibiotics, clinicians also conduct antibiotic susceptibility tests to determine bacterial responses to drugs. As a step toward a culture-free antibiotic susceptibility test using Raman spectroscopy, we train a binary CNN classifier to differentiate between methicillin-resistant and -susceptible isolates of S. aureus. This model achieves 89.1$\pm$0.1% identification accuracy (Figure 3a). Because the consequences for misdiagnosing MRSA as MSSA are often more severe than the reverse misdiagnosis, the binary decision can be tuned for higher sensitivity (low false negative rate), as shown in the receiver operating characteristic (ROC) curve in Figure 3b (dotted line denotes performance of random guessing). The area under the curve (AUC) is 0.953, meaning that a randomly selected positive example (i.e., Raman sample from patient with MRSA) will be predicted to be more likely to be MRSA than a randomly selected negative example (i.e., sample from patient with MSSA) with probability 0.953.

Extension to clinical patient isolates

To demonstrate that this approach can be extended to new clinical settings, we test our model on two groups of 25 clinical isolates derived from patient samples, for a total of 50 patients, Within each patient group, samples include 5 isolates from each of the 5 most prevalent[39] empiric treatment groups (see Supplementary Table 2 and Supplementary Figure 4). We first consider isolates from 25 patients collected from Palo Alto VA Medical Center in 2018. We augment our reference dataset with this clinical dataset comprised of 400 spectra per clinical isolate. To account for changes in the relative prevalence of species and antibiotic resistances over time, the model may be fine-tuned on a small dataset that is representative of current patient populations. We use a leave-one-patient-out cross-validation (LOOCV) strategy for fine-tuning, where we assign 1 patient in each class to the test set (5 patients total) and use the other 4 for fine-tuning (20 patients total), fine-tuning on 10 randomly sampled spectra per patient isolate — we repeat this process 5 times, so all 25 patient isolates appear in the held-out test set once. We then use 10 randomly sampled spectra from each patient isolate in the test set to reach an infection identification for that patient isolate. The sampling procedure for identification is repeated for 10,000 trials, and we report the average accuracy and standard deviation, and display a trial representing the modal result in Figure 4a (full experiment details can be seen in Supplementary Note 1). A CNN pre-trained on the reference dataset serves both as initialization for the fine-tuned model and as a baseline, achieving 89.0$\pm$3.6% ($\pm$ calculated as standard deviation across 10,000 sampling trials) species identification accuracy, a statistically significant improvement over logistic regression and support vector machine baselines (see Methods for details). When the CNN is fine-tuned on clinical data and then evaluated on the held-out patients, the identification accuracy is improved to 99.0$\pm$1.9% (Supplementary Figure 5). Samples for the clinical tests were prepared separately for each patient, so we conclude that the measured performance is not due to batch effects from sample preparation or measurement conditions.

Because patient samples may contain very low numbers of bacterial cells without culturing (e.g. 1 CFU/mL or fewer in blood[40]), only a few individual bacterial spectra per patient may be available to make a diagnosis. As seen in Figure 4c, just 10 cellular spectra are enough to reach high identification accuracy. The rate of correct identification using 10 spectra is 99.0%, within 1% of the performance with 400 spectra (100.0%). While acquiring spectra from 400 individual bacterial cells would likely necessitate culturing, we achieve high accuracy on spectra from 10 individual bacterial cells, commensurate with typical levels of bacterial cells present in uncultured samples[40, 41].

For a proof-of-concept antibiotic susceptibility test on clinical isolates, we collect Raman spectra on 5 additional clinical MRSA isolates and test the binary MRSA/MSSA classifier that is pre-trained on the reference MRSA and MSSA isolates. Using the same LOOCV process, we fine-tune the binary classifier on the clinical spectra. A representative result is shown in Figure 4b; any misclassifications of MSSA as MRSA are labeled as “suboptimal”, indicating that Vancomycin (prescribed for MRSA) is also effective on MSSA but is not considered optimal treatment and may introduce adverse patient effects. On average, the pre-trained binary classifier achieves 61.7$\pm$7.3% accuracy and the fine-tuned binary classifier achieves 65.4$\pm$6.3% accuracy (Supplementary Figure 5).

Finally, to test the robustness of the fine-tuning approach over multiple clinical datasets, we use our second patient group of 25 isolates, collected from Stanford Hospital from February 2019 to March 2019. We conduct additional fine-tuning of the model that is pre-trained on the reference dataset and fine-tuned on the original clinical dataset. The treatment group identification accuracy on the new clinical dataset using only 10 spectra per patient is 99.7$\pm$1.1% Figure 4 d, e, with improved performance for both S. aureus and P. aeruginosa, demonstrating the potential for continuous improvement of the trained model.

Discussion

In this work, we apply state-of-the-art deep learning techniques to noisy Raman spectra to identify clinically relevant bacteria and their empiric treatment. A CNN model pre-trained on our dataset can easily be extended to new clinical settings through fine-tuning on a small number of clinical isolates, as we have shown on our clinical dataset. We envision that fine-tuning processes such as the one demonstrated here could be important components for continuously evaluating and improving deployed models. Our model, applied here to the identification of clinically relevant bacteria, can be applied with minimal modification to other identification problems such as materials identification, or other spectroscopic techniques such as nuclear magnetic resonance, infrared, or mass spectrometry.

This study uses measurement times of 1 s, corresponding to SNRs that are an order of magnitude lower than typical reported bacterial spectra — while still achieving comparable or improved identification accuracy on more isolate classes than typical Raman bacterial identification studies. A common strategy for reducing measurement times is surface-enhanced Raman scattering (SERS) using plasmonic structures, which can increase the signal strength by several orders of magnitude[11, 42, 43]. SERS spectra can be highly variable and difficult to reproduce, particularly on cell samples[44, 8], making it difficult to develop a reliable diagnostic method based on SERS. However, with a dataset capturing the breadth of variation in SERS spectra, a CNN could enable a platform that processes blood, sputum, or urine samples in a few hours.

Compared to other culture-free methods[45] including single-cell sequencing[46, 47, 48, 49] and fluorescence or magnetic tagging[50], Raman spectroscopy has the unique potential to be a technique for identifying phenotypes that does not require specially designed labels, allowing for easy generalizability to new strains.

To achieve treatment recommendations as fine-grained as those from culture-based methods, larger datasets covering more resistant and susceptible clinical isolates, greater diversity in antibiotic susceptibility profiles, cell states, and growth media and conditions would be needed. Though collecting such datasets is beyond an academic scope, requiring highly automated sample preparation and data acquisition processes, there is promise for clinical translation. Similarly, studies applying the Raman-CNN system to identify pathogens in relevant biofluids such as whole blood, sputum, and urine are a promising future direction to demonstrate the validity of the method as a diagnostic tool. When combined with such an automated system, the Raman-CNN platform presented here could rapidly scan and identify every cell in a patient sample and recommended an antibiotic treatment in one step, without needing to wait for a culture step. Such a technique would allow for accurate and targeted treatment of bacterial infections within hours, reducing healthcare costs and antibiotics misuse, limiting antimicrobial resistance, and improving patient outcomes.

Methods

Dataset

The reference dataset consists of 30 bacterial and yeast isolates, including multiple isolates of Gram-negative and Gram-positive bacteria, as well as Candida species. We also include an isogenic pair of S. aureus from the same strain, in which one variant contains the mecA resistance gene for methicillin (MRSA) and the other does not (MSSA)[51] (see Supplementary Table 1 for full isolate information). The reference training dataset consists of 2000 spectra each for the 30 reference isolates plus isogenic MSSA at 3 measurement times. The reference fine-tuning and test datasets each consist of 100 spectra for each of the 30 reference isolates. The first clinical dataset consists of 30 patient isolates distributed across 5 species, with 400 spectra per isolate. The second clinical dataset consists of 25 patient isolates distributed across the same 5 species, with 100 spectra per isolate. Due to degradation in optical system efficiency, the measurement times for the reference fine-tuning and test and second clinical datasets were increased from 1 s to 2 s in order to keep SNR consistent across datasets. Antibiotic susceptibility was performed by first genotypic testing for methicillin by detecting mecA using PCR (PMID: 19741081). Then phenotypic antimicrobial susceptibility testing was performed on the Microscan Walkaway instrument (Beckman Coulter, Brea, CA) and VITEK️ 2 (Biomérieux, Inc., Durham, NC).

Dataset variance

For our datasets, we observe that intra-sample variance is high, as demonstrated by the pairwise spectral difference analysis summarized in Supplementary Figure 2. For 19 out of 30 isolates, spectra from at least one other isolate are more similar on average than spectra from the same isolate, on average. For example, when we rank isolates in order of similarity to E. faecalis 2 (Supplementary Figure 2c), there are 8 other isolates where the average difference between a spectrum from E. faecalis 2 and a spectrum from the other isolate is smaller than the average difference between two spectra from E. faecalis 2. When intra-sample variance is high, a large number of spectra per sample may help to better represent the full data distribution and lead to higher predictive performance.

Sample preparation

Bacterial isolates were cultured on blood agar plates each day before measurement. Plates were sealed with Parafilm and stored at 4°C for 20 minutes to 12 hours before sample preparation. Storage times varied to allow for multiple measurement times per day; however all other sample preparation conditions were kept consistent between samples. Differences in storage time were not found to result in spectral changes greater than spectral changes due to strain or isogenic differences. All clinical isolates were prepared in separate samples with consistent sample preparation conditions. Because test samples were prepared separately from samples used for training, we conclude that classifications are not due to batch effects such as differences in sample preparation. We prepared samples for measurement by suspending 0.6 mg of biomass from a single colony in 10 $\mu$L of sterile water (0.4 mg in 5 $\mu$L water for Gram-positive species) and drying 3 $\mu$L of the suspension on a gold-coated silica substrate (Figure 1a and b). Substrates were prepared by electron beam evaporation of 200 nm of gold onto microscope slides that were pre-cleaned using base piranha. Samples were allowed to dry for 1 hour before measurement.

Raman measurements

We measured Raman spectra across monolayer regions of the dried samples (Figure 1a) using the mapping mode of a Horiba LabRAM HR Evolution Raman microscope. 633 nm illumination at 13.17 mW was used with a 300 l/mm grating to generate spectra with 1.2 cm*-1* dispersion to maximize signal strength while minimizing background signal from autofluorescence. Wavenumber calibration was performed using a silicon sample. The 100X 0.9 NA objective lens (Olympus MPLAN) generates a diffraction-limited spot size, $\sim$1 $\mu$m in diameter. A 45x45 discrete spot map is taken with 3 $\mu$m spacing between spots to avoid overlap between spectra. The spectra are individually background corrected using a polynomial fit of order 5 using the subbackmod Matlab function available in the Biodata toolbox (see Supplementary Figure 1 for examples of raw and corrected spectra). The majority of spectra are measured on true monolayers and arise from 1̃ cell due to the diffraction-limited laser spot size, which is roughly the size of a bacteria cell. However, a small number of spectra may be taken over aggregates or multilayer regions. We exclude the spectra that are most likely to be non-monolayer measurements by ranking the spectra by signal intensity and discarding the 25 spectra with highest intensity, which includes all spectra with intensities greater than two standard deviations from the mean. We measured both monolayers and single cells, and found that monolayer measurements have SNRs of 2.5$\pm$0.7, similar to single-cell measurements (2.4$\pm$0.6), while allowing for the semi-automated generation of a large training dataset. The spectral range between 381.98 and 1792.4 cm*-1* was used, and spectra were individually normalized to run from a minimum intensity of 0 and maximum intensity of 1 within this spectral range. SNR values are calculated by dividing the total intensity range by the intensity range over a 20-pixel wide window in a region where there is no Raman signal.

CNN architecture & training details

The CNN architecture is adapted from the Resnet architecture[37] that has been widely successful across a range of computer vision tasks. It consists of an initial convolution layer followed by 6 residual layers and a final fully connected classification layer — a block diagram can be seen in Figure 1. The residual layers contain shortcut connections between the input and output of each residual block, allowing for better gradient propagation and stable training (refer to reference 37 for details). Each residual layer contains 4 convolutional layers, so the total depth of the network is 26 layers. The initial convolution layer has 64 convolutional filters, while each of the hidden layers has 100 filters. These architecture hyperparameters were selected via grid search using one training and validation split on the isolate classification task. We also experimented with simple MLP (multi-layer perceptron) and CNN architectures but found that the Resnet-based architecture performed best.

We first train the network on the 30-isolate classification task, where the output of the CNN is a vector of probabilities across the 30 classes and the maximum probability is taken as the predicted class. The binary MRSA/MSSA and binary isogenic MRSA/MSSA classifiers have the same architecture as the 30-isolate classifier, aside from the number of classes in the final classification layer. We use the Adam optimizer[52] across all experiments with learning rate 0.001, betas (0.5, 0.999), and batch size 10. Classification accuracies are reported across 5 randomly selected train and validation splits. We first pre-train the CNN on the reference training dataset, then fine-tune on the reference fine-tuning dataset to account for measurement changes due to degradation in optical system efficiency. For each of the 5 splits, we split the fine-tuning data into 90/10 train and validation splits, train the CNN on the train split, and use the accuracy on the validation split to perform model selection. We then evaluate and report the test accuracy on the test dataset which is gathered from independently cultured and prepared samples. The binary MRSA/MSSA classifier is trained and fine-tuned using the same procedure. The binary isogenic MRSA/MSSA classifier is trained using a similar procedure on data from a single measurement series.

All error values reported for tests on the reference dataset are standard deviation values across 5 splits.

While a high number of samples is good for ensuring dataset variation, deep learning approaches can still benefit from having a high number of examples per sample. When intra-sample variance is high, as we observe for our datasets, a large number of spectra per sample may better represent the full distribution and lead to higher predictive performance.

For the clinical isolates, we start by pre-training a CNN on the empiric treatment labels for the 30 reference isolates. We then use the following leave-one-patient-out cross-validation (LOOCV) strategy to fine-tune the parameters of the CNN. There are a total of 25 patient isolates across 5 species. In each of the 5 folds, we assign 1 patient in each species to the test set, 1 patient in each species to the validation set, and the remaining 3 patients in each species to the training (i.e., fine-tuning) set. We then use the clinical training set (consisting of isolates from 15 patients) to fine-tune the CNN parameters, and use accuracy on the validation set (5 patient isolates) to do model selection. The test accuracy for each fold is evaluated on the test set (5 patient isolates) using the method described below.

Clinical identification data analysis

To reach an identification for patient isolates, 400 spectra are measured across a sample from each patient isolate. 10 of these spectra are chosen at random to be classified. The most common class out of the 10 spectral classifications is then chosen as the identification for each patient isolate, with ties broken randomly. All error values reported for tests on the clinical dataset are standard deviations across 10,000 trials of random selections of 10 spectra, with an upper accuracy bound of 100%. For the second clinical dataset, we perform the same procedure, except that we choose 10 out of 100 spectra for each patient isolate, and use a model that is both pre-trained on the reference dataset and fine-tuned on the first clinical dataset.

Baselines

In all experiments where logistic regression (LR) and support vector machine (SVM) baselines were used, we first used PCA to reduce the input dimension from 1000 to 20 — this hyperparameter was determined by plotting test accuracies for different settings on one training and validation split for the 30 isolate task and picking a value near where the test accuracy saturated. Using only the first 20 principal components not only decreases computation costs, but also increases accuracy by reducing the amount of noise in the data. For each fold of the cross validation procedure, we use grid search to choose the regularization hyperparameter for each model achieving the best validation accuracy and report the corresponding test accuracy. Using both the training and fine-tuning reference datasets to train the baseline models, LR and SVM achieve 57.5% and 56.8% on the 30-class task and 89.0% and 88.3% on the empiric treatment task, respectively. Using only the fine-tuning reference dataset, LR and SVM achieve 75.7% and 74.9% on the 30-class task and 93.3% and 92.2% on the empiric treatment task, respectively. The latter performance is higher because the baseline models do not benefit from additional training data as the CNN does, but rather benefit from training data the most closely matches the measurement conditions of the test data.

Two-sample test of sample means

We use the Welch’s two-sample $t$-test to test whether the differences in mean clinical accuracy for the CNN and the SVM and LR baselines were statistically significant. Welch’s $t$-test is a variation of the Student’s $t$-test that is used when the two samples may have unequal variances. In each case, we start by computing the pooled standard deviation as

[TABLE]

We then compute the standard error of the difference between the means as

[TABLE]

Finally, we can compute the test statistic as

[TABLE]

and then compute the p-value using the corresponding Student’s $t$-distribution. For our computations, $n_{\text{CNN}}=n_{\text{LR}}=n_{\text{SVM}}=10000$, $\mu_{\text{CNN}}=89.0$, $\mu_{\text{LR}}=81.8$, $\mu_{\text{SVM}}=82.9$, $\sigma_{\text{CNN}}=3.6$, $\sigma_{\text{LR}}=6.0$, and $\sigma_{\text{SVM}}=5.9$. In comparing the CNN with LR, we computed a $t$-statistic of 102.9 and in comparing the CNN with SVM, we computed a $t$-statistic of 88.3. In both cases, we reject the null hypothesis that the means are equal at the 1e-6 p-level.

Data availability

All data needed to replicate these results are available at https://github.com/csho33/bacteria-ID.

Code availability

All code needed to replicate these results is available at https://github.com/csho33/bacteria-ID.

Biological materials availability

Unique isolates are available from the authors upon reasonable request.