Comparison of spontaneous Raman spectrometers –

AI usage disclosure

We used ChatGPT to help write code and help clarify and streamline text that we wrote. We also used Claude to help write code and help clarify and streamline text that we wrote.

Motivation

At Arcadia, we’re interested in using spontaneous Raman spectroscopy as an agnostic, high-throughput tool for biological phenotyping. To implement this approach effectively, we evaluated a variety of Raman spectroscopy systems to determine which would best align with our research requirements. By testing the same research-relevant specimens across all instruments, we could directly compare their performance characteristics, including resolution, sensitivity, and signal-to-noise ratios under real experimental conditions. To facilitate data processing and characterization, we also created an open-source Python package, ramanalysis, to process, normalize, and help interpret the spectral data from different manufacturers. We hope that others getting started with similar systems can use this notebook and the open-source Python package associated with it to assist in their analyses. Our analysis has several caveats based on the acquisition parameters and samples we used, and as such, isn’t intended to generate a definitive ranking or absolute comparison of the instruments themselves.

Instrumentation

Our goal for this study was to test a variety of spontaneous Raman systems with various capabilities to determine which instrument is best suited for our samples and applications. In total, we tested five different systems:

• Horiba MacroRAM
• Renishaw inVia Qontor
• Wasatch WP 785X
• Wasatch WP 532X
• OpenRAMAN

Prior to the collection of a full dataset, we conducted preliminary measurements on a subset of samples to calibrate and optimize the acquisition parameters for each spectrometer. We provide these acquisition parameters in the tables below. For each sample, the laser power listed in the tables is either provided by the instrument or measured at the sample surface. Unless stated otherwise, we focused the laser spot by eye when obtaining measurements and collected spectra from cells growing on agar plates. Note that the OpenRAMAN is a spectrometer we’ve previously used and described in our research (Avasthi et al., 2024; Braverman et al., 2025).

Sample preparation

Acetonitrile

We chose acetonitrile, a common organic solvent, as a reference material for calibration. It has multiple easily resolvable, well-documented peaks in regions of Raman spectra relevant to our sample set. We used acetonitrile, ≥ 99.5% from Avantor/VWR, and analyzed it in a capped quartz cuvette.

Chlamydomonas strains

We compared three wild-type strains of Chlamydomonas grown on two types of solid media for a total of six treatments. There were two strains of Chlamydomonas reinhardtii, CC-124 (mating type -) and CC-125 (mating type +), and one strain of Chlamydomonas smithii, CC-1373 (mating type +).

Media components

We made all solid media with 1.5% agar. We made “water” plates with ultrapure water from a Milli-Q EQ 7008 system equipped with a 0.22 μm filter. Individual media components are as follows:

• TAP (tris-acetate-phosphate) media
  We purchased complete liquid media from UTEX.
• M-N (minimal media lacking nitrogen) media
  Made in-house following the above protocol.

Protocols

We maintained cells on TAP or M-N media plates with 1.5% agar under a 12/12 h light-dark cycle¹ at ambient temperature for 5–15 days before use. Cells were densely streaked to cover about half of the plate, creating a lawn. To acquire spectra from cells, we focused the laser on the portion of the plate with cells. To acquire background spectra from media and agar, we focused the laser on the area without cells.

¹ Note that we didn’t account for the timing within the light phase of the cycle when measurements were taken.

Layout of Chlamydomonas plates.

To prepare slides with cells for microscopy (this was done only for the Renishaw InVa Qontor), we drew a wax circle on a #1.5 glass coverslip to create a hydrophobic barrier for each sample, making six total. We pipetted 3 µL of either liquid TAP or M-N media within the wax circles, resulting in three coverslips with each media type. We then transferred a small number of cells to the liquid droplet using a clean pipette. After transfer, clumps of cells were removed with forceps, resulting in a more continuous distribution of cells in each droplet. We repeated this for each of the three strains. Then, using forceps, we flipped each coverslip, placed it onto an ethanol-cleaned glass slide, and sealed each edge with melted VALAP (1:1:1 mixture of Vaseline, lanolin, and paraffin wax) using a paintbrush.

Analysis & results

Characterizing spectral resolution with acetonitrile

We acquired and analyzed measurements of the same sample of acetonitrile on each instrument to assess system performance. By analyzing the spectra and comparing them to established references, we’re able to characterize each instrument’s approximate spectral resolution and accuracy (the accuracy of the observed peak positions with respect to reference peaks).

We’ll start by loading the acetonitrile spectrum from each instrument as well as the standard Raman peaks for acetonitrile (Shimanouchi, 1972). While five standard peaks are listed in the NIST database as medium, strong, or very strong, the peak at 2,999 cm^–1 can be hard to detect, particularly with 532 nm excitation. For this reason, we’ll only characterize the four most prominent peaks and fit each with a Voigt model² (a convolution of a Gaussian distribution with a Lorentzian distribution) with a window size of 25.

² A Voigt profile is commonly used in peak fitting applications within spectroscopy because it accurately models spectral line shapes by accounting for both Doppler broadening (Gaussian component, due to thermal motion) and Lorentzian broadening (from collisional or lifetime effects). This makes it ideal for fitting experimental spectra where multiple broadening mechanisms contribute to the observed peak shape.

from analysis.load_spectra import load_acetonitrile_spectra
from lmfit.models import VoigtModel
from ramanalysis.calibrate import ACETONITRILE_PEAKS_CM1

# Load acetonitrile spectra and their corresponding spectrometers
acetonitrile_spectra, acetonitrile_dataframe = load_acetonitrile_spectra()

# Set parameters for peak finding
num_peaks = 4
window_size = 25
wavenumber_range_cm1 = (750, 3100)

Figure 1 shows the results of the peak fitting. Because the quartz cuvette of acetonitrile yields relatively strong and clean spectra, very little preprocessing is applied to the spectra. If you expand the code cell below, you’ll see that the only preprocessing steps applied are cropping the spectra to the wavenumber range (750–3,100 cm^–1) and min-max normalization. Only the OpenRAMAN spectrum is smoothed via median filtering to suppress hot pixels and read noise from the detector.

Note

