{'fraction_non_alphanumeric': 0.041188598063822156, 'fraction_numerical': 0.016851918250268914, 'mean_word_length': 4.74633015709503, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nondestructive microbial discrimination using single-cell Raman spectra and random forest machine learning algorithm

Summary Raman microspectroscopy is a powerful tool for obtaining biomolecular information from single microbial cells in a nondestructive manner. Here, we detail steps to discriminate prokaryotic species using single-cell Raman spectra acquisitions followed by data preprocessing and random forest model tuning. In addition, we describe the steps required to evaluate the model. This protocol requires minimal preprocessing of Raman spectral data, making it accessible to non-spectroscopists, yet allows intuitive visualization of feature importance. For complete details on the use and execution of this protocol, please refer to Kanno et al. (2021).


SUMMARY
Raman microspectroscopy is a powerful tool for obtaining biomolecular information from single microbial cells in a nondestructive manner. Here, we detail steps to discriminate prokaryotic species using single-cell Raman spectra acquisitions followed by data preprocessing and random forest model tuning. In addition, we describe the steps required to evaluate the model. This protocol requires minimal preprocessing of Raman spectral data, making it accessible to non-spectroscopists, yet allows intuitive visualization of feature importance. For complete details on the use and execution of this protocol, please refer to Kanno et al. (2021).

BEFORE YOU BEGIN Background
Despite substantial advances in cell culture methods, the vast majority of microorganisms remain uncultured (Rinke et al., 2013). To address this microbial dark matter, culture-independent single-cell approaches have recently been developed and applied to a wide range of microbial systems regardless of whether they are laboratory or environmental samples. Raman microspectroscopy, which is based on inelastic light scattering from molecules, has been shown to be powerful for studying single microbial cells in a nondestructive manner through their molecular fingerprints (i.e., Raman spectra) (Lorenz et al., 2017). Artificial intelligence technology is also being actively used to analyze microbial data, which can often become very large. Combining these two relatively new technologies will enable accurate discrimination of diverse microbial species simply by using Raman spectra of microbial cells (Ho et al., 2019;Lu et al., 2020). We have demonstrated a machine learning-based, nondestructive discrimination of prokaryotic species at the single-cell level. The algorithm employed here, random forest (RF), has the advantage of readily visualizing the importance of the features that contribute to the discrimination. The protocol below describes step-by-step details for confocal Raman microspectroscopic measurements, data preprocessing, and RF model construction involved in microbial discrimination. REAGENT

MATERIALS AND EQUIPMENT
Note: Autoclaved PBS can be stored at 4 C for more than six months.
Note: Adjust the pH and/or salt concentration of PBS if cell lysis occurs. Alternatively, try another buffer such as bicarbonate buffer and HEPES buffer.

STEP-BY-STEP METHOD DETAILS
Set up your confocal Raman microscope

Timing: 1-2 h
This section describes the setting of the confocal Raman microscope you use prior to measurements of single microbial cells. We used a laboratory-built confocal Raman microscope with the 632.8 nm light from a He-Ne laser as the Raman excitation source ( Figure 1) (Huang et al., 2011(Huang et al., , 2012Matsuda et al., 2019).
Alternatives: A commercial Raman microscope available from Renishaw, HORIBA, Bruker, WI-Tec, etc. can be used as well. We recommend an inverted microscope because it facilitates bright-field observation of cell suspension dropped in a glass bottom dish.
Note: The same protocol can be applied for 532 nm excitation. Modifications may, however, be required in data preprocessing (see step 11 below), because 532 nm excitation often causes stronger cellular autofluorescence than 632.8 nm excitation.
CRITICAL: The following parameters or components of the confocal Raman microscope ( Figure 1) are of particular importance because they substantially affect measured Raman spectra: laser wavelength; an objective lens; a pinhole for confocal detection; and a grating of the spectrometer.

Laser wavelength
The choice of the laser wavelength concerns with the resonance Raman effect and laser-induced damage of the cell. 532 nm excitation is most commonly used, but in this protocol, we adopted 632.8 nm excitation, which falls near the optical window of biology.

Objective lens
A phase-contrast objective is generally used to observe prokaryotic cells. We chose an oil-immersion phase-contrast objective (Nikon, 1003, NA 1.3, CFI Plan Fluor DLL). According to the Rayleigh criterion, the lateral (XY) resolution is given by 0.61l/NA, where l is the wavelength of light and NA is the numerical aperture of the objective. Using this relation, the lateral resolution of our Raman microscope is calculated to be 0.3 mm, while experimentally it is estimated to be $0.45 mm. Confocal pinhole A 100-mm pinhole is used to get rid of out-of-focus light; the resulting depth (Z) resolution is estimated to be $4.5 mm.

Grating
The use of a 600 grooves/mm grating enables us to record the entire spectral window (approximately 660-3,020 cm À1 ) covering both fingerprint and CH-stretching regions with a spectral resolution of $5 cm À1 .
1. Wait typically more than 30 min until the laser output and the temperature of the CCD detector (set to À70 C) are stabilized. Re-align the optical path if necessary. 2. Set the laser power at the sample point to a desired value (3 mW in the present case) by adjusting a variable neutral-density (ND) filter.
Note: The desired value of the laser power depends on purposes. If you do not care about the viability of the measured cell, you could use a higher power (e.g., $10 mW (Ho et al., 2019)) and improve the signal-to-noise ratio (S/N) of the spectrum.
3. Set the center wavenumber of the grating in order to determine the wavenumber region to be measured.
Note: The calibrated wavenumber region varies slightly from day to day, resulting in a shift in the horizontal axis of the measured spectra. This shift can be corrected for by using emission lines of a standard neon lamp (model: 6032, Newport) recorded on different days. To make the wavenumber regions as close as possible among different measurements, adjust them so that the peak position of a selected Raman band of the reference material comes at roughly the same CCD pixel.
Note: In principle, do not change the grating setting during the experiment.

Single-cell Raman measurements
Timing: 1.5 h for 40 cells/species This section outlines steps to acquire Raman spectra of single cells of a microbial species. CRITICAL: When the cell density is too high (e.g., 5 3 10 7 cells/mL), multiple cells may be laser-trapped simultaneously (see Methods video S1). Adjust the cell density so that there are not too many ( Figure 2B, right) or too few cells in the field of view. For example, for E. coli cultured for 14-16 h (about 5 3 10 9 cells/mL) in LB medium, add a 0.5 mL aliquot to 200 mL of PBS, resulting in a cell density of $1.3 3 10 7 cells/mL.  (Xie et al., 2002), and record its Raman spectrum (see Methods video S2).
Note: We used laser power of 3 mW at the sample point and an exposure time of 30 s per cell. The total timing depends on the exposure time and how many cells you measure for each species. Troubleshooting 2.
Note: Make sure that the measured spectrum is free from the glass signal ( Figure 3B, red). When a cell very close to the cover glass (see Figure 3A-ii) is trapped and measured, the resulting Raman spectrum ( Figure 3C, lower trace) shows a Raman band of glass at $915 cm À1 , which is absent in the one measured at the proper height from the cover glass ( Figure 3C, upper trace).
6. Record the PBS spectrum ( Figure 3B, blue) 10 times for averaging purposes. This section describes how to preprocess the recorded Raman spectra prior to training and testing in machine learning. The overall workflow is illustrated in Figure 4, which is adapted with permission from (Kanno et al., 2021).
Alternatives: This protocol uses Igor Pro 9.0.0.10 (WaveMetrics) to perform all of the following data preprocessing, but other data analysis software such as Origin (https://www.originlab. com/) can also be used.
7. Remove spiky artifact caused by cosmic ray, if any. 8. Subtract the PBS spectrum (average of 10 spectra) from each spike-free, single-cell Raman spectrum ( Figure 4A). 9. Use emission lines of a standard neon lamp to correct for a daily shift of the horizontal axis (Figure 4B). a. For each measurement day, fit a particular emission line to the Lorentzian function and determine the wavenumber of its peak ( Figure 5). b. Delete (typically a few) data points at both ends of each spectrum so that the peak wavenumbers of the above neon emission line match that of a reference spectrum within G0.5 pixel that is pre-determined from all of the data from different days. 10. Delete the so-called silent region (approximately 1,800-2,775 cm À1 ) and retain the fingerprint (660-1,800 cm À1 ) and CH-stretching (2,775-3,020 cm À1 ) regions ( Figure 4C).
Note: No essential Raman bands appear in the silent region with few exceptions (e.g., ChC and ChN stretching bands (Yamakoshi et al., 2011;Zhao et al., 2017) and overtone and combination bands (Horiue et al., 2020)). CRITICAL: What function is used to fit the baseline could significantly affect the performance of machine learning. Other functions than a linear function are available in this package. However, using a more complex function such as higher-order polynomials may often lead to overfitting. To avoid overfitting, linear approximation is a simple and effective means, although it may not be the best.
c. Drag on the graph and then click the '+/-' button to select/deselect the region(s) to be fitted.
A linear fit will immediately be shown as a blue line ( Figure 6A).
Note: Fit only the higher and lower wavenumber edges of each region, as indicated by grey parts in Figure 4D and light-blue parts in the Graph in Figure 6A.
d. Click 'Subtract' for a selected spectrum and 'All in One' for all spectra displayed on the Graph, yielding baseline-subtracted Raman spectra ( Figure 6B). 12. Perform vector normalization on each region ( Figure 4E), by dividing the intensity at each pixel by the square root of the sum of the squared intensities of the spectrum. Troubleshooting 3. 13. Combine the normalized spectra in the fingerprint and CH-stretching regions to obtain a preprocessed single-cell Raman spectrum ( Figure 4F).

Tuning hyper-parameters for random forest model construction
Timing: 1 h This section details the procedures for optimizing the parameters needed for constructing RF machine learning models using the preprocessed data (40 spectra per species) in Python ( Figure 7A). First, a data file to be imported into Python is created. Subsequently, an optimal set of parameters is sought for by using grid search in Python, with a script in the Jupyter Notebook format, which is available on Zenodo (see key resources table). We used Jupyter Notebook in JupyterLab, a webbased interactive development environment.
14. Prepare a data file containing preprocessed single-cell Raman spectra in text format using Microsoft Excel or other software. An example is shown in Figure 8. The columns 'waves' and 'species' specify the names of spectra and class labels (i.e., microbial species), respectively. 15. Import RandomForestClassifier and GridSearchCV from the scikit-learn package into Python (refer to the Jupyter notebook 'GridSearch_10models.ipynb').   Note: N differs depending on the size of the dataset. In the case of 10-fold cross validation (N = 10), which we employed, the dataset is split into 10 folds, and 10 patterns of training (9 folds) and test (1 fold) sets are generated and exported as CSV files. You will find the resulting files as 'train_01.csv', 'test_01.csv', ., 'train_10.csv', and 'test_10.csv' in the working directory. The same number indicates a train/test set.
18. Execute grid search for hyper-parameter optimization with M-fold (e.g., M = 5) cross-validation on each training set.
Note: The hyper-parameters, n_estimators and max_features ( Figure 7B), that are frequently adopted in the N models are used to search for the optimal parameters that achieve high classification accuracy in the test datasets.
a. Select a model used for grid search.
b. Set parameters for the search (see the example script below). The combinations of the numbers given in [ ] of 'n_estimators' and 'max_features' are attempted.
Note: If you set too many parameters to search, it will take an enormous amount of time. We recommend that you first attempt a broad search, then narrow it down to the finest details, and search again. train_concat.to_csv("train_"+"%02.f"%(i)+".csv", index=False) test_concat.to_csv("test_"+"%02.f"%(i)+".csv", index=False)  Note: 'n_estimators' specifies the number of decision trees, while 'max_features' is one of the parameters that restrict the shape of the decision trees and sets an upper limit of the number of features used for node splitting in each decision tree ( Figure 7B).

Note:
In steps 17 and 18b, a random number needs to be specified ('random_state'). This ensures that the same result will be obtained whenever the same training samples and parameters are used.
Optional: In addition to 'max_features', there are other parameters, such as 'max_depth', 'min_samples_split', and 'min_samples_leaf' in the scikit-learn library, that restrict the shape of the tree for reducing overfitting.
Optional: For simplicity, dimensionality reduction (Guo et al., 2021) is not employed in this protocol, which may be useful in some cases for reducing the number of features used for training.

Evaluation of random forest discrimination model
Timing: 1 h 19. Using the optimized hyper-parameters, construct a RF model for each combination of training and test sets and evaluate the constructed model using the test set ( Figure 7A). For details, refer to 'RF_10models.ipynb'. a. Set the hyper-parameters derived from the above grid search. c. For each train/test set, output the feature importance and class probability in the CSV format and the confusion matrix, and calculate the accuracy score. d. Calculate the score of feature importance as the total reduction of the criterion brought by that feature (Gini importance) in each validation. Troubleshooting 4.
Note: An average of the scores over N validations (N = 10 in the present case) is used to represent each feature importance. The means of the feature importance ( Figure 10A) and class probability ( Figure 10B) are obtained from the exported data using Microsoft Excel or similar spreadsheet software.
e. Display the accuracy for each train/test set.
f. Display the mean accuracy and standard deviation of 10 models. Troubleshooting 5.
20. Generate a confusion matrix to visualize the performance of the classification model, in which the sum (or average) of the N patterns of results is displayed in a table layout (Figure 11). a. Import all of the 10 arrays obtained from the above classification. c. Set up displaying of a confusion matrix. The parameters 'normalize=True' and 'normalize=-False' specify that the values in the matrix are displayed as percentage or as they are, respectively. Note: The order of the labels to be included in the 'species' should be the same as that of the 'labels' above (see step 19b).
Note: Recently, a method for visualizing feature importance in convolutional neural network (CNN) has been reported (Lu et al., 2020), but it requires complicated calculations and own coding. In contrast, this protocol using RF can extract feature importance more easily with standard commands available in Python.

OPEN ACCESS
Note: In constructing each decision tree, a subset of the spectral data sampled by bootstrapping (random resampling of data with replacement) is used, thereby leaving unsampled data (i.e., OOB samples). The performance of each tree can be validated by using these OOB samples. As can be seen from Figure 7C, the OOB error rate in general decreases with increasing the number of trees ('n_estimators'). This curve is used to derive an optimal value of 'n_estimators' (Figure 9); it converges for the 'n_estimator' value larger than $200 in the present case. Detailed instructions on how to display the OOB error rate can be found on the scikit-learn library web site: https:// scikit-learn.org/stable/auto_examples/ensemble/plot_ensemble_oob.html?highlight=oob (our script file deposited on Zenodo also contains the code for OOB error rate).
21. Assign the wavenumbers with high feature importance scores obtained in step 19d by referring to one of the many previous papers that show tables of the Raman band assignments for major biological molecules; see, for example, (Carey, 1982;Kuhar et al., 2018;Matsuda et al., 2019;Naumann, 2001).

EXPECTED OUTCOMES
A highly accurate discrimination model for microbial cells can be constructed by using single-cell Raman spectra and the RF algorithm. It yields information on the features (Raman shift) that make significant contributions to discrimination between classes. Notable advantages of this protocol are that it can achieve sufficient discrimination accuracy with less computational load and that it is straightforward to extract the importance of features.
Depending on your research purposes, you can change the class label (e.g., from species to a higher or lower taxonomic level), as well as microbial species to be studied. For example, the present RF model was shown to work very well for bacteria/archaea binary classification, where the class of each sample was changed from species to domain. The preprocessed Raman spectral data (Figure 4F) can serve as input data for other types of analysis, such as multivariate analysis (e.g., principal component analysis and clustering analysis) and comparative analysis of peak intensities.

LIMITATIONS
With this protocol, we did not evaluate discrimination accuracies of microbial species whose singlecell Raman spectra are obtained with different Raman microscopes. Therefore, we suggest that all single-cell Raman data be acquired with the same instrument and settings.
Furthermore, we did not take into consideration the influence of different physiological states (e.g., exponential vs. stationary state) on discrimination results. Our preliminary studies showed that the intensities of some Raman bands (those of nucleic acids, in particular) vary significantly with cell's physiological state and hence affect the discrimination of microbes by RF.

TROUBLESHOOTING
Problem 1 Fail to trap a microbial cell under the microscope.

Potential solution
Increase laser power for more effective trapping. Small particles suspended in the buffer (e.g., some organic/inorganic precipitate) may sometimes be misidentified as microbial cells. Such non-cellular objects tend to be difficult to trap.

Problem 2
Acquired single-cell Raman spectra suffer from very low S/N.

Potential solution
Capture a microbial cell that is not too far from the cover glass. For long rod-shaped microorganisms, trap a cell by focusing the laser at either end of the cell and not around the center. Avoid mechanical shocks to the microscope to prevent trapped cells from escaping. Trap only one cell because trapping multiple cells results in fluctuating cell positions and weak signal intensity. Note that if the size of a microbe is very small as in the case of Nitrososphaera viennensis, a coccus with a diameter of 0.6-0.9 mm, its Raman spectrum is very likely to have poor S/N because the number of Raman scatterers within the focal volume is small compared to normal-sized cells.

Problem 3
Fail to perform vector normalization.

Potential solution
It is often difficult to normalize data with very different noise levels. To cope with this problem, perform spectral smoothing prior to linear baseline fitting ( Figure 6A). Subtract the baseline from both unsmoothed and smoothed spectra. Then use the norm of the baseline-subtracted, smoothed spectrum to vector-normalize the baseline-subtracted, unsmoothed spectrum.

Problem 4
Features having high importance scores are localized in regions where there are no Raman bands.

Potential solution
This could occur due to inappropriate baseline subtraction or normalization. For the former, adjust the region to be fitted (step 11c and Figure 6A). For the latter, refer to troubleshooting 3.
Overfitting may also cause this issue. RF models may classify microbial cells not by the differences in molecular vibrations characteristic of the cells ('molecular fingerprints'), but by the features that happen to exist only in certain sample group. In such cases, validate the model using a dataset different from that used to create the model (e.g., data acquired on different measurement days) as a test dataset. If groups are very similar (e.g., between strains), you should be aware of the possibility of overfitting.

OPEN ACCESS
Problem 5 Classification accuracy is not high.

Potential solution
Increase the number of decision trees, namely, 'n_estimators'. A caveat is that the use of too many trees can often lead to overfitting. If the number of Raman spectra differs substantially from class to class, specify the 'class_weight' parameter in the RamdomForestClassifier.

RESOURCE AVAILABILITY
Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Shinsuke Shigeto (shigeto@kwansei.ac.jp).

Materials availability
This study did not generate new unique reagents.

Data and code availability
The datasets generated during this study are available at Mendeley Data: https://doi.org/10.17632/ 8cd34fckgt.1.
This paper includes part of the Python code generated during this study. The source code is available at Zenodo: https://doi.org/10.5281/zenodo.7122583.
Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.