Combining Acoustic Bioprinting with AI-Assisted Raman Spectroscopy for High-Throughput Identification of Bacteria in Blood

Identifying pathogens in complex samples such as blood, urine, and wastewater is critical to detect infection and inform optimal treatment. Surface-enhanced Raman spectroscopy (SERS) and machine learning (ML) can distinguish among multiple pathogen species, but processing complex fluid samples to sensitively and specifically detect pathogens remains an outstanding challenge. Here, we develop an acoustic bioprinter to digitize samples into millions of droplets, each containing just a few cells, which are identified with SERS and ML. We demonstrate rapid printing of 2 pL droplets from solutions containing S. epidermidis, E. coli, and blood; when they are mixed with gold nanorods (GNRs), SERS enhancements of up to 1500× are achieved.We then train a ML model and achieve ≥99% classification accuracy from cellularly pure samples and ≥87% accuracy from cellularly mixed samples. We also obtain ≥90% accuracy from droplets with pathogen:blood cell ratios <1. Our combined bioprinting and SERS platform could accelerate rapid, sensitive pathogen detection in clinical, environmental, and industrial settings.


Gold nanorod synthesis and characterization
Hexadecyl(trimethyl)ammonium bromide (CTAB) and sodium oleate (NAOL) coated gold nanorods were synthesized following previously described protocols [1]. The nanorods were cleaned by centrifuging 1.5 mL aliquots twice at (9000 rpm, 20 min), allowing for one wash after synthesis as this has been shown to be adequate to maintain cell viability while preventing nanorod aggregation [2]. Samples were concentrated down to 10 µL to be mixed with cell samples and diluted to a final volume of 200 µL. Absorption spectra were recorded using a Cary 5000 UV-vis-NIR spectrometer. Scanning electron microscopy images were taken using FEI Magellan 400 XHR Scanning Electron Microscope (SEM). Transmission electron microscopy images were taken using FEI Tecnai G2 F20 X-TWIN Transmission Electron Microscope (TEM).

Scanning electron microscopy (SEM) of printed samples
For scanning electron microscope (SEM) imaging, printed samples were imaged after completion of all Raman Spectroscopy. Samples were prepared by evaporating a ∼10 nm layer of 60:40 gold to palladium to allow for better visualization of cells under electron beam illumination. SEM images were taken using FEI Magellan 400 XHR Scanning Electron Microscope.
Bacteria culturing and preparation E. coli, ATCC 25922, and S. epidermidis, ATCC 12228, were grown from frozen stocks on Trypticase Soy Agar 5% Sheep Blood 221239 BD plates. A single colony was seeded in 10 mL Lysogeny broth (LB) culture medium and incubated at 37°C shaking at 300 rpm for 15 hrs using Thermo Scientific MaxQ 4450 incubator. 1.5 mL of culture was washed with water three times at 6000 rpm for 3 min using a mySPINTM 6 Mini Centrifuge. Samples were then concentrated down to 100 µL volumes. The cell count was collected using a Bright-Line Hemacytometer using a 1:5000 dilution of the cell culture stock solution. Stock solutions contained on average ∼1e10 cells/mL.

Preparation of red blood cell solutions
For data collected with purified mouse red blood cells: CD-1 (1CR) purified Mouse Red Blood Cells (RBCs) K2EDTA Gender Unspecified Pooled samples MSE00RBK2-0104095, were purchased from BioIVT in 5mL volumes. RBCs were diluted in a 1:9 v/v mixture of Invitrogen UltraPure 0.5 M EDTA, Invitrogen 15575020, to a final dilution of 1:5000 and cell counts were collected using a Nexcelom Cellometer X2 cell counter.
For data collected with mouse whole blood: CD-1 (1CR) Mouse Whole Blood K2EDTA Gender Unspecified Pooled samples MSE00WBK2-0000627, were purchased from BioIVT in 5mL volumes. The whole blood were diluted in a 1:9 v/v mixture of Invitrogen UltraPure 0.5 M EDTA, Invitrogen 15575020, to a final dilution of 1:5000 and cell counts were collected using a Bright-Line Hemacytometer. Whole blood contained ∼1e10 cells/mL RBCs.