Another strategy to deal with hot pixels from our acquisition would have been to acquire and subtract a dark spectrum (a spectrum with the laser powered off) to capture only the detector’s fixed pattern noise.

import numpy as np
import pandas as pd
import plotly.graph_objects as go
from analysis.plotting import darken, get_custom_colorpalette, get_default_plotly_layout
from plotly.subplots import make_subplots
from ramanalysis.peak_fitting import find_n_most_prominent_peaks

# Store peak fitting results
peak_fitting_records = []

# Create figure
color_palette = get_custom_colorpalette()
num_rows = acetonitrile_dataframe.groupby(["λ_nm"]).count()["instrument"].max().item()
num_cols = acetonitrile_dataframe["λ_nm"].nunique()
fig = make_subplots(
    rows=num_rows,
    cols=num_cols,
    shared_xaxes=True,
    shared_yaxes=True,
    horizontal_spacing=0.02,
    vertical_spacing=0.05,
    column_titles=["785 nm", "532 nm"],
)

# Map each (instrument, wavelength) combo to a subplot
subplot_indices = {  # (instrument, λ_nm) -> (col, row)
    ("horiba", 785): (1, 1),
    ("renishaw", 785): (1, 2),
    ("wasatch", 785): (1, 3),
    ("openraman", 532): (2, 2),
    ("wasatch", 532): (2, 3),
}

for i, dataframe_row in acetonitrile_dataframe.iterrows():
    # Get the raw spectrum and the instrument and wavelength with which it was recorded
    wavelength_nm = dataframe_row["λ_nm"]
    instrument = dataframe_row["instrument"]
    raw_spectrum = acetonitrile_spectra[i]
    subplot_col, subplot_row = subplot_indices[(instrument, wavelength_nm)]

    # Apply median filtering on OpenRAMAN data to filter out hot pixels from the detector
    if instrument == "openraman":
        raw_spectrum = raw_spectrum.smooth(kernel_size=3)
    # Crop and normalize spectra for peak detection
    preprocessed_spectrum = raw_spectrum.between(*wavenumber_range_cm1).normalize()

    # Plot each acetonitrile spectrum
    name = f"{instrument.capitalize()}" if instrument != "openraman" else "OpenRAMAN"
    fig.add_trace(
        go.Scatter(
            x=preprocessed_spectrum.wavenumbers_cm1,
            y=preprocessed_spectrum.intensities,
            name=name,
            hoverinfo="skip",
            marker={"color": color_palette[(instrument, wavelength_nm)]},
        ),
        row=subplot_row,
        col=subplot_col,
    )

    # Find and get the indices of the most prominent peaks in the spectrum
    peak_indices = find_n_most_prominent_peaks(
        preprocessed_spectrum.intensities, num_peaks=num_peaks
    )
    for i_peak in peak_indices:
        # Create window around each peak
        peak_cm1, peak_height = preprocessed_spectrum.interpolate(i_peak)
        window_range = (i_peak - window_size / 2, i_peak + window_size / 2)
        window = np.linspace(*window_range, num=500)
        # Interpolate the window along the spectrum
        x_data_cm1, y_data = preprocessed_spectrum.interpolate(window)

        # Fit the model to each peak
        model = VoigtModel()
        initial_fit_params = model.make_params(amplitude=1, center=peak_cm1, sigma=1, gamma=0.5)
        fit_result = model.fit(y_data, initial_fit_params, x=x_data_cm1)

        # Find the offset from acetonitrile reference peaks
        peak_center = fit_result.params["center"].value
        reference_peak_index = np.abs(ACETONITRILE_PEAKS_CM1 - peak_center).argmin()
        reference_peak = ACETONITRILE_PEAKS_CM1[reference_peak_index]
        offset = peak_center - reference_peak

        # Record the peak fitting results
        peak_fitting_results = {**dataframe_row.to_dict(), **fit_result.params.valuesdict()}
        peak_fitting_results["reference"] = reference_peak
        peak_fitting_results["offset"] = offset
        peak_fitting_records.append(peak_fitting_results)

        # Plot and annotate each peak
        annotation_text = f"{peak_fitting_results['fwhm']:.1f} cm<sup>-1</sup>"
        hovertext_fields = ["fwhm", "height", "offset"]
        hovertext = [
            f"{field}: {peak_fitting_results[field]:.3g}<br>" for field in hovertext_fields
        ]
        fig.add_trace(
            go.Scatter(
                x=[peak_center],
                y=[peak_fitting_results["height"]],
                showlegend=False,
                hovertext="".join(hovertext),
                text=annotation_text,
                mode="markers+text",
                textposition="top center",
                marker={"color": darken(color_palette[(instrument, wavelength_nm)])},
            ),
            row=subplot_row,
            col=subplot_col,
        )

        # Plot fitted Voigt profile for each peak
        x_data = np.linspace(peak_center - window_size, peak_center + window_size, 500)
        model_params = [peak_fitting_results[p] for p in ["amplitude", "center", "sigma", "gamma"]]
        y_data = model.func(x_data, *model_params)
        fig.add_trace(
            go.Scatter(
                x=x_data,
                y=y_data,
                showlegend=False,
                hoverinfo="skip",
                marker={"color": darken(color_palette[(instrument, wavelength_nm)])},
            ),
            row=subplot_row,
            col=subplot_col,
        )

# Create pandas DataFrame of peak fitting results
acetonitrile_peaks_dataframe = pd.DataFrame.from_records(peak_fitting_records).drop("gamma", axis=1)

# Configure plotly layout
layout = get_default_plotly_layout()
fig.update_layout(layout, height=450)
fig.update_xaxes(range=[wavenumber_range_cm1[0] - 100, wavenumber_range_cm1[1] + 100])
fig.update_yaxes(range=[-0.05, 1.3])
fig.update_xaxes(title_text="Wavenumber (cm<sup>-1</sup>)", row=num_rows, col=1)
fig.update_xaxes(title_text="Wavenumber (cm<sup>-1</sup>)", row=num_rows, col=2)
fig.update_yaxes(title_text="Intensity (normalized)", row=2, col=1)
fig.show()

Figure 1: Acetonitrile spectra from each instrument. The four most prominent peaks for each spectrum are detected and fit using a Voigt profile. Hovering the cursor over each peak lets you see the measured FWHM, height, and offset from the reference peak value.

By analyzing the same acetonitrile sample across instruments, we can control for sample-dependent variables, allowing us to make meaningful comparisons of the relative performance of each spectrometer based on full width at half maximum (FWHM)³ measurements. The reference peak at ~918 cm^–1 is comparatively isolated, so the FWHM won’t be impacted by other components of the spectrum. Let’s choose this peak to compare FWHM values, noting that the observed FWHM is influenced by both instrumental factors (e.g., quality of optical components) and sample characteristics (e.g., molecular dynamics and environmental conditions).

³ While the FWHM serves as a useful indicator of spectral quality, it’s important to note that true spectral resolution is defined by the minimum separation at which two adjacent peaks can be distinguished rather than by individual peak widths alone.

acetonitrile_peaks_dataframe.query("reference == 918").set_index(["instrument", "λ_nm"]).round(1)

		amplitude	center	sigma	fwhm	height	reference	offset
instrument	λ_nm
horiba	785	10.9	918.5	4.7	11.7	0.8	918	0.5
renishaw	785	6.2	921.8	3.2	8.2	0.7	918	3.8
wasatch	785	6.5	921.4	4.8	11.9	0.5	918	3.4
openraman	532	4.0	921.4	8.7	21.0	0.2	918	3.4
wasatch	532	1.9	927.1	5.9	14.5	0.1	918	9.1

Table 1: Peak fitting results from acetonitrile spectra collected by each instrument.

The peak fitting results in Table 1 show that the Renishaw has the lowest FWHM of any instrument at 8.2 cm^-1. The Horiba and the 785 nm Wasatch are close behind at 11.7 cm^-1 and 11.9 cm^-1, respectively. Additionally, we found that the Horiba is the best calibrated in this spectral range as it has the smallest offset in peak position from the reference. While spectral accuracy is important, we’re less concerned here, given that we can always recalibrate spectra as part of a post-processing workflow using reference materials. Furthermore, for most of our research, we’re more interested in relative differences among spectra than absolute differences.

Let’s continue the comparison by exploring spectra from biological samples.

Characterizing SNR with C. reinhardtii spectra

We’re shifting our focus from a standard sample (acetonitrile) to biological samples, specifically unicellular algae, to determine how the different spectrometers perform for our experiments of interest. We collected Raman spectra of Chlamydomonas reinhardtii (CC-124) — a model organism we frequently use in the lab — from each instrument. Since biological samples introduce more complexity than pure compounds, we need to assess the signal-to-noise ratio (SNR) of these spectra to determine which instruments offer the highest sensitivity and are best suited for detecting subtle spectral features. This characterization will help us identify the most sensitive spectrometer for future experiments where detecting faint biochemical signals is critical. In addition to SNR measurements, we’ll make (relatively crude) peak-finding calculations to get a rough sense of the number of features we can resolve from each spectrum. We’ll use the find_peaks function from SciPy, which finds the local maxima in a 1D signal, and filter this set of local maxima to the peaks we’re interested in resolving by defining parameters such as the relative prominence and the distance between peaks.

We’ll also need to perform baseline subtraction to derive meaningful SNR measurements, as raw Raman spectra often contain fluorescence and other background signals that obscure meaningful peaks. In a separate analysis, we evaluated different baseline estimation algorithms. We found that adaptive smoothness penalized least squares (asPLS) performed the best at removing background while preserving key spectral features (Zhang et al., 2020). We use the implementation of asPLS in pybaselines with the smoothing parameter changed to lam = 1e4, as we found empirically that this value provided better baseline estimation for our dataset.

from analysis.load_spectra import load_cc124_tap_spectra
from pybaselines.whittaker import aspls
from ramanalysis import RamanSpectrum
from scipy.signal import find_peaks

cc124_tap_spectra, cc124_tap_dataframe = load_cc124_tap_spectra()


# Define a function for calculating peak SNR
def compute_peak_snr(
    spectrum: RamanSpectrum,
    noisy_region_cm1: tuple[float, float],
) -> float:
    """Peak SNR measurement."""
    noise = spectrum.normalize().between(*noisy_region_cm1).intensities.std()
    peak_snr = 1 / noise
    return peak_snr


# Parameters for baseline subtraction
lam = 1e4
wavenumber_range_cm1 = (520, 3200)

# Parameters for SNR calculation
noisy_region_cm1 = (1700, 2500)

# Parameters for peak finding
relative_prominences = {
    ("horiba", 785): 0.5,
    ("renishaw", 785): 0.1,
    ("wasatch", 785): 0.1,
    ("openraman", 532): 0.5,
    ("wasatch", 532): 0.1,
}
peak_separation = 10

Figure 2 shows the peak SNR measurements of each instrument and the approximate number of peaks.

# Store peak SNR measurement results
snr_measurements = []

# Create figure
num_rows = cc124_tap_dataframe.groupby(["λ_nm"]).count()["instrument"].max().item()
num_cols = cc124_tap_dataframe["λ_nm"].nunique()
fig = make_subplots(
    rows=num_rows,
    cols=num_cols,
    shared_xaxes=True,
    shared_yaxes=True,
    horizontal_spacing=0.02,
    vertical_spacing=0.05,
    column_titles=["785 nm", "532 nm"],
)

# Map each (instrument, wavelength) combo to a subplot
subplot_indices = {  # (instrument, λ_nm) -> (col, row)
    ("horiba", 785): (1, 1),
    ("renishaw", 785): (1, 2),
    ("wasatch", 785): (1, 3),
    ("openraman", 532): (2, 2),
    ("wasatch", 532): (2, 3),
}

for i, dataframe_row in cc124_tap_dataframe.iterrows():
    # Get the raw spectrum and the instrument and wavelength with which it was recorded
    wavelength_nm = dataframe_row["λ_nm"]
    instrument = dataframe_row["instrument"]
    raw_spectrum = cc124_tap_spectra[i]
    subplot_col, subplot_row = subplot_indices[(instrument, wavelength_nm)]

    # Crop and normalize spectra for baseline estimation
    normalized_spectrum = raw_spectrum.between(*wavenumber_range_cm1).normalize()
    # Apply median filtering on OpenRAMAN data to filter out hot pixels from the detector
    if instrument == "openraman":
        normalized_spectrum = normalized_spectrum.smooth(kernel_size=3)

    # Estimate and subtract baseline
    baseline_estimate, _params = aspls(normalized_spectrum.intensities, lam=lam)
    baseline_subtracted_intensities = normalized_spectrum.intensities - baseline_estimate
    baseline_subtracted_spectrum = RamanSpectrum(
        wavenumbers_cm1=normalized_spectrum.wavenumbers_cm1,
        intensities=baseline_subtracted_intensities,
    ).normalize()

    # Crop out saturated region from Wasatch 532 nm spectrum
    if (instrument == "wasatch") and (wavelength_nm == 532):
        baseline_subtracted_spectrum = baseline_subtracted_spectrum.between(0, 1800).normalize()

    # Calculate peak SNR
    snr = compute_peak_snr(baseline_subtracted_spectrum, noisy_region_cm1)
    snr_measurements.append(snr)

    # Find prominent peaks
    peaks, _peak_properties = find_peaks(
        baseline_subtracted_spectrum.intensities,
        relative_prominences[(instrument, wavelength_nm)],
        distance=peak_separation,
    )
    peaks_cm1, peak_heights = baseline_subtracted_spectrum.interpolate(peaks)

    # Plot each raw (normalized) spectrum
    fig.add_trace(
        go.Scatter(
            x=normalized_spectrum.wavenumbers_cm1,
            y=normalized_spectrum.intensities,
            hoverinfo="skip",
            showlegend=False,
            marker={"color": color_palette[(instrument, wavelength_nm)]},
            opacity=0.5,
        ),
        row=subplot_row,
        col=subplot_col,
    )

    # Plot each baseline estimate
    fig.add_trace(
        go.Scatter(
            x=normalized_spectrum.wavenumbers_cm1,
            y=baseline_estimate - 0.05,
            hoverinfo="skip",
            showlegend=False,
            marker={"color": darken(color_palette[(instrument, wavelength_nm)])},
            opacity=0.5,
        ),
        row=subplot_row,
        col=subplot_col,
    )

    # Plot baseline subtracted spectra
    name = f"{instrument.capitalize()}" if instrument != "openraman" else "OpenRAMAN"
    fig.add_trace(
        go.Scatter(
            x=baseline_subtracted_spectrum.wavenumbers_cm1,
            y=baseline_subtracted_spectrum.intensities,
            name=name,
            hoverinfo="skip",
            marker={"color": color_palette[(instrument, wavelength_nm)]},
        ),
        row=subplot_row,
        col=subplot_col,
    )

    # Plot peaks
    fig.add_trace(
        go.Scatter(
            x=peaks_cm1,
            y=peak_heights,
            name=name,
            mode="markers",
            marker={"color": color_palette[(instrument, wavelength_nm)]},
            showlegend=False,
        ),
        row=subplot_row,
        col=subplot_col,
    )

    # Plot SNR and num peaks annotation
    text = f"Peak SNR: {snr:.1f}<br>Num peaks: {peaks.size:.0f}"
    fig.add_trace(
        go.Scatter(
            x=[2000],
            y=[0.6],
            showlegend=False,
            text=text,
            mode="text",
            textposition="middle right",
            name=name,
        ),
        row=subplot_row,
        col=subplot_col,
    )

# Create pandas DataFrame of peak fitting results
cc124_tap_dataframe["snr"] = snr_measurements

# Configure plotly layout
fig.update_layout(layout, height=450)
fig.update_xaxes(range=[wavenumber_range_cm1[0] - 100, wavenumber_range_cm1[1] + 100])
fig.update_yaxes(range=[-0.05, 1.3])
fig.update_xaxes(title_text="Wavenumber (cm<sup>-1</sup>)", row=num_rows, col=1)
fig.update_xaxes(title_text="Wavenumber (cm<sup>-1</sup>)", row=num_rows, col=2)
fig.update_yaxes(title_text="Intensity (normalized)", row=2, col=1)
fig.show()

Figure 2: Peak SNR measurements and detected peaks in C. reinhardtii spectra from each instrument. Each subplot shows the baseline-subtracted spectrum on which the peak SNR measurements are made. Behind the baseline-subtracted spectra, the raw C. reinhardtii spectra and estimated baselines are shown with semi-transparency. Note that the estimated baselines are slightly recessed from the raw spectra for visualization. Clicking on the instrument name in the legend toggles the display of the baseline-subtracted spectra to view the raw spectrum and estimated baseline more easily.

These measurements indicate that the Renishaw and 785 nm Wasatch systems provide the highest peak SNR values. The OpenRAMAN exhibits the lowest SNR, which aligns with expectations given its relatively low-cost components. In particular, the OpenRAMAN has by far the weakest laser power of any system and likely the most read noise from the detector; in contrast to the other systems, the camera is not cooled. Perhaps not surprisingly, the two instruments with the highest SNR also have the highest number of discernible spectral features.

The excitation wavelength greatly affects the usable signal. While it also has a relatively high SNR, the number of features in the 532 nm Wasatch spectrum is much lower. This suggests that 532 nm may be less optimal for analyzing these particular types of samples. It should be noted that the number of peaks detected in each spectrum is only an approximation.

Now, we’ll move on to applications requiring large spectral datasets, such as batch processing and machine learning (ML) techniques. These approaches will help us evaluate the reproducibility of each instrument and assess their sensitivity to subtle variations in differentiating strains and growth conditions of our cell samples, for example.

Batch preprocessing of Chlamydomonas spectra

The C. reinhardtii spectra characterized above come from a much larger dataset of Chlamydomonas spectra. Three strains of Chlamydomonas, representing two species, were cultured in two different media conditions, as described in Section 3. With each instrument, we acquired Raman spectra from plates (or, in the case of the Renishaw, glass slides) of cells representing each of these six treatments.

Let’s load the full dataset of spectra and associated metadata — a random sample of metadata from ten spectra from this dataset is shown in Table 2.

from analysis.load_spectra import load_chlamy_spectra

# Load all the spectra of Chlamydomonas and the corresponding metadata
chlamy_spectra, chlamy_dataframe = load_chlamy_spectra()

# Show a sample of spectra metadata
print(f"Total number of Chlamydomonas spectra: {len(chlamy_dataframe)}")
chlamy_dataframe.sample(10, random_state=57).reset_index(drop=True)

Total number of Chlamydomonas spectra: 536

	instrument	λ_nm	species	strain	medium
0	renishaw	785	C. reinhardtii	CC-124	TAP
1	renishaw	785	C. reinhardtii	CC-125	TAP
2	wasatch	532	C. smithii	CC-1373	TAP
3	openraman	532	C. reinhardtii	CC-125	TAP
4	openraman	532	C. smithii	CC-1373	M-N
5	wasatch	532	C. smithii	CC-1373	TAP
6	wasatch	785	C. reinhardtii	CC-125	M-N
7	wasatch	785	C. reinhardtii	CC-124	M-N
8	openraman	532	C. reinhardtii	CC-125	M-N
9	wasatch	785	C. reinhardtii	CC-125	M-N

Table 2: Metadata corresponding to a random sample of ten spectra from the Chlamydomonas dataset comprised of spectra from each instrument.

Since these spectra capture biochemical variations between species and environmental influences, we aim to apply machine-learning techniques to distinguish between groups. However, before diving into the modeling, we need to perform batch preprocessing to correct for background fluorescence, noise from the detectors, and other experimental inconsistencies. This ensures that our spectral differences reflect biological variability (as much as possible) rather than artifacts from acquisition conditions, improving the reliability of our classification models. The preprocessing pipeline is comprised of the following steps:

1. Denoising: Smooth spectra with Savitzky-Golay filter.
2. Cropping: Crop spectra to the region of interest (300–1,800 cm^–1).
3. Baseline subtraction: Fit and remove the fluorescence baseline using the asPLS algorithm.
4. Normalization: Apply min-max normalization to scale the spectra for consistency across samples.

Because data acquired with the Renishaw instrument was collected from glass slides, we need to perform an additional background subtraction. Glass has a distinct Raman spectrum that overlaps with spectral regions of interest from our biological samples (Fikiet et al., 2020), requiring a separate correction to prevent it from confounding our analysis. Figure 3 shows the results of our preprocessing pipeline on the Chlamydomonas dataset.

Note

When we collected data on the Horiba instrument, our focus wasn’t on ML classification tasks, so we didn’t design the dataset with that goal in mind. As a result, we don’t have a sufficient number of samples to support meaningful classification analysis. Given these limitations, we excluded this dataset from our classification efforts.

from scipy.signal import savgol_filter

# Preprocessing parameters
savgol_window_length = 9
savgol_polynomial_order = 3
fingerprint_region_cm1 = (300, 1800)

# Load spectrum from Renishaw glass slide control
txt_filepath = "data/Renishaw_Qontor/glass_slide_background.txt"
glass_slide_control_spectrum = RamanSpectrum.from_renishaw_txtfile(txt_filepath)

# Store preprocessed spectra
preprocessed_chlamy_spectra = []

for i, dataframe_row in chlamy_dataframe.iterrows():
    # Get the raw spectrum and its associated metadata
    raw_spectrum = chlamy_spectra[i]
    wavelength_nm = dataframe_row["λ_nm"]
    instrument = dataframe_row["instrument"]
    strain = dataframe_row["strain"]
    medium = dataframe_row["medium"]

    # Preprocessing steps
    # -------------------
    # 0.5) Apply median filtering on OpenRAMAN data to filter out hot pixels from the detector
    if instrument == "openraman":
        raw_spectrum = raw_spectrum.smooth(kernel_size=3)
    # 1) Denoising
    smoothed_intensities = savgol_filter(
        raw_spectrum.intensities,
        window_length=savgol_window_length,
        polyorder=savgol_polynomial_order,
    )
    # 1.5) Glass slide background subtraction for the Renishaw
    if instrument == "renishaw":
        smoothed_intensities -= glass_slide_control_spectrum.intensities
    # 2) Crop spectral range to the fingerprint region
    cropped_spectrum = RamanSpectrum(
        wavenumbers_cm1=raw_spectrum.wavenumbers_cm1,
        intensities=smoothed_intensities,
    ).between(*fingerprint_region_cm1)
    # 2.5) Additional cropping for 532 nm Wasatch
    if (instrument == "wasatch") and (wavelength_nm == 532):
        cropped_spectrum = cropped_spectrum.between(0, 1700)
    # 3) Baseline subtraction
    baseline_estimate, _params = aspls(cropped_spectrum.intensities, lam=lam)
    baseline_subtracted_intensities = cropped_spectrum.intensities - baseline_estimate
    # 4) Compose into RamanSpectrum and normalize
    preprocessed_spectrum = RamanSpectrum(
        wavenumbers_cm1=cropped_spectrum.wavenumbers_cm1,
        intensities=baseline_subtracted_intensities,
    ).normalize()
    preprocessed_chlamy_spectra.append(preprocessed_spectrum)

# Create figure
num_rows = len(chlamy_dataframe.groupby(["instrument", "λ_nm"]).nunique()) + 1
num_cols = 2
fig = make_subplots(
    rows=num_rows,
    cols=num_cols,
    row_heights=[0.5, 0.5, 0.1, 0.5, 0.5],  # add space to separate 532 vs 785 subplots
    shared_xaxes=True,
    shared_yaxes=True,
    horizontal_spacing=0.02,
    vertical_spacing=0.05,
    column_titles=["Raw", "Preprocessed"],
)
# Map each (instrument, wavelength) combo to a subplot
subplot_rows = {  # (instrument, λ_nm) -> row
    ("renishaw", 785): 1,
    ("wasatch", 785): 2,
    ("openraman", 532): 4,
    ("wasatch", 532): 5,
}

for (instrument, wavelength_nm), dataframe_group in chlamy_dataframe.groupby(
    ["instrument", "λ_nm"]
):
    # Only show a sample of spectra to avoid rendering issues
    num_samples = min(50, len(dataframe_group))
    for i, dataframe_row in dataframe_group.sample(num_samples).iterrows():
        # Get the raw and preprocessed spectra and its associated metadata
        raw_spectrum = chlamy_spectra[i]
        preprocessed_spectrum = preprocessed_chlamy_spectra[i]
        strain = dataframe_row["strain"]
        medium = dataframe_row["medium"]
        subplot_row = subplot_rows[(instrument, wavelength_nm)]

        # Plot each raw spectrum
        fig.add_trace(
            go.Scatter(
                x=raw_spectrum.wavenumbers_cm1,
                y=raw_spectrum.normalize().intensities,
                hoverinfo="skip",
                legendgroup=strain + medium,
                showlegend=False,
                marker={"color": color_palette[(strain, medium)]},
                opacity=0.5,
            ),
            row=subplot_row,
            col=1,
        )

        # Plot each preprocessed spectrum
        fig.add_trace(
            go.Scatter(
                x=preprocessed_spectrum.wavenumbers_cm1,
                y=preprocessed_spectrum.intensities,
                hoverinfo="skip",
                legendgroup=strain + medium,
                showlegend=False,
                marker={"color": color_palette[(strain, medium)]},
                opacity=0.5,
            ),
            row=subplot_row,
            col=2,
        )

# Plot fake traces for legend
ordered_groups = [
    ("CC-124", "TAP"),
    ("CC-125", "TAP"),
    ("CC-1373", "TAP"),
    ("CC-124", "M-N"),
    ("CC-125", "M-N"),
    ("CC-1373", "M-N"),
]
for strain, medium in ordered_groups:
    fig.add_trace(
        go.Scatter(
            x=[None],
            y=[None],
            mode="lines",
            name=f"{strain} | {medium}",
            line={"color": color_palette[(strain, medium)], "width": 3},
            legendgroup=strain + medium,
            showlegend=True,
        )
    )

# Configure plotly layout
fig.update_layout(layout, height=600)
fig.update_xaxes(title_text="Wavenumber (cm<sup>-1</sup>)", row=num_rows, col=1)
fig.update_xaxes(title_text="Wavenumber (cm<sup>-1</sup>)", row=num_rows, col=2)
for (instrument, wavelength_nm), row in subplot_rows.items():
    title = f"{instrument.capitalize()}<br>{wavelength_nm} nm"
    fig.update_yaxes(title_text=title, row=row, col=1)
fig.show()

Figure 3: Raw (left) and preprocessed (right) spectra from Chlamydomonas samples. Each row of subplots shows the spectra from a different instrument. Spectra from TAP-grown cells are shown in shades of blue; spectra from nitrogen-depleted cells are shown in shades of red. You can double-click on a particular strain, medium combination in the legend to isolate it. At most, 50 spectra from each group are shown to prevent rendering issues.

Our preprocessing pipeline worked well for the raw spectra, even though the instruments and excitation wavelengths vary greatly. For example, the fluorescence background for the 532 nm and 785 nm excitations had very different shapes but were managed well by the baseline subtraction we applied. As expected, the most significant factor driving spectral differences was the excitation wavelength, followed by which instrument we used. After controlling for these technical variables, the remaining differences in our processed spectra should mainly reveal the unique biochemical fingerprints of each algal strain and how they responded to different growth conditions.

The 532 nm spectra match well with peaks reported in some of the literature for C. reinhardtii CC-125, showing three major peaks ~1,006, ~1,156, and ~1,520 cm^–1 that have been ascribed to carotenoid pigments (Pandey et al., 2022). The 785 nm spectra generally match those from additional literature on Chlamydomonas (Samek et al., 2010), and we think a majority of the peaks are due to vibrational modes of pigments, based on studies with similarly pigmented species (Ishihara and Takahashi, 2023). Finally, we noticed that cells cultured in nitrogen-depleted media don’t have the expected lipid peaks as seen in other research studying Raman-detectable responses of other algae to nitrogen depletion (He et al., 2012). This change is also expected in C. reinhardtii (CC-125) (Msanne et al., 2012). None of our spectra had strong peaks in the 2,800–3,000 cm^–1 C-H stretch region, where lipid peaks are prominent, contrasting to the literature on other strains of Chlamydomonas (Sharma et al., 2015). The absence of these peaks may be attributed to insufficient nitrogen depletion conditions in our cultures, though further investigation is required to know for sure.

Linear discriminant analysis of Chlamydomonas spectra

Now that we have preprocessed the Raman spectra from each instrument, we can apply linear discriminant analysis (LDA) to explore how the spectra cluster based on strain differences and growth media (Figure 4). This analysis will help us identify the spectral features that best separate these groups, providing insight into how well Raman spectra can distinguish between strains and media compositions. We’ll first conduct principal component analysis (PCA) as a preprocessing step for LDA to denoise the data, reduce dimensionality, and remove multicollinearity between features. We’ve selected ten principal components, which collectively explain approximately 80% of the total variance in our dataset. We saw clear signs of overfitting when testing with a higher number of components.

from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

# Set number of PCA and LDA components
num_pca_components = 10
num_lda_components = 2

# Define groups and group labels
group_labels = {
    ("CC-124", "M-N"): 0,
    ("CC-125", "M-N"): 1,
    ("CC-1373", "M-N"): 2,
    ("CC-124", "TAP"): 3,
    ("CC-125", "TAP"): 4,
    ("CC-1373", "TAP"): 5,
}

import plotly.express as px
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

ordered_groups = [
    ("renishaw", 785),
    ("wasatch", 785),
    ("wasatch", 532),
    ("openraman", 532),
]

for instrument, wavelength_nm in ordered_groups:
    # Filter DataFrame for instrument and excitation wavelength
    grp = chlamy_dataframe.query("instrument == @instrument & λ_nm == @wavelength_nm")

    # Compose training data
    X = []
    y = []
    for (_strain, _medium), label in group_labels.items():
        indices_in_group = grp.query("strain == @_strain & medium == @_medium").index.values
        spectra_in_group = [preprocessed_chlamy_spectra[i] for i in indices_in_group]
        for spectrum in spectra_in_group:
            X.append(spectrum.intensities)
            y.append(label)
    X = np.array(X)
    y = np.array(y)

    # Construct LDA pipeline
    pipeline = Pipeline(
        [
            ("scaler", StandardScaler()),
            ("pca", PCA(n_components=num_pca_components, random_state=42)),
            ("lda", LinearDiscriminantAnalysis(n_components=num_lda_components, solver="svd")),
        ]
    )
    # Fit and transform LDA
    spectral_components = pipeline.fit(X, y).transform(X)
    # Extract explained variance ratios
    _pca_variances = pipeline.named_steps["pca"].explained_variance_ratio_
    lda_variances = pipeline.named_steps["lda"].explained_variance_ratio_

    # Add LDA components to DataFrame for plotting
    source = pd.DataFrame(
        {
            "LD-1": spectral_components[:, 0],
            "LD-2": spectral_components[:, 1],
            "label": y,
        }
    )
    group_labels_inverted = {k: v for v, k in group_labels.items()}
    source["Group"] = source["label"].map(group_labels_inverted)

    # LDA plots
    fig = px.scatter(
        source,
        x="LD-1",
        y="LD-2",
        color="Group",
        color_discrete_map=color_palette,
        marginal_x="box",
        marginal_y="box",
    )

    # Configure plotly layout
    fig.update_layout(
        layout,
        title=f"{instrument.capitalize()} | {wavelength_nm} nm",
        xaxis_title=f"LD-1 ({lda_variances[0]:.1%})",
        yaxis_title=f"LD-2 ({lda_variances[1]:.1%})",
        height=500,
        width=600,
        xaxis={"scaleanchor": "y"},
    )
    fig.show()

(a) LDA of Chlamydomonas spectra from the Renishaw (785 nm).

(b) LDA of Chlamydomonas spectra from the Wasatch (785 nm).

(c) LDA of Chlamydomonas spectra from the Wasatch (532 nm).

(d) LDA of Chlamydomonas spectra from the OpenRAMAN (532 nm).

Figure 4: Linear discriminant analysis (LDA) of Chlamydomonas spectra.

Although the LDA results exhibit substantial variance across instruments, a consistent pattern emerges in which one of the linear discriminants differentiates spectra based on growth media. This is evident in the data from the 785 nm instruments (Figure 4, A–B) in which spectra from M-N cultures predominantly cluster above the spectra from TAP-grown cells (LD-2). Interestingly, the 785 nm Wasatch data also show evidence of species-level differentiation along LD-1. While there doesn’t appear to be a clear separation between strains or species in the spectra from the 532 nm instruments (Figure 4, C–D), we do see variation from growth media in the 532 nm Wasatch (LD-1). Despite the OpenRAMAN and 532 nm Wasatch having the same excitation wavelengths, we were able to extract more spectral features from the Wasatch system (Figure 2), which may account for the differences in clustering between these two systems.

Species classification

Given that we see some species-level variation in the LDA, let’s determine how accurately we can classify each species using different machine-learning (ML) classifiers. We’ll implement and compare a few classification algorithms we expect to perform reasonably well for this task: random forest, logistic regression, and linear discriminant analysis (already tested qualitatively). We’ll include a dummy classifier as a control to validate that the models outperform random guessing. For a more reliable performance estimate, we’ll employ stratified k-fold cross-validation (with k = 5), which preserves the class distribution in each fold, ensuring our test results better reflect how these models would perform on new, unseen spectra. The results of our species classification task are shown in Figure 5.

from analysis.classification import BatchClassifier
from sklearn.dummy import DummyClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold
from tqdm.notebook import tqdm

classifiers = [
    RandomForestClassifier,
    LogisticRegression,
    LinearDiscriminantAnalysis,
    DummyClassifier,
]
species_labels = {
    "C. reinhardtii": 0,
    "C. smithii": 1,
}
num_splits = 5

classification_results_dataframe = pd.DataFrame()
for (instrument, wavelength_nm), grp in chlamy_dataframe.groupby(["instrument", "λ_nm"]):
    # Compose training data
    X = []
    y = []
    for _species, label in species_labels.items():
        indices_in_group = grp.query("species == @_species").index.values
        spectra_in_group = [preprocessed_chlamy_spectra[i] for i in indices_in_group]
        for spectrum in spectra_in_group:
            X.append(spectrum.intensities)
            y.append(label)
    X = np.array(X)
    y = np.array(y)

    # Initialize K-folds cross-validation and batch classifier
    kfolds = StratifiedKFold(n_splits=num_splits, shuffle=True, random_state=37)
    batch_classifier = BatchClassifier(classifiers=classifiers, random_state=37)

    # Train classifiers and predict
    for train, test in tqdm(kfolds.split(X, y), total=num_splits):
        X_train, X_test, y_train, y_test = X[train], X[test], y[train], y[test]
        scores = batch_classifier.fit(X_train, X_test, y_train, y_test).reset_index()
        scores["instrument"] = instrument
        scores["wavelength_nm"] = wavelength_nm
        classification_results_dataframe = pd.concat([scores, classification_results_dataframe])

# Map each (instrument, wavelength) combo to a subplot
subplot_rows = {  # (instrument, λ_nm) -> row
    ("renishaw", 785): 1,
    ("wasatch", 785): 2,
    ("openraman", 532): 3,
    ("wasatch", 532): 4,
}
subplot_titles = [f"{k[0].capitalize()} | {k[1]} nm" for k in subplot_rows.keys()]

# Create figure
num_rows = len(chlamy_dataframe.groupby(["instrument", "λ_nm"]).nunique())
fig = make_subplots(
    rows=num_rows,
    cols=1,
    shared_xaxes=True,
    subplot_titles=subplot_titles,
)

classification_results_dataframe = classification_results_dataframe.reset_index(drop=True)
for (instrument, wavelength_nm), scores in classification_results_dataframe.groupby(
    ["instrument", "wavelength_nm"]
):
    # Aggregate mean + std from each k-fold
    # (unfortunately plotly can't automatically do this from long-form DataFrames like seaborn)
    source = scores.groupby("model")["balanced_accuracy"].agg(["mean", "std"]).reset_index()
    subplot_row = subplot_rows[(instrument, wavelength_nm)]

    # Create bar plot to show classification accuracy for each instrument
    fig.add_trace(
        go.Bar(
            y=source["model"],
            x=source["mean"],
            showlegend=False,
            error_x={"array": source["std"]},
            marker={"color": color_palette[(instrument, wavelength_nm)]},
            orientation="h",
        ),
        row=subplot_row,
        col=1,
    )

# Configure plotly layout
fig.update_layout(layout, height=650)
fig.update_xaxes(title_text="Balanced accuracy", row=num_rows, tickformat=".0%")
fig.show()

Figure 5: Species classification of Chlamydomonas spectra. Error bars indicate the standard deviation in the balanced accuracy of the k = 5 folds.

Overall, the classifiers were fairly accurate in differentiating Chlamydomonas species based on spectral data. Importantly, all instruments produced spectra that yielded above-random classification performance. The 785 nm Wasatch achieved the highest classification accuracy (97%) using logistic regression, consistent with its highest-ranking clustering in the LDA analysis. Interestingly, the OpenRAMAN outperformed the 532 nm Wasatch, possibly due to its effectiveness in distinguishing C. smithii grown in M-N media, which appeared distinctly separated in its LDA projection (Figure 4, B). Model selection didn’t heavily impact outcomes, though logistic regression generally showed slight advantages across instruments — only for the OpenRAMAN did random forest classification surpass it.

Discussion

Recall that our objective wasn’t to establish which spectrometer was “best” but rather to determine which instrument(s) and acquisition parameters (particularly excitation wavelength) best aligned with our specific applications and budget. The Wasatch systems emerged as particularly well-suited to our needs, demonstrating excellent SNR characteristics, revealing numerous spectral features (especially at 785 nm), producing effective clustering in LDA, and showing superior species differentiation capabilities. The results of the multidimensional analysis are clearly sensitive to the choice of instrument and acquisition parameters, as evidenced by clustering patterns varying substantially between different instruments. This dependence highlights the importance of implementing appropriate controls and orthogonal analysis methods to validate LDA-derived results prior to drawing biological conclusions.

Beyond these performance metrics, there are a number of important considerations not captured in our analysis. Chief among them were ease of operation and automation potential — critical factors for high-throughput data collection. In this regard, the Wasatch systems also excelled, with their open-source software framework facilitating the efficient collection of large spectral datasets. We also didn’t evaluate the spectrometers’ spatial characteristics (e.g., spot size, depth of field) as these parameters were less important for our applications, and not all instruments had imaging capabilities. Had this been a greater priority, we might have leaned more towards the Renishaw, as this system demonstrated superior performance in this area. Our experiments demonstrated, however, that a probe-based Raman spectroscopy system is sufficient for a wide range of biological applications. Although cost was another factor we considered, we chose to focus the scope of this work purely on the analytical components.

Our secondary objective was to showcase a range of analytical techniques, particularly well-suited to Raman spectral datasets, and to present them in an interactive and easily reproducible format. We’ve demonstrated various approaches — from basic preprocessing to dimensionality reduction and classification methods — that can extract meaningful biological information from spectral data. We hope these techniques and the open-source package that facilitated some of the analysis, ramanalysis, prove valuable to other researchers and perhaps inspire new analytical strategies in their own work. This analysis reflects our belief that Raman spectroscopy, despite its powerful capabilities for non-destructive molecular profiling, remains underutilized as a tool for phenotyping in biological research.

References

Avasthi P, Braverman B, Bigge BM, Essock-Burns T, Lane R, Mets DG, Patton AH, York R. (2024). Raman spectroscopy enables rapid and inexpensive exploration of biology. https://doi.org/10.57844/arcadia-085e-3ecf

Braverman B, Mets DG, Sharma S. (2025). DIY Raman spectroscopy for biological research. https://doi.org/10.57844/arcadia-cd7e-443b

Fikiet MA, Tuschel D, Ermolenkov VV, Lednev IK. (2020). Clarifying glass luminescence at near-infrared excitation. https://doi.org/10.1177/0003702819879109

He X, Allen J, Black P, Baldacchini T, Huang X, Huang H, Jiang L, Lu Y. (2012). Coherent anti-stokes raman scattering and spontaneous raman spectroscopy and microscopy of microalgae with nitrogen depletion. https://doi.org/10.1364/BOE.3.002896

Ishihara J, Takahashi H. (2023). Raman spectral analysis of microbial pigment compositions in vegetative cells and heterocysts of multicellular cyanobacterium. https://doi.org/10.1016/j.bbrep.2023.101469

Msanne J, Xu D, Konda AR, Casas-Mollano JA, Awada T, Cahoon EB, Cerutti H. (2012). Metabolic and gene expression changes triggered by nitrogen deprivation in the photoautotrophically grown microalgae chlamydomonas reinhardtii and coccomyxa sp. C-169. https://doi.org/10.1016/j.phytochem.2011.12.007

Pandey S, Archana G, Bagchi D. (2022). Micro-raman spectroscopy of the light-harvesting pigments in chlamydomonas reinhardtii under salinity stress. https://doi.org/10.1016/j.saa.2022.121613

Samek O, Jonáš A, Pilát Z, Zemánek P, Nedbal L, Třı́ska J, Kotas P, Trtı́lek M. (2010). Raman microspectroscopy of individual algal cells: Sensing unsaturation of storage lipids in vivo. https://doi.org/10.3390/s100908635

Sharma SK, Nelson DR, Abdrabu R, Khraiwesh B, Jijakli K, Arnoux M, O’Connor MJ, Bahmani T, Cai H, Khapli S, others. (2015). An integrative raman microscopy-based workflow for rapid in situ analysis of microalgal lipid bodies. https://doi.org/10.1186/s13068-015-0349-1

Shimanouchi T. (1972). Tables of molecular vibrational frequencies consolidated volume i. https://doi.org/10.6028/NBS.NSRDS.39

Zhang F, Tang X, Tong A, Wang B, Wang J, Lv Y, Tang C, Wang J. (2020). Baseline correction for infrared spectra using adaptive smoothness parameter penalized least squares method. https://doi.org/10.1080/00387010.2020.1730908