Raman measurements. All Raman measurements were performed using the home-built confocal Raman microscope used in our previous studies 20, and details are reported in the Supplementary Information. In brief, the Raman setup uses a Kr ion laser at the wavelength of 647.1 nm coupled with a modified microscope (BX41, Olympus, Japan) equipped with a dipping water objective (63x/1.0 NA). The Raman photons are confocally detected by a home-built spectrograph equipped with a EMCCD (1600×200 pixels). For cell measurements, Raman maps containing 64×64 (4096) spectra were collected by scanning the entire cell area and each point spectrum was measured using 100 ms acquisition time. The optical lateral resolution was around 355 nm FWHM and the axial resolution was around 1270 nm. The laser power measured on the sample stage was set to 35 mW. Pre-processing. Single spectra were first independently corrected using custom software (LabVIEW, National Instruments Corp., TX) as follows: a) cosmic rays removal by singular value decomposition (SVD) ; b) camera offset subtraction; c) CCD response correction (intensity and etaloning) using a tungsten halogen light with known emission (Avalight-HAL, Avantes BV, NL)); d) wavenumber calibration using the zero-wavenumber laser line, toluene and argonmercury emission (CAL-2000, Ocean Optics, Germany); e) denoising by SVD . Before global cluster analysis a fully automated pre-processing pipeline was written in Python (NumPy, SciPy, pandas, scikit-learn, numba, Matplotlib, and seaborn modules) and applied to the entire hyperspectral dataset of the 319 cells using identical parameters and including the following nine steps (see Fig. ): a) preliminary spectral truncation (550 and 3100 cm -1 ); b) map truncation (62×64 (3968) pixels) to remove the first two scanlines affected by high autofluorescence; c) spatial median filter (window size 3x3 pixels) to reduce the signal heterogeneity between neighboring pixels; d) background identification by a 4-level kmeans clustering step (performed in the 2800-3030 cm -1 region) identifying as background the cluster with the lowest mean intensity. The spectrum of each cell pixel was then corrected by subtracting the median spectrum of the background. Pixel belonging to the background (>600.00 spectra) were then excluded from the final dataset. Also small (<64 pixels) isolated connected regions of foreground pixels were removed, and morphological erosion (with a 1-pixel disc) was performed on the foreground pixels to optimize cell image quality; e) Whittaker-Eilers smoother at pixel level (secondorder penalties, λ=100); f) Simultaneously with the previous step, each spectrum was put on an equidistant wavenumber grid (step size of 2 cm -1 ) over the entire spectral region; g) final spectral truncation (600-1800; 2800-3030 cm -1 ); h) baseline correction using the rubberband method; i) unit vector normalization (L 2 norm). Before the global clustering analysis, four cells were excluded from the subsequent analysis (see Table ).