Preparation of mixtures for printing
Printing was completed using 200 µL of solution. All samples were diluted to a final volume of 200 µL in a 1:9 v/v mixture of Invitrogen UltraPure 0.5 M EDTA, Invitrogen 15575020, and Millipore water, unless otherwise noted. For samples with cells and no nanorods, a single concentrated cell solution or a mixture of cell solutions was diluted in aqueous EDTA to a final concentration of 1e9 cells/mL of each cell type in a given mixture. This concentration was chosen to ensure a majority of printed droplets contained at least 1 cell. For samples of cells mixed with nanorods, cleaned, concentrated nanorod solution was first mixed with concentrated cell stock solution, for our final concentration of 1e9 cells/mL per cell type, and then subsequently diluted with our aqueous EDTA solution. All solutions are mixed by inverting our microcentrifuge tubes a minimum of 10 times.
EDTA was chosen as our sample buffer to avoid crystallization upon drying seen with PBS( Supplementary Fig. 6). On top of that, EDTA provides two further advantages for our printed samples. When EDTA-containing droplets dry on a hydrophobic substrate, a central region of aggregated EDTA, cells, and GNRs dries in a much smaller area than that of a full droplet, forcing the cells and GNRs into a much smaller volume, ensuring better coverage of the cells with GNRs ( Supplementary Fig. 9). Furthermore, we hypothesize that the EDTA induces aggregation among the GNRs due to the electrostatic interaction between any residual CTAB on our GNRs and the EDTA [3], as demonstrated in Supplementary Fig. 7. We hypothesize that this clustering, when coupled with the addition of cells, allowed for greater quantities of nanorods to coat the cells, and led to the creation of SERS "hot spots" amongst the aggregated GNRs, providing strong enhancements [4,5].
Finally, we show that the addition of the EDTA and nanorods adds minimal Raman background noise (Supplementary Fig. 8,9,14), making it appropriate for our work in Raman cellular identification. Finally, to further minimize coffee-ring effects from nanorods upon droplet drying, we used vapor deposition to coat our gold-coated slides with a hydrophobic silane layer (3-Aminopropyl)triethoxysilane (APTMS) which allows for a more close packing of our GNRs, providing greater and more uniform enhancement on our cells [6][7][8].

Fabrication of silanized, gold-coated glass slides
The gold substrates used in this work were prepared by evaporating a 5 nm adhesion layer of titanium, followed by 200 nm of gold at a rate of 1 Angstrom/second using a KJ LEsker e-beam evaporator onto piranha cleaned borosilicate glass slides. The gold-coated glass slides were then cleaned with an oxygen plasma, using a Diener Pico Oxygen Plasma Cleaner, for 3 min at 100 W power and ∼2 mbar of pressure, and silanized with 3-(aminopropyl)trimethoxysilane, APTMS, using vapor deposition in order to make the surface more hydrophobic and allow for greater aggregation of the gold nanorods on the cells [8][9][10]. Slides were placed in a 1 liter flask in the presence of 100 µL of APTMS, Sigma-Aldrich 281778-5ML. The flask was then placed in a water bath at 40°C and allowed to react for 1 hr, after which the slides were removed from the flask and placed on a hot plate heated to 40°C for 10 min to allow for the evaporation of loosely bound molecules. We demonstrate that this APTMS layer also provides minimal Raman background noise, making it a great candidate for quick and easy substrate modification for biological Raman analysis ( Supplementary Fig. 8).

Acoustic printing
Acoustic printing was completed using our custom-built ultrasonic, immersion transducer with a center frequency of 147 MHz and a focal distance of 3.5 mm (unless otherwise noted) as determined using a network analyzer, Hewlett Packard 8751A, and through pulse echo measurements taken on an oscilloscope, Keysight InfiniiVision DSOX3054A. The transducer was bound to a quartz, spherical focusing lens. The transducer was mounted on x,y,z manual translation stages, facing downwards, held 3.5 mm above a 303 stainless steel ejection plate with a 1 mm hole. For printing experiments, fluid was pipetted into the gap between the tip of the focusing lens and the ejection plate, held in place through surface tension. During printing experiments, droplets were ejected downwards through this 1 mm hole onto our chosen substrates ( Supplementary Fig. 2).
To generate our droplets, our transducer was powered by a waveform generator, Keysight 33600A Series Trueform Waveform Generator. The waveform generator was connected to a synthesized RF signal generator, Fluke 6062A, which in turn is connected to a power amplifier, Minicircuits ZHL-03-5WF+. Our waveform generator produces a square-wave burst with a repetition frequency of 1 kHz, when operating continuously, at our desired pulse width of 5.5 µs and at a voltage of 1.5 volts, enough to trigger our RF synthesizer. The RF synthesizer generates a sinusoidal wave at 147 MHz and at our desired voltage, which then gets amplified before reaching the transducer. Droplets printed from deionized water were ejected with 0.096 µJ of energy, droplets printed from samples of S. epi and E. coli with and without GNRs were printed with 0.139 µJ of energy, and droplets printed with RBCs with and without GNRs and from mixtures of RBCs, S. epi, E. coli, and GNRs were all printed with 0.386 µJ of energy.
To ensure stable ejection, we monitored ejection using a camera, Allied Vision Guppy Pro F-125 1/3" CCD Monochrome Camera, coupled with a 20x objective pointed at the bottom of our ejection plate. This camera was mounted opposite a strobing LED, also triggered by our waveform generator. We also monitored the acoustic echo using an inline oscilloscope, Keysight InfiniiVision DSOX3054A. To set up our printer, we pipette in 200 µL of fluid, turn on power to our transducer, and ensure that the transducer is in focus by manually adjusting the focal distance of our transducer until we maximize the echo as observed on the oscilloscope. We then vary the output voltage of the RF synthesizer until we stably eject a single droplet without any additional satellite droplets, as observed through our camera feed. We were then ready to pattern print arrays of droplets ( Supplementary Fig. 3, 4).

