The following bacteria were used in the study; they are Salmonella typhimurium (S. typhimurium) (LT2, Sa 11030), Shigella flexneri (S. flexneri), Listeria monocytogenes (L. monocytogenes) (Lin), Vibrio cholerae (V. cholerae) (Non-toxigenic strain, 93097), Staphylococcus aureus (S. aureus) (ATCC 25923), and C. botulinum. The strains were provided by the State Key Laboratory for Infectious Disease Prevention and Control, National Institute for Communicable Diseases Control and Prevention, and Chinese Center for Disease Control and Prevention (ICDC; Beijing, China).

C. botulinum strains were stored in TPGY broth mixed with glycerol at –70^°C; other strains were housed in Luria-Bertani broth mixed with glycerol at –70^°C. To multiply bacterial cells, C. botulinum strains were grown in TPGY broth for 24 h at 37^°C anaerobically; other strains were grown in Luria-Bertani broth for 24 h at 37^°C. The culture media were purchased from Beijing Land Bridge Technology Co., Ltd. (Beijing, China). All spectra were collected during the stationary phase to avoid the influence of different growth periods of bacteria. In this study, their initial and final OD₆₀₀ values are approximately equal based on the results of ultraviolet spectrophotometer (Varian, USA).

Each culture was vortexed quickly and centrifuged for 8 min at 4,000 rpm at 4^°C using Centrifuge 5418 R (Eppendorf, Germany). The supernatant was discarded after centrifugation. Subsequently, 5 mL of 0.9% NaCl solution was added and stirred. This procedure was repeated one time. The sediment was suspended in 1 mL of 0.9% NaCl solution and transferred to a new microcentrifuge tube (1.5 mL). The suspension was centrifuged for 3 min at 13,500 rpm at room temperature using microcentrifuge (Pic017; Heraeus, Germany). This procedure was repeated two times. Finally, the sediment was resuspended in 200 μL of 0.9% NaCl solution. To record the spectra, about 5 μL of solution was dropped onto the aluminum foil and dried at room temperature for 5 min for Raman spectra measuring (Zhang et al., 2021).

All Raman experiments were performed using confocal and high performance Raman microscope (XploRA PLUS, HORIBA, Japan). Bacteria samples were detected with a 532-nm laser under a 600 g/mm grating, the laser power is 1 mW, a spectral resolution of 0.6 cm^–1, and 50X objective lens. The integration time was 10 s. Each bacteria sample was collected six times. Each bacterial spectrum was collected and consisted of 603 points in the range of 600–1,800 cm^–1. All the spectra were processed via Origin software (OriginLab, USA). The background was removed using baseline correction (method: second derivative, baseline mode: user defined, number of baseline points: 16). LabSpec 6.3 (HORIBA, Japan) was utilized to optimize acquisition parameters and collect sufficient Raman spectra for statistical analysis. A single bacteria cell was mapped since the advantage of confocal Raman. Eighty samples of each type of bacteria were cultured, Raman spectra were collected, respectively, and a total of 480 Raman spectra were collected.

By increasing the integration time and the number of scans, the spectra jumper by external interference can be eliminated, and the signal strength of the Raman spectrum can be promoted. However, this method takes longer to collect spectra. Therefore, it is necessary to choose an appropriate integration time and number of scans. Clostridium botulinum was chosen as the test object. Based on the same experimental conditions, the integration time was set to 5, 10, and 15 s, respectively, and the Raman spectra of Clostridium botulinum under different integration times were compared. Based on the same experimental conditions, the number of scans was set to 3, 6, and 9, respectively, and the Raman spectra of Clostridium botulinum under different scan times were compared.

In each bacterium, a sample was randomly selected and tested 10 times in parallel to examine the repeatability of the method.

The dataset was preprocessed by four single pretreatment methods, respectively, to explore the impact of various methods on the model. The pretreatment methods included Savitzky–Golay algorithm smoothing (SG), SNV, MSC, and SG 1st Der. Among them, SG can remove the noise of the spectral. SNV and MSC can reduce the influence of scattering on the original spectrum. SG 1st Der can eliminate the interference of the baseline and background. All the preprocessing procedures were carried out in the Unscrambler X 10.4 (CAMO, Norway).

Principal component analysis (PCA) is an unsupervised and multivariate technique that projects a set of correlated features onto a set of uncorrelated features using an orthogonal transformation. It preserves the most of information. In this study, PCA was used for two purposes. The first one is to plot the distribution of data based on the PCA scores to investigate the difference of various samples. The second one is to reduce the dimensionality of the dataset and extract feature. Computational costs of classifiers were reduced, and overfitting was prevented.

Decision tree (DT) gives the various outcomes from a series of decisions based on a flowchart-like diagram. It is a strong tool for planning strategy, research analysis, and making decision. One of primary advantage is easy to follow and understand. There are four popular algorithms in DT model; they are CHAID, exhaustive CHAID, CART, and QUEST algorithm. All of them were applied to build DT models.

ANN is inspired based on modern neuroscience research. The large amount of processing unit is utilized to build a complex model. Human brain neural network structure and function are imitated in the ANN model. ANN is widely applied for spectral analysis and identification (Park and Lek, 2016). There are two popular algorithms in ANN model; they are multilayer perceptron (MLP) and radial basis function (RBF). All of them were applied to build models.

FDA is a supervised technique using a discriminant function to assign data to different groups. FDA is often used to build a model combined with PCA. The PCs from PCA are used in FDA to define and predict classes.

The hold-out method is used as a method of cross-validation in DT, ANN model. In the hold-out method, the dataset is divided into two parts, 70% of the samples are allocated to the training set, and there are 320 samples in the training set. About 30% of the samples are assigned to the testing set, and there are 160 samples in the testing set. They were used to test the predictions of the models (Liu et al., 2018). Leave-one-out cross-validation (LOOCV) technique is used for the cross-validation of FDA model, and it can examine the quality of the classifier and avoid overfitting. In the LOOCV, one spectrum in the dataset is allocated to the validation data and the remaining spectra are allocated to the training data. This process is repeated until each spectrum is used once as the validation data (Tie et al., 2020). Therefore, there are 479 samples in the training set, and there is a sample in the testing set. The tests of equality of group means are used for evaluating the potential of every independent variable before constructing the model. CRR, area under the ROC, cumulative gains, and lift charts were used for comparing the performance of classifiers under different pretreatment methods. All statistic procedures were operated in the IBM SPSS Statistics 25 and the Unscrambler X 10.4 software.

In this study, the CRM technique was used for the detection and identification of six in total foodborne pathogenic bacteria species, namely S. typhimurium, S. flexneri, L. monocytogenes, V. cholerae, S. aureus, and C. botulinum.

The mean normalized CRM spectra of six bacterial species are shown in Figure 3. Major spectral bands and peaks assignment are shown in Table 1. Some spectral bands are common in all species, and every spectrum showed bands at ca. 746, 1,140, and 1,245 cm^–1 assignable to cytosine and uracil, = C-O-C = (unsaturated fatty acids in lipids), amide III (random), and thymine, respectively. There are also many differences in the spectral images. For example, only L. monocytogenes showed a band at ∼1,094 and 1,574 cm^–1 assignable to CC skeletal and COC stretch (Zheng et al., 2020) from glycosidic link, CN stretching of amide II (Fan et al., 2011), respectively; only S. aureus showed a band at ∼1,302 cm^–1 assignable to amide III (Jarvis et al., 2004). A band was located at 1,658 cm^–1 in all species except V. cholerae. There are few differences among C. botulinum and S. flexneri. However, it is possible to differentiate them by using the ratio of the peak intensities. A total of 1,322/1,350 and 1,463/1,587 were used for differences between C. botulinum and S. flexneri. A significant difference of peak intensity ratios was shown as depicted in Table 2. Therefore, the different species can be characterized based on the abundant and unique spectral information of CRM.

As shown in Figure 3, mean normalized spectra of different bacterial species have a high level of similarity; although they can be differentiated based on several peaks and peak intensity ratio, visual detection of these minor differences is quite time-consuming and may lead to misdiagnosis (Lin et al., 2019). Thus, it is essential to analyze the data combined with multivariate analysis techniques for revealing minor spectral differences. In this study, PCA, DT, ANN, and FDA were used for further study.

In this study, the PCA method was used to construct classification models of bacterial species based on the four kinds of pretreatment methods, and it has reduced the dimensionality of the dataset and extracted feature.

Table 3 shows the results of PCA under different pretreatment methods. It can be seen that the number of PCs and corresponding cumulative contribution rate has changed along with difference of pretreatment methods. Figures 4A–E represents PCA score plots under various preprocess methods. These results represented the influence of preprocessed spectra on PCA models and the possibility of differentiation for bacterial species with a high accuracy. Figure 4A shows that all bacterial species can be differentiated. In Figures 4B,D, the samples of S. aureus and S. typhimurium became more scattered, the categories were closer to each other, and it was difficult to distinguish different categories. In Figure 4C, although the spatial distance of each group became larger, the problem of overlap was more serious. In Figure 4E, the plot had placed samples of the same types closer by SG 1st Der, and overlapping samples were reduced between C. botulinum and S. flexneri compared to Figure 4C.

Figure 4F represents the most important variables for classification in the original data. The major positive loadings were located at ca. 1,322, 1,353, 1,465, and 1,590 cm^–1. The largest loadings in the positive direction were located at ca. 1,353 cm^–1 (Amide III (protein)) and ca. 1,590 cm^–1 (Guanine ring stretching). The major negative loadings were located at ca. 755, 1,001, 1,123, 1,221, and 1,302 cm^–1. The largest loadings in the negative direction were located at ca. 755 cm^–1 (Cytosine or uracil) and ca. 1,302 cm^–1 (Amide III).

As mentioned above, we can conclude that preliminary classification can be achieved by PCA model combined with the preprocess methods. It is worth noting that PCA is an unsupervised technique, not meant for classification purpose. DT, ANN, and FDA, these supervised and multivariate techniques, were used for the construction of classifiers. They have more advantages than PCA because the information among the classes is taken into account (Wichmann et al., 2020).

To evaluate the potential of every independent variable before constructing the classification model, the tests of equality of group means are used for this study. Supplementary Tables 1–5 give the summary of the tests of equality of group means with spectra preprocessed by various pretreatment methods. In each test, the grouping variable was treated as the factor, and the results of a one-way ANOVA are given for the independent variable. Smaller Wilk’s lambda values suggest that the variable is better at discrimination of groups. If the significance value (Sig) of the independent variable is less than 0.05, it indicates that the variable is significant for the creation of the models. As shown in Supplementary Tables 1–5, Wilk’s lambda value of major PCs is smaller, and every Sig value is far less than 0.5. It indicates that every variable is significant for the construction of the models and there is not any evident anomaly (Chiang et al., 2001).