Toward Gene-Correlated Spatially Resolved Metabolomics with Fingerprint Coherent Raman Imaging

Raman spectroscopy has long been known to provide sufficient information to discriminate distinct cell phenotypes. Underlying this discriminating capability is that Raman spectra provide an overall readout of the metabolic profiles that change with transcriptomic activity. Robustly associating Raman spectral changes with the regulation of specific signaling pathways may be possible, but the spectral signals of interest may be weak and vary somewhat among individuals. Establishing a Raman-to-transcriptome mapping will thus require tightly controlled and easily manipulated biological systems and high-throughput spectral acquisition. We attempt to meet these requirements using broadband coherent anti-Stokes Raman scattering (BCARS) microscopy to spatio-spectrally map the C. elegans hermaphrodite gonad in vivo at subcellular resolution. The C. elegans hermaphrodite gonad is an ideal model system with a sequential, continuous process of highly regulated spatiotemporal cellular events. We demonstrate that the BCARS spatio-spectral signatures correlate with gene expression profiles in the gonad, evincing that BCARS has potential as a spatially resolved omics surrogate.

: An example of spatial segmentation of 10 sections of the hermaphrodite gonad.  Figure S2: Clustermap generated from node frequency data from each section. Unlike in all other analyses, where the correlation was calculated between each section using node data per section (generating a 10x10 matrix), here, the correlation was calculated between each node. This 300x300 matrix was supplied to Seaborn's clustermap function to obtain this.  Figure S4: Worms imaged at multiple wavenumbers. The chosen wavenumbers specify the contrasts indicated above the images. Arrows point to mature nuclei that are not expected in the pre-loop gonad. These may be from sheath cells that were also imaged due to the finite z-plane depth of BCARS imaging. (a-c) Up and Down worms. (d-f) One worm. The scale bar (20 µm) shown in (f) applies to all images in this figure. Note that the contrast at 785 cm −1 is due to the vibration of the O-P-O phosphodiester bond. This wavenumber is the best-known marker band for DNA and RNA (nucleic acids), but phospholipids 37 and potentially other molecules that incorporate phosphate groups can also show some signal at this wavenumber. In (a), the large, more mature nuclei appear darker than the small nuclei. We note that both large and small nuclei are darker than their surrounding cytoplasm, likely due to the cytoplasm containing large amounts of RNA (and hence, O-P-O bonds). The large nuclei also incorporate more fluid, potentially as nucleoplasm, as suggested by Fig. 2(c), where the intercellular fluid branch also selects pixels within the large nuclei. Dilute fluid would contain less signal at 785 cm −1 , making these nuclei darker at this wavenumber. Figure   undergoing a presumably continuous maturation as they pass through the gonad. Therefore, as the cell state resolving power of the analysis improves, we expect the section-section correlation matrix to approach a diagonal matrix with 1 on the diagonal and -1 in other locations. (a) Directly correlating the spatially averaged mean Raman spectra across sections gives the lowest dynamic range. (b) D.S. = "Difference Spectrum". Subtracting the mean spectrum of the entire gonad from each section average improves the dynamic range. Section 7 includes sections before and after the loop, and its mean spectrum is comparable to the mean of the whole worm. When the mean is subtracted, the difference spectrum is close to 0, leading to this exceptional behaviour. See Fig. 1e. (c) The SPADE analysis shows the highest dynamic range and cell state resolving power.

cle-1
Pixel density (nor.) /Gene expression levl (nor.) Figure S8: Gene-correlated node pixels and the corresponding BCARS spectra. The scatter plots of normalized gene-correlated pixel density and gene expression level by gonad sections for (a) rme-1 and (b) cle-1. The corresponding mean BCARS spectra and differential BCARS spectra relative to the mean spectrum of the whole gonad for (c)(e) rme-1 and (d)(f) cle-1 correlated node pixels. S10 Figure S9: Coloring the Gonad SPADE trees from Fig. 3 using data from Down and One. (a) Down trees from 1-5 are shaded qualitatively similarly to those from Up, but the proportions are different, owing to a change in z-plane. Section 3 is different from all others due to the unexpected inclusion of two mature nuclei in the image (Fig. S4). Section 6 of Down is outside the image and hence excluded from the analysis. Sections 7-10 are very similar to those from Up. (b) One trees are only vaguely similar to the Up ones. While some similar nodes are highlighted, most pixels are assigned to a small number of nodes. This is characteristic of projection error, potentially due to systematic error in the retrieved spectra in the two images. S12 Figure S11: UMAP analysis of the three worms using the same procedure as SPADE. Comparison of projection error with direct BCARS spectra and difference spectra with the subtraction of an appropriate mean. Figure S11 shows the same analysis performed for the gonads of the three worms using UMAP instead of SPADE. All spectra from the gonad of the Up worm were used to train 671 two reducers: "Direct", and "Difference". For Direct, the spectra were supplied as-is. For The red, two-pronged shape in the background in both analyses is the respective training 678 set (all the spectra from the gonad of the Up worm) represented using the respective reducer.

679
The cyan pixels come from the respective gonad sections marked above each image. The UpGonadMean was subtracted from the spectra from both the Up and Down worms before 690 reduction by the Difference reducer. However, for spectra from One, the "OneGonadMean", 691 the mean of the spectra in the gonad of the One worm, was subtracted from each spectrum 692 before reduction. The One column now shows much less projection error, with pixels being 693 spread out across the set, similar to Up and Down.

694
This suggests that the systematic error across images in BCARS is small and could 695 potentially be corrected with simple modifications to the existing workflow, such as the use 696 S14 of difference spectra.
Of note is that the UMAP reduction shows a continuous and significant change between 698 every section in all three columns in (b), matching the intuition that the "cells" in these 699 sections are undergoing continuous changes. A similar result was obtained when SPADE was 700 performed with z-score standardized spectra (where the intensity at each wavenumber was 701 mean-subtracted and is divided by its variance) instead of direct spectra (data not shown).

702
This suggests that the BCARS microscopy data has much higher cell state resolving power 703 than that shown in the analysis in this work.

704
The present analysis may be justified as follows. SPADE inherently ignores data den-705 sity as information, due to the density-dependent downsampling step. UMAP inherently 706 attempts to fit a manifold to equalise the data density on it. In other words, SPADE inher-707 ently discards data density information (density of points in spectral space), while UMAP 708 exploits it. It is possible to "spread out" the data using techniques like z-score, so that 709 SPADE is more sensitive to small changes in the Raman spectrum, while such a technique 710 would have essentially no effect on UMAP.

711
However, when the z-score is applied, the branches in the SPADE tree are harder to 712 interpret, and do not provide branches corresponding strictly to the well-understood cell 713 components. This would complicate the analysis in Fig. 2. Future work can incorporate 714 techniques to resolve this, such as annotating areas generated by UMAP using known infor-715 mation about the intensity and origin of Raman peaks at different wavenumbers, or a similar 716 approach with SPADE nodes.