Pattern printing
Pattern printing was completed using a custom 3D printed substrate holder mounted to two perpendicularly stacked Thorlabs DDS100M 100mm brushless DC linear translation stages controlled by two Thorlabs K-Cube brushless DC servo drivers. Our substrate is mounted ∼1 mm below our ejection plate to minimize droplet translation before it reaches the substrate. A MATLAB, Mathworks, 2018b, script was used to pattern print droplets onto our substrate by controlling both our motorized stages and our waveform generator to trigger droplet ejection at specific substrate locations.

Raman spectroscopy
Raman spectra was collected using the Horiba XploRa confocal Raman microscope. The excitation wavelength for all measurements was 785 nm. The Raman shift from 400 cm -1 to 2000 cm -1 was collected using 600 gr/mm grating. For baseline Raman spectra shown in Fig. 1, laser light was directed to and Raman scattered light was collected from the sample using a 100x LWD, 0.6 NA objective with spot size of 0.83 µm, with laser power at the sample of ∼6.71 mW, and acquisition time of 180 s. For spectra collected from each entire droplet, laser light was directed to and Raman scattered light was collected from the sample using a 10x, 0.25 NA objective with spot size of 2 µm, with laser power at the sample of ∼10.6 mW, and acquisition time of 15 s ( Supplementary Fig. 11). Bacterial NR mixtures were measured within ∼2 hr of sample preparation.

Spectral data processing
Python (Jupyter Notebook) was used to process spectral data. For spectra preprocessing, samples were first thresholded to a minimum intensity of 150 a.u. to remove any spectra with weak signal that likely were collected on the substrate without the presence of cells. We then transformed our data by taking log 10 (y) [11,12] and smoothed the spectra using wavelet denoising [13,14]. To perform our smoothing, we used the denoise wavelet function from the scikitimage Python library: denoise wavelet(y, method='BayesShrink', mode='soft', wavelet levels=1, wavelet='coif3', rescale sigma='True'). We then performed a baseline removal by using a polynomial fit with degree 10. The specific package used and code line is: peakutils.baseline(y, deg=10, max it =1000, tol=0.0001). Note, the need for a higher degree polynomial arises from a typical instrumental background that is difficult to fit with lower degree fits.
Following this baseline correction, Spectra were then individually normalized across all wavenumbers by subtracting the spectral mean and dividing by the standard deviation using the NumPy Python library [15], where y is the array of intensity values across all wavenumbers for each spectra: (ynumpy.mean(y))/numpy.std(y) (Supplementary Fig. 31).
For classification of samples, we further pre-process data by reducing dimensionality of our spectra from 508 to the number of components necessary to account for 90% of our sample variance using the PCA algorithm from Scikit-learn [16] (Supplementary Fig. 19, 21). Classification was performed using a Random Forest Classifier. We first tuned our classifier hyperparameters using a cross-validated grid search to generate optimized parameters. To do this we use Scikit-learn StratifiedShuffleSplit [16] function to randomly split our sample 20 times into an 80:20 train:test split and created a parameter grid for our number of estimators: {50, 100, 150, 200, 250, 300}, max features: {auto, sqrt, log2}, and max depth: {2, 5, 10, 15, 20}. We created our Random Forest Classifier using Scikit-learn, with a min samples split=2, and then performed our grid search using Scikit-learn GridSearchCV, with refit = True, n jobs = 3, and verbose = 190. We then perform a stratified K-fold cross validation (Scikit-learn StratifiedKFold [16] with shuffle=True) of our classifier's performance across 10 splits using these optimal parameters. Finally, we use the Scikit-learn confusion matrix [16] function to plot our results. Intermediate t-SNE projections were plotted using Scikit-learn manifold.TSNE with a perplexity = 10 [16] (Supplementary Fig. 20, 22).
Raman wavenumber importance was performed using Voigt profile perturbation across all spectral wavenumbers. To achieve this, all spectra were first preprocessed as described above. Our spectra of interest (600 spectra across all 6 cellular classes) were partitioned into an 80:20 train/test split using Scikitlearn StratifiedShuffleSplit [16] with 10 splits. We reduce the dimensionality of our training set to 10 components using the PCA algorithm from the Scikitlearn [16]. We train our Random Forest Classifier on our training spectra using optimized hyperparameters determined using a cross-validated grid search as previously described. We iteratively perturb our test set at each wavenumber to determine the relative importance of each wavenumber to accurate spectra classification. To do this, we iterate over each wavenumber in each normalized spectrum in our test set (120 spectra per split). For each wavenumber, we perturb our test spectra with a Voigt profile curve [17,18], varying the intensity of the Voigt function 5 times at each wavenumber for each spectrum to get a large sample set. To generate our Voigt profiles, we first take all spectra in our entire sample set (600 spectra) and shift the intensity at a given wavelength (w) to guarantee that the intensities are positive. We then randomly shuffle all intensities and randomly select one to be used to generate our Voigt profile (voigt intensity). We generate this profile with our half-width at half-maximum (HWHM) of the Lorentzian profile, γ = 2, and the standard deviation of the Gaussian profile, σ = α / np.sqrt(2*np.log(2)), where α = 5. From here, we create our voigt profile = np.real(wofz((x -w + 1j*γ)/σ/np.sqrt(2))) / σ/np.sqrt(2*np.pi), where x is the entire range of wavenumbers, and scale this distribution to range from [0,1]. We utilize the wofz function from the SciPy Python library to implement the Faddeeva function as the Voigt profile is related to the real part of the Faddeeva function. We also utilize various mathematical functions from the NumPy Python library to generate our profile. Voigt profile width was chosen to match peak widths seen in our spectra. To perturb our spectra with this profile, we take each spectra and shift intensities by the minimum, so that all intensities are positive (pos spectrum). We then perturb each wavenumber by our Voigt profile to generate a modified spectrum = pos spectrum*(1-voigt profile) + voigt intensity*voigt profile. We transform this perturbed spectrum with our established PCA and classify it using our trained Random Forest Classifier. We then plot a confusion matrix for each wavenumber using Scikit-learn confusion matrix [16] and generate an accuracy and f1 score using Scikit-learn classification report [16], across all 6000 trials per wavenumber. Finally, we use our confusion matrix to generate the per class performance. See Supplementary Fig. 23 for more.
Supplementary Note 2: Acoustic printing for handling biological samples Acoustic printing works by using ultrasonic waves to eject a droplet from a free surface of fluid. A radio frequency (RF) burst signal is used to excite a transducer at its resonant frequency, generating ultrasonic waves that exert force on the fluid surface [51,52]. When the focus of the transducer is aligned with the liquid-air interface and the intensity of the acoustic field is high enough, the generated radiation pressure will overcome the surface tension and the sound wave gives rise to a mound of fluid from the surface [52,53]. If the energy of the incident wave exceeds the threshold energy, a droplet breaks free from the fluid surface at a velocity of a few meters per second due to the Rayleigh-Taylor instability [51,54]. The droplet diameter has been shown to closely match the diffraction-limited focal width at the liquid-air interface, and as such, the droplet diameter is inversely proportional to the frequency of the transducer, with 5 MHz and 300 MHz ultrasonic waves generating droplet diameters of 300 µm and 5 µm, respectively [52]. (see Fig. 1a). ADE droplet ejection has been well characterized and has great tunability for handling a variety of biological samples [55]. Furthermore, the focused ultrasonic waves completely control the size, speed, and directionality of the ejected droplet and allow for ADE from an open liquid surface. Given that the dimensions of this open liquid surface are much larger than the diameter of the focal spot size, ADE is considered a nozzle-less technology [56]. This holds true for our downwards setup utilizing an ejection plate, given that our focal spot size is ∼2 orders of magnitude smaller than the 1 mm diameter hole [51,52]. As a nozzle-less technology, ADE has unparalleled advantage in biological sample handling as compared with other commercial piezo or thermal inkjet printers that rely on physical flow focusing. In particular, ADE eliminates system clogging and compromised cell viability or biomarker structure due to shear forces generated by nozzles. Additionally, ADE relies on ultrasonic waves to generate droplets, as such, the transducer never has to contact the ejection medium, but rather can propagate through a matched coupling medium, eg. through the bottom of an acoustically "transparent" multiwell plate, with minimal loss of acoustic energy, mitigating risks of sample contamination and loss of sterility [56]. ADE has also gained traction for versatility in setup and ability for high-throughput droplet generation. For a single acoustic ejector, the limiting factor for droplet ejection rate is the dissipation of capillary waves propagating radially outwards on the fluid surface after ejection [51,57]. Advances in ADE have led to improvements in fabrication methods of the focusing lenses and ejector arrays including: spherical lenses in silicon, spherical PZT shells, and fresnel acoustic lenses [51,58]. These advancements have lead to the development of high-throughput ejector arrays greater than 1000 print heads and ejection rates of 25 kHz, allowing for ejection of a 10 mL of fluid in under an hour [58]. Furthermore, these advances have expanded the utility of ADE for biological samples handling to include cellular acoustic printing [55,59,60], biological crystallography [61][62][63], high-throughput screening (HTS) of biological agents [64], and for sample preparation in MALDI [65], highlighting the vast potential for biological analysis with ADE. Supplementary Fig. 1 Photographs of droplets printed with a range of acoustic frequencies. Droplets were printed with 4.8 MHz, 17 MHz, 44.75 MHz, and 147 MHz and had droplet diameters of 300 µm, 84 µm, 44 µm, and 15 µm respectively, highlighting the tunability of acoustic droplet ejection. Scale bars are 500 µm, 200 µm, 100 µm, and 25 µm respectively.
Supplementary Fig. 2 a, Stroboscopic photograph showing a droplet being ejected downwards through the 1 mm hole on our ejection plate. Droplet was ejected from a pool of deionized water using a transducer operating at 147 MHz, Photo was taken with a 12 µs phase delay after the burst was triggered. Scale bar is 100 µm. b, Stroboscopic images of the time evolution of downward droplet ejection through the 1 mm hole at an acoustic frequency of 147 MHz. Droplet shown here is 15 µm in diameter and ∼2 pL in volume. Droplet was ejected with 0.096 µJ of energy with a pulse width of 5.5 µs, and was ejected downwards at ∼3.5 m/s. Scale bar is 100 µm. All images were captured with an image exposure time of 40 ms and a droplet ejection rate of 1 kHz. As such, each frame is composed of 40 droplet ejections.
Supplementary Fig. 3 Photo a, and rendering b, showing acoustic printing setup, including camera with 20x objective, baseplate with 1 mm diameter hole, strobing LED, gold-coated glass slide mounted onto a motorized XY stage, acoustic transducer, and printing fluid (teal).
Supplementary Fig. 4 Schematic showing the acoustic droplet ejection setup. The printing fluid (teal) rests between the focused acoustic transducer and the ejection plate with the 1 mm hole, held in place through surface tension. The droplets are ejected downwards onto a gold-coated glass slide mounted onto a motorized XY stage (stacked single axis stages). The burst signals to the transducer are generated from a function generator, routed through an RF synthesizer, and finally through a power amplifier before reaching the transducer. Ejection and movement of the mounted slide are controlled synchronously using MATLAB code.
Supplementary Fig. 5 Absorption spectra of GNRs used for Raman spectroscopy. Inset shows TEM of GNRs. Scale bar is 50 nm.
Supplementary Fig. 6 Bacterial interrogation across multiple nanorod syntheses and resonance frequencies. We synthesized 4 different batches of GNRs and evaluated bacterial droplets with each batch. a, UV-Vis measurements of the 4 rods showing a range of resonance frequencies between 770 nm and 960 nm. Our chosen nanorods are those listed as NR4. b, Average spectral intensities and standard deviations collected from droplets printed with each of our four GNR batches mixed with S. epi bacteria diluted to a final concentration of 1e9 cells/mL in a 1:9 mixture v/v of Invitrogen UltraPure 0.5M EDTA, Invitrogen 15575020, and Mili-Q purified water onto a silanized, gold-coated glass slide. Spectra were collected with a 10x objective lens with a 0.25 NA and ∼10.6 mW power. Each droplet was interrogated with a time study of 5 time points, with each exposure lasting 15 seconds. 10 droplets were analyzed from each GNR batch for a total of 50 data points across each group. The data shows the primary S. epi peaks around 731 and 1317 cm -1 , with varying max intensities, highlighting both the robustness of our system for spectral collection as well as the potential for improvements through tuning of GNR aspect ratio. Milli-Q purified water. Photograph was taken 5 min after the solution was mixed together. We note that these images highlight that, as expected, the EDTA appears to aggregate the GNRs into clusters. As the S. epi bacteria seem to cause clustering regardless of the presence of the EDTA due to their surface charge, the difference between the sample with and without the EDTA is less noticeable than in the samples of only GNRs.
Supplementary Fig. 9 Raman of background signals. a, Spectra were collected of a goldcoated glass slide, a gold-coated glass slide with an APTMS silane layer, and of a droplet printed onto a gold-coated glass slide with APTMS. Droplets were printed from Invitrogen UltraPure 0.5M EDTA, Invitrogen 15575020, mixed in a 1:9 ratio v/v with Mili-Q purified water. The spectra show that most of the background signal comes from the gold substrate with little additional background from our APTMS deposition and the additional EDTA used in our cell solutions. b, Identical spectra to that shown in a overlaid with a spectrum taken of S. epi bacteria and GNRs suspended in EDTA solution at a concentration of 1e9 cells/mL. The plot highlights that the spectral signal intensity from our bacteria is much higher than that of the background. All spectra were collected with a 10x objective lens with a 0.25 NA and ∼10.6 mW power for 15 s. Supplementary Fig. 10 Droplets were printed from cellular dilution mixture without any cells. Droplets were printed from Invitrogen UltraPure 0.5 M EDTA, Invitrogen 15575020, mixed in a 1:9 ratio v/v with Mili-Q purified water onto a silanized, gold-coated glass slide. Spectra were collected with a 10x objective lens with a 0.25 NA and ∼10.6 mW power for 15 s. The SEM clearly shows minimal spread of the EDTA solution onto the gold-coated slide. The spectra show minimal, and consistent background signals from the EDTA solution. Scale bar is 5 µm.
Supplementary Fig. 11 SEMs of grids printed with cell and GNR mixtures. SEMs show 16 droplets imaged out of a grid of over 400 droplets. Mixtures were printed from cells with and without GNRs diluted in a 1:9 mixture of EDTA solution and Milli-Q water to a final concentration of 1e9 cells/mL. All grids were acoustically printed using a transducer operating at 147 MHz with a 5. Supplementary Fig. 15 Background signal from gold nanorods (GNR). Droplets were printed from a sample containing GNRs suspended without any cells in a 1:9 v/v mixture of Invitrogen UltraPure 0.5 M EDTA, Invitrogen 15575020, and Mili-Q purified water onto a silanized, gold-coated glass slide. Spectra were collected of each droplet with a 10x objective lens with a 0.25 NA and ∼10.6 mW power for 15 s. The SEM shows a droplet containing a few clusters of GNRs clearly distinguishable, highlighting both the presence of the GNRs and the absence of a coffee-ring of nanorods. The plot shows the mean and standard deviation of 100 droplets printed from a sample of S. epi bacteria with GNRs in EDTA solution (in blue, identical to that from Fig. 3), and the mean and standard deviation of 20 droplets printed from the GNR solution without cells. The spectra show that while the GNRs have a background signal, hypothesized to be from any remaining CTAB present in the solution after rinsing the rods, the cells have a clearly distinguishable signal separate from that of the GNRs, similar to our baseline S. epi spectra. Scale bar is 5 µm.
Supplementary Fig. 16 Spectra were collected from a sample of droplets printed from GNRs mixed with mouse RBCs at a final concentration of 1e9 cells/mL diluted in a 1:9 ratio v/v of Invitrogen UltraPure 0.5 M EDTA, Invitrogen 15575020, and Mili-Q purified water onto a silanized, gold-coated glass slide. Spectra were collected with a 10x objective lens with a 0.25 NA and ∼10.6 mW power for 15 s. The first spectra is taken while the focal spot is centered on the droplet while the other is taken when on the silanized gold substrate to the side of the droplet highlighting that our signal is coming directly from the droplet and not from any background material on the slide. Fig. 17 Plot showing the mean and standard deviation of SERS spectra taken from droplets printed from three cell lines (S. epi, E. coli, and RBCs) with and without GNRs. Spectra were collected from 100 and 15 droplets, respectively. The plots highlight both the enhancements generated with the presence of GNRs as well as the variations in peak spectra intensity due to the variations in surface charge density on each cell line and subsequently the cells' varying attraction to the positive surface charge of the GNRs, resulting from the CTAB surfactant on their surface [2].

Supplementary
Supplementary Fig. 18 Data from droplets print from mouse whole blood mixed in with GNRs. Mouse whole blood was purchased with the addition K2EDTA as an anti-coagulating agent for processing and shipment. Blood was then diluted 1:9 with our EDTA solution (1:9 ratio v/v of Invitrogen UltraPure 0.5 M EDTA, Invitrogen 15575020 in DI water) for a final RBC count of 1.1 e9 cells/mL, closely approximating the RBC count of samples printed with RBC mixtures. GNRs at a concentration matching that all other cellular samples was added to the solution and then samples were printed onto a gold-coated slide and Raman spectra were collected. a, SEMs showing stable printing of solutions from diluted mouse whole blood mixed with GNRs. Top image shows a 6x6 grid of droplets as a subset of the hundreds printed. The remaining two SEMs show two representative droplets from the sample with 1 and 7 RBC, respectively. Scale bars are 50 µm and 5 µm. b, Mean spectra standard deviation of 100 spectra collected from samples printed from mouse RBCs and mouse whole blood, both diluted in aqueous EDTA to a final RBC count of 1e9 cells/mL mixed with GNRs. Spectra are shown after baseline correction c, Normalized mean spectra from data shown in b. Normalized means show that while the whole blood exhibits similar peaks to that of the RBCs, there are many more minor peaks in the spectrum, highlighting expected sample complexity from the whole blood sample. Supplementary Fig. 19 Plot of the percentage of variance attributed to each principal component and the cumulative explained variance over 50 components. The green line indicates the number of PCA components necessary to capture 90% of all explained variance in our samples. For all 300 spectra from our single cell-line droplets, we demonstrate that we can account for at least 90% of all variance with 24 components generated from all 508 wavenumber features in our spectra. Supplementary Fig. 20 Plots showing a 2-component, t-distributed stochastic neighbor embedding projection (t-SNE) with perplexity = 10 across our 3 single cell-line classes. Data is plotted a, with data inclusive of all wavenumber features and b, after performing a 24component PCA for dimensionality reduction. Plots show relative clustering of our classes and minimal variation to clustering after dimensionality reduction. Supplementary Fig. 21 Plot of the percentage of variance attributed to each principal component and the cumulative explained variance over 50 components. The green line indicates the number of PCA components necessary to capture 90% of all explained variance in our samples. For all 600 spectra from our 3 single cell-line droplet classes and our 3 cell mixture classes, we demonstrate that we can account for at least 90% of all variance with 30 components generated from all 508 wavenumber features in our spectra. Supplementary Fig. 22 Plots showing a 2-component, t-distributed stochastic neighbor embedding projection (t-SNE) with perplexity = 10 across all 6 of our classes. Data is plotted a, with data inclusive of all wavenumber features and b, after performing a 30-component PCA for dimensionality reduction. Plots show relative clustering of our classes and minimal variation to clustering after dimensionality reduction.
Supplementary Fig. 23 Feature Validation. We perform feature validation on our spectra to determine which wavenumbers and spectral bands are most important for our classifier. We take our 600 spectra across all 6 cellular classes and split the samples using a stratified shuffle split into an 80:20 train/test split. For each spectrum in our training set, we iteratively perturb the spectrum at each wavenumber. After each perturbation, we calculate the classification accuracy and compare with our baseline accuracy. Every wavenumber of each spectrum in the test set is perturbed 5 times and all results are averaged for our final feature extraction. Spectra were perturbed with a normalized Voigt profile. Line width chosen to roughly match peak widths seen in our spectra. a, plot showing an example Voight curve (blue), unperturbed example spectrum from our dataset (red), and perturbed spectrum (green). b, Voigt profile intensities were chosen through random sampling of all spectra in our training set. Plot shows an example of a spectrum from our sample set with 100 different perturbations. c, Heatmap highlighting feature validation performed to determine relative weight of spectral wavenumbers in our Random Forest classification. Heatmap is overlaid with a plot of mean and standard deviation of the perturbed classification accuracy (red) and f1 score (blue) calculated across all trials. Mean accuracy is plotted in green. Wavenumbers with lower accuracies are shown to be critical features, as random perturbations in these regions are highly correlated with decreases in classification accuracy. d, Plot of the mean classification accuracy broken down into accuracies across each of our cellular and mixture classes. Supplementary Fig. 24 Classification using spectral feature bands of interest evaluated across 600 spectra collected from single cell-line droplets of S. epi, E. coli, and mouse RBCs mixed with GNRs, and our 3 cell mixtures. a, Heatmap presented in Supplementary Fig.  22, highlighting feature importance calculations performed to determine relative weight of spectral wavenumbers in our random forest classification. Heatmap is overlaid with 3 bands representing key spectral bands used by our classifier. We further demonstrate that these bands are primarily responsible for our classification accuracies by preprocessing our spectra by removing spectral features outside these bands (420-522 cm -1 , 700-775 cm -1 , 1200-1454 cm -1 ). We then reduced the dimensionality of our remaining features using an 8-component PCA as previously reported. b, Plot of the percentage of variance attributed to each principal component and the cumulative explained variance over 50 components. The green line indicates the number of PCA components necessary to capture 90% of all explained variance in our samples. For this sample set taking only specific wavenumber bands from our spectra, we demonstrate that we can account for at least 90% of all variance with only 8 components generated from all 508 wavenumber features in our spectra. c, Finally, we use our previously described random forest classifier on our samples and perform a stratified Kfold cross validation of our classifier's performance across 10 splits. Results are plotted on a normalized confusion matrix. We show that we achieve ≥ 81% classification accuracy across all samples as compared with the ≥ 87% classification accuracy achieved when evaluating the entire spectra window from 400-1700 cm -1 . These results further validate our feature recognition model. Furthermore, they pave the way for future development of low cost POC systems by demonstrating that the use of low-cost spectrometers with limited spectral windows may be possible for such diagnostic work. Supplementary Fig. 25 Plot showing the mean and standard deviation of SERS spectra taken from droplets printed from our 6 droplet classes: three single-cell line classes (S. epi, E. coli, and RBCs) and three mixture classes (equal mixtures of S. epi and RBCs, E. coli and RBCs, and S. Epi, E. coli, and RBCs) all diluted to a final concentration of 1e9 cells/mL of each cell type in our aqueous EDTA solution and mixed with GNRs. Spectra were collected from 100 droplets for each class.

Supplementary
Supplementary Fig. 27 We further analyze samples printed from mixtures of S. epi bacteria, RBCs, and GNRs to understand how our ML algorithms would handle droplets printed with a greater number of RBCs as compared to S. epi. For the 100 droplets analyzed in Figure 4 from the main text, we count the number of S. epi bacteria and RBCs and a, plot ratio of S. epi:RBCs in each of the droplets. Of note, 8 droplets had 0 RBCs and therefore are not included on the plot. b, There were 10 droplets with a greater number of RBCs than S. epi in each droplet. For these droplets, we train our ML algorithm on the remaining 590 spectra in our 6 cellular classes. We then test it on these 10 spectra. The pre-processed spectra are plotted above, with the color indicating the classifier prediction. We got a 90% classification accuracy, with 1 spectra falsely classified as containing both S. epi and E. coli bacteria along with RBCs. c, We show two representative SEMs from these 10 droplets, with the RBCs and S. epi false colored in red and blue, respectively. The scale bar is 5 µm.
These results are indicative of what we might see in a more clinically relevant sample with a greater concentration of RBCs than bacteria. We see that even in our limited dataset, we get comparable classification accuracies to samples with greater bacterial numbers, showing promise for our platform's ability to detect lower bacterial concentrations.
Supplementary Fig. 30 Future vision for Raman hyperspectral imaging. a, Calculated plot of the classification accuracy for each independent sample class, highlighting wavenumber importance for each individual class. Wavenumbers with lower accuracies are shown to be critical features as random perturbations are highly correlated with decreases in classification accuracy. b,Calculated plot of the mean classification accuracy across all 6 sample classes, calculated across all perturbation trials. Plot is shown as a rainbow spectrum to correspond with simulated hyperspectral images shown in c. Plot is identical to that shown in the main text Fig. 4d. c, Simulated image showing a grid array of 25 droplet SEM. 23 droplets contain only RBCs and GNRs while 2 droplets contain a mix of S. epi bacteria, RBCs, and GNRS. The image on the left shows the droplets containing only RBCs lighting up when interrogated with 816 nm light, while the image on the right shows droplets containing a mixture of bacteria and RBCs lighting up when interrogated with light at 833 nm. Chosen wavelengths correspond to accuracy dips associated with each droplet class. Scale bar is 10 µm.
Supplementary Fig. 31 Spectral Preprocessing. Plot showing a sample spectra taken from our dataset of spectra collected from droplets printed with E. coli bacteria and GNRs. Plot shows (from top to bottom) the raw spectrum, the spectrum after a log 10 transformation, spectrum after smoothing using a wavelet denoising, spectrum with baseline correction, and finally the normalized spectrum.