To examine the correspondence between Raman spectra and proteomes in E. coli, we reproduced 15 environmental conditions for which absolute quantitative proteome data are already available (Schmidt et al., 2016) and measured Raman spectra of E. coli cells under those conditions (Figure 2A and B). The culture conditions we adopted include (i) exponential growth phase in minimal media with various carbon sources, (ii) exponential growth phase in rich media, (iii) exponential growth phase with various stressors, and (iv) stationary phases (Appendix 1—table 1). We measured Raman spectra of single cells sampled from each condition and focused on the fingerprint region of biological samples, where the signals from various biomolecules concentrate (spectral range of 700–1800 cm⁻¹, Figure 2B and Appendix 1—figure 2). The Raman spectra were classified on the basis of the environmental conditions using a simple linear classifier, linear discriminant analysis (LDA) (Goodacre et al., 1998; Huang et al., 2004; De Bie et al., 2005; Figure 2C–E and Appendix 1—figure 1). This classifier calculates the most discriminatory axes by maximizing the ratio of between-condition variance to within-condition variance and reduces the dimensions of Raman data to $m-1$, where $m=15$ is the number of conditions (see ‘Experimental methods, data acquisition, and data analyses’ in Materials and methods and Section 2.1 in Appendix).

Figure 2

Download asset Open asset

The result shows that Raman spectral points from different environmental conditions are distinguishable in the $(m-1)$-dimensional LDA space (Figure 2C–E). For example, the first and second LDA axes clearly distinguish the conditions ‘LB’ and ‘stationary3days’ (Figure 2C), and the third axis distinguishes ‘Glucose42C’ and ‘GlycerolAA’ (Figure 2D). Notably, the first principal axis LDA1 correlated with growth rate significantly (Pearson correlation $r=0.81\pm 0.09$, Appendix 1—figure 2). Visualizing the Raman LDA data by embedding them on a two-dimensional plane using t-distributed stochastic neighbor embedding (t-SNE) (van der Maaten and Hinton, 2008) confirms that the points for each condition form a distinctive cluster (Figure 2F). These results imply that positions in the Raman LDA space reflect condition-dependent differences in cellular physiological states.

We next asked whether these Raman spectral differences in the LDA space could be linked to the different proteome profiles (Appendix 1—figure 1). To examine this, we hypothesized linear correspondence between the $n$-dimensional proteome column vector ${\hat{\mathit{p}}}_j$, where $n=2058$ is the number of protein species in the proteome data, and the low-dimensional ($(m-1)$-dimensional) Raman column vector ${\widehat{\mathit{𝒓}}}_j$ in condition $j$,

$B$ is an $n\times m$ matrix that connects ${\widehat{\mathit{𝒑}}}_j$ and ${\widehat{\mathit{𝒓}}}_j$. We calculated ${\widehat{\mathit{𝒓}}}_j$ as the average of the low-dimensional LDA Raman data of single cells in condition $j$ since the proteomes were measured for cell populations (Table 1).

We conducted leave-one-out cross-validation (LOOCV) to verify this linear correspondence (Figure 2G). We excluded one condition (here, $j$) as a test condition and estimated $B$ as $B_{-j}^{\mathrm{est}}$ by simple ordinary least squares (OLS) regression using the data of the rest of the conditions. We thereby estimated the proteome in condition $j$ as ${\hat{\mathit{p}}}_j^{\mathrm{e}\mathrm{s}\mathrm{t}}=B_{-j}^{\mathrm{e}\mathrm{s}\mathrm{t}}\cdot \left[\begin{array}{c}1\\ {\hat{\mathit{r}}}_j\end{array}\right]$.

The proteome profile estimated using the first four major LDA axes (LDA1–LDA4) agreed well with the actual proteome data under most conditions (Figure 2H and Appendix 1—figure 3; see ‘Raman-proteome statistical correspondence’ in Materials and methods for the estimation with all the LDA axes). Changing the condition to exclude, we estimated the proteomes for all the 15 conditions and calculated the overall estimation error by the Euclidean distance $\sum _j\Vert {\hat{\mathit{p}}}_j^{\mathrm{e}\mathrm{s}\mathrm{t}}-{\hat{\mathit{p}}}_j{\Vert }^2$. The result shows that the overall estimation error is significantly small ($p=0.00005$ by permutation test; Fisher, 1935; Pitman, 1937; Phipson and Smyth, 2010). Adopting other distance measures does not change the conclusion (Appendix 1—tables 2 and 3). These results, therefore, validate the assumption of linear correspondence between cellular Raman spectra and proteomes and confirm that condition-dependent changes in proteomes can be inferred from the corresponding low-dimensional Raman spectra.

Since the dimensionality of the proteome data is significantly higher than that of the Raman data, the result above suggests that changes in proteome profiles are constrained in low-dimensional space. The regression matrix $B$ considered above determines how the proteomes relate to the Raman LDA axes. Therefore, analyzing $B$ should provide some insights into constraints on condition-dependent changes in the proteomes (Appendix 1—figure 1).

The $n\times m$ matrix $B$ is represented as $B=\left[\begin{array}{cccc}{\mathit{b}}_0& {\mathit{b}}_1& \cdots & {\mathit{b}}_{m-1}\end{array}\right]$, where the $(k+1)$-th column ${\mathit{b}}_k={\left(b_{1k}{b}_{2k}\cdots b_{nk}\right)}^{\mathrm{\top }}$ ($0\le k\le m-1$) is the collection of coefficients of all $n$ proteins for the $k$-th LDA axis (Table 1). In the case of $k=0$, the coefficients are constant terms. We first asked whether any shared features might exist in the coefficients of $B$ depending on biological functions of corresponding proteins. We then classified the proteins according to functional annotations of Clusters of Orthologous Group (COG) classes (Tatusov et al., 1997; Tatusov et al., 2003; Galperin et al., 2015) and found that, for many proteins belonging to the ‘information storage and processing’ (ISP) COG class, the coefficients corresponding to different LDA axes are approximately proportional to the constant terms, i.e., $b_{lk}\approx c_k{b}_{l0}$, where $l$ is the index of an ISP COG class protein species and $c_k$ is the proportionality constant common to many ISP COG class protein species for the $k$-th LDA axis (Figure 3A). The ISP COG class contains various proteins involved in processing genetic information such as translation, transcription, DNA replication, and DNA repair (Schmidt et al., 2016). Simple calculations show that these proportionality relationships imply that proteins in the ISP COG class conserve their mutual abundance ratios, i.e., stoichiometry, irrespective of environmental conditions (see ‘Characterizing an SCG by analyzing the Raman-proteome correspondence matrix’ in Materials and methods).

Figure 3

Download asset Open asset

Since this is an implication from the Raman-proteome correspondence, we next examined the stoichiometry conservation only with the proteome data, evaluating the expression levels with Pearson correlation coefficients for all the pairs of the conditions for each COG class (Figure 3B). For the ISP COG class, the correlation coefficients were close to 1, whereas those for the other COG classes were significantly smaller depending on condition pairs. We also evaluated the coordination of gene expression patterns within each COG class using cosine similarity and obtained consistent results (Appendix 1—figure 4). Therefore, stoichiometry conservation is stronger in the ISP COG class than in the other COG classes. Remarkably, neither shared transcription factors nor chromosome locations can account for the observed stoichiometry conservation of many protein pairs. Indeed, although the ISP COG class shows highly coordinated expression patterns (Figure 3C) compared to the non-ISP COG class (Figure 3D), the gene loci are not chromosomally clustered in either example. Additionally, the similarity/dissimilarity of expression patterns cannot easily be inferred from transcription factor regulation patterns. These results imply multi-level regulation of their abundance.

We consulted other public quantitative proteome data of Mycobacterium tuberculosis (Schubert et al., 2015), Mycobacterium bovis (Schubert et al., 2015), and Saccharomyces cerevisiae (Lahtvee et al., 2017) under environmental perturbations and consistently found strong stoichiometry conservation of the ISP COG class (Appendix 1—figure 4). Furthermore, the same trend was observed for the genotype-dependent expression changes in E. coli proteomes (Schmidt et al., 2016; Appendix 1—figure 4).

Inspired by the existence of a large class of proteins that conserves their stoichiometry, we considered a systematic way to extract SCGs without relying on artificial functional classification of COG (Appendix 1—figure 1). Focusing only on the proteome data, we evaluated stoichiometry conservation for all the pairs of proteins in the proteome by calculating the cosine similarity of expression patterns (i.e. all $\cos {\theta }_{{\mathit{p}}_i{\mathit{p}}_j}:=\left({\mathit{p}}_i/\Vert {\mathit{p}}_i{\Vert }_2\right)\cdot \left({\mathit{p}}_j/\Vert {\mathit{p}}_j{\Vert }_2\right)$ in Figure 4A and Table 1, where each element of the $m$-dimensional vector ${\mathit{𝒑}}_i$ denotes the expression level of protein species $i$ under one of the $m$ conditions), and extracted groups in each of which the component proteins exhibit coherent expression change patterns by setting a high threshold of cosine similarity ($\ge 0.995$, Figure 4B; see ‘Direct characterization of SCGs in omics data’ in Materials and methods for details).

Figure 4

Download asset Open asset

The largest SCG (SCG 1) included many proteins in the ISP COG class (91 out of 191 SCG 1 members), such as ribosomal proteins and RNA polymerase, and also proteins in the other COG classes (Figure 4B, Appendix 1—table 4). We call this largest SCG homeostatic core, as it constitutes the largest stoichiometry-conserving unit in cells. We found that the abundance of each protein in the homeostatic core (SCG 1) increased approximately linearly with the growth rate in each condition (Figure 4C). This relationship is reminiscent of the growth law: The total ribosomal contents for translation increase linearly with growth rate (Neidhardt and Magasanik, 1960; Scott et al., 2010; Bremer and Dennis, 2008). The linear increase in the abundance of each protein in Figure 4C indicates that the growth law is valid even at the single-gene level for a large class of ribosomal and non-ribosomal proteins in the homeostatic core (Appendix 1—figure 5) (see Section 3.1 in Appendix).

Though not evenly distributed, the gene loci of the proteins in the homeostatic core are scattered throughout the chromosome (Figure 4D). Therefore, localization of gene loci to a single or a small number of operons is not likely a cause of the observed stoichiometry conservation.

The proteins in the second largest SCG (SCG 2) are expressed at high levels in the fast growth conditions, especially in the ‘LB’ condition (Figure 4C). The SCG 2 includes many proteins in the metabolism COG class (21 out of 26 SCG 2 members) (Appendix 1—table 5), and their abundance increases approximately exponentially with growth rate (Figure 4C). We also identified other condition-specific small SCGs, such as a group most expressed in the ‘GlycerolAA’ condition (SCG 3) (Appendix 1—table 6), a group mainly expressed in the ‘Fructose’ condition (SCG 4) (Appendix 1—table 7), and a group most expressed in the stationary phase conditions (SCG 5) (Appendix 1—table 8; Figure 4C).

To understand the overall strength of stoichiometry conservation of the proteins in the different SCGs, we calculated the sum of cosine similarity, $d_i=\sum _j\cos {\theta }_{{\mathit{p}}_i{\mathit{p}}_j}$, for each protein species $i$, where $\cos {\theta }_{{\mathit{p}}_i{\mathit{p}}_j}$ is cosine similarity between the $m$-dimensional expression level vectors of protein $i$ and protein $j$ (Figure 4A), and the sum is taken over all the protein species (see ‘Global proteome structures based on stoichiometric balance’ in Materials and methods). We refer to $d_i$ as ‘stoichiometry conservation centrality’ (Table 1).

The proteins in the homeostatic core had high centrality scores (Figure 5A). Therefore, these proteins tend to have more connections with other proteins in terms of stoichiometry conservation. On the other hand, the proteins in the condition-specific SCGs tend to have low centrality scores among all the proteins (Figure 5A), which suggests that their stoichiometry conservation is localized within each SCG.

Figure 5

Download asset Open asset

The stoichiometry conservation centrality is biologically relevant because it correlates with gene essentiality. Fractions of essential genes almost monotonically decrease with the ranks of centrality score (Figure 5B and Appendix 1—figure 6). We also noted that genes with high centrality scores have more orthologs determined by OrthoMCL-DB (Chen et al., 2006) across the three domains of life (Figure 5C and Appendix 1—figure 6). Likewise, genes with many orthologs tend to have higher centrality scores (Figure 5C and Appendix 1—figure 6). Therefore, the stoichiometry conservation in cells correlates with the evolutionary conservation of proteins.

To determine if the correlation of stoichiometry conservation centrality with gene essentiality and evolutionary conservation is general, we analyzed the transcriptome data from other organisms and found comparable correlations in Schizosaccharomyces pombe (Appendix 1—figure 6). In addition, we found that fractions of coding genes almost monotonically decreased with ranks of centrality score in the S. pombe data (Appendix 1—figure 6).

We further analyzed two kinds of Homo sapiens transcriptome data. One is a human cell atlas, in which expression of both coding and non-coding genes in 15 fetal organs was quantified (Cao et al., 2020), and the other is genome-wide Perturb-seq data (Replogle et al., 2022), in which genetically perturbed transcriptomes were measured mainly for coding genes. Our analysis of the human cell atlas data revealed that, while the overall distribution of stoichiometry conservation centrality was broad (Figure 5D, top), the centrality distribution of coding genes was skewed to higher values (Figure 5D, bottom) as observed for the E. coli proteome. Fractions of coding genes almost monotonically decreased with ranks of centrality (Figure 5E) as seen in the S. pombe data (Appendix 1—figure 6). Essentiality of each gene in human cells was quantified by an index called CRISPR score, which measures the fitness cost imposed by CRISPR-based inactivation of the gene (Wang et al., 2015). Genes with lower CRISPR scores are considered more essential. Our analysis revealed that genes with higher stoichiometry conservation centrality scores tend to have lower CRISPR scores, thus more essential (Figure 5F). Similarly, genes with lower CRISPR scores tend to have higher stoichiometry conservation centrality scores. Furthermore, genes with higher centrality scores have more orthologs across the three domains of life and vice versa (Figure 5G). Comparable correlations of stoichiometry conservation with essentiality and evolutionary conservation were also found in the genome-wide Perturb-seq data (Figure 5H and I). Together, these results suggest that correlations of stoichiometry conservation centrality with gene essentiality and evolutionary conservation are general and preserved from E. coli to human cells regardless of the type of perturbation (see ‘Relevance of centrality of csLE structure to biological functions’ in Materials and methods for details).

Although the previous analysis revealed the biological relevance of stoichiometry conservation centrality, it is a one-dimensional quantity and cannot capture the global architecture of omics profiles. To gain further insights into genome-wide stoichiometry-conserving relationships among genes, we next analyzed the proteomes using a method similar to Laplacian eigenmaps (LE) (Appendix 1—figure 1; Belkin and Niyogi, 2003). We consider a symmetric $n\times n$ matrix $A$ whose $(i,j)$ entry is $\cos {\theta }_{{\mathit{p}}_i{\mathit{p}}_j}$ (Figure 5J, Table 1). The entire proteome structure can be represented using the eigenvectors of normalized $A$. Major differences of this method from the ordinary LE are that we consider an edge for all node pairs and that we adopt cosine similarity for weighting edges. This method places the proteins with higher cosine similarity closer in the resulting $(m-1)$-dimensional space (see ‘Global proteome structures based on stoichiometric balance’ in Materials and methods and Section 2.1 in Appendix); we call this linear method cosine similarity LE (csLE).

In this $(m-1)$-dimensional csLE space ${\mathrm{\Omega }}_{\mathrm{L}\mathrm{E}}$, the stoichiometry conservation centrality of the proteins decreased from center to periphery (Figure 5K), which confirms that it indeed measures the extent to which each protein is close to the center in the entire stoichiometry conservation architecture. Furthermore, the proteins formed polyhedral distributions with the cluster of the proteins in the homeostatic core at the center and the clusters of the proteins in the other condition-specific SCGs at distinct vertices (Figure 5L). This distribution is consistent with the fact that the condition-specific SCGs are the components whose expression patterns are distant from the homeostatic core and also between each other.

Given that the analysis of the LDA Raman-proteome regression coefficients $B$ (Figure 3A) eventually led us to identify the stoichiometry conservation architecture in the proteome data (Figure 5), the low-dimensional proteome structure in ${\mathrm{\Omega }}_{\mathrm{LE}}$ might be related to major changes in cellular Raman spectra in the LDA space and provide insight into the Raman-proteome correspondence. To investigate this, we considered representing the proteomes on the basis of the Raman LDA axes (Appendix 1—figure 1).

The coefficients in the $n\times m$ regression matrix $B$ must satisfy the proportionality $b_{ik}/b_{i0}=b_{jk}/b_{j0}$ for all $k$-th LDA axes ($1\le k\le m-1$) for the pair of protein $i$ and protein $j$ that perfectly conserve their stoichiometry, as previously mentioned in the analysis of the ISP COG class (Figure 3A; see ‘Characterizing an SCG by analyzing the Raman-proteome correspondence matrix’ in Materials and methods and Section 2.1 in Appendix). Noting this property, we constructed another $(m-1)$-dimensional proteome space ${\mathrm{\Omega }}_{\mathrm{B}}$, assigning each protein species $i$ a coordinate $\left({\beta }_i^{\mathrm{L}\mathrm{D}\mathrm{A}1}{\beta }_i^{\mathrm{L}\mathrm{D}\mathrm{A}2}\cdots {\beta }_i^{\mathrm{L}\mathrm{D}\mathrm{A}(m-1)}\right)$, where ${\beta }_i^{\mathrm{L}\mathrm{D}\mathrm{A}k}:=b_{ik}/b_{i0}$ is the normalized coefficient of gene $i$ corresponding to the $k$-th LDA axis. As in $(m-1)$-dimensional ${\mathrm{\Omega }}_{\mathrm{LE}}$, a pair of proteins with strong stoichiometry conservation is expected to position closely in this $(m-1)$-dimensional proteome space ${\mathrm{\Omega }}_{\mathrm{B}}$. Note that the proximity of the coordinates ${\beta }_i^{\mathrm{LDA}k}$ of different proteins $i$ in ${\mathrm{\Omega }}_{\mathrm{B}}$ is equivalent to the approximate proportionality of different proteins $i$ in Figure 3A, demonstrated for the ISP COG class using the proportionality constants (normalized coefficients) $c_k$ common to different proteins.

We then found that the distribution of the proteins in ${\mathrm{\Omega }}_{\mathrm{B}}$ closely resembled the one in ${\mathrm{\Omega }}_{\mathrm{LE}}$ when visualized using the first few major axes (Figure 5L and Figure 6A). This similarity is nontrivial because ${\mathrm{\Omega }}_{\mathrm{LE}}$ is constructed only from the proteome data, whereas ${\mathrm{\Omega }}_{\mathrm{B}}$ depends on the $(m-1)$-dimensional Raman LDA space (Figure 2C–E).

Figure 6

Download asset Open asset

We remark that each axis of ${\mathrm{\Omega }}_{\mathrm{B}}$ is directly linked to the corresponding Raman LDA axis. Consequently, the orthants in ${\mathrm{\Omega }}_{\mathrm{B}}$ where the condition-specific protein species reside agree with those in the Raman LDA space where the cellular Raman spectra under corresponding conditions reside (Appendix 1—figure 10) (see ‘Global omics structures characterized by Raman-omics correspondences’ in Materials and methods and Section 2.1 in Appendix). Indeed, we find such orthant agreement between the proteins in the condition-specific SCGs (SCG 2–SCG 5) and the cellular Raman spectra under the corresponding conditions (Figure 6B and C). This straightforward correspondence between ${\mathrm{\Omega }}_{\mathrm{B}}$ and the Raman LDA space allows us to examine the relationship between changes in cellular Raman spectra and omics components’ stoichiometry conservation architecture by comparing the two proteome structures in ${\mathrm{\Omega }}_{\mathrm{B}}$ and ${\mathrm{\Omega }}_{\mathrm{LE}}$.

To understand rigorously what the similarity of the proteome structures in ${\mathrm{\Omega }}_{\mathrm{B}}$ and ${\mathrm{\Omega }}_{\mathrm{LE}}$ signifies (Figure 6C and D), we clarified the mathematical relation between the coordinates of the proteins in these two spaces (Figure 6E and Appendix 1—figure 1; see Sections 2.1 and 2.2 in Appendix for details). We then characterized the two mathematical conditions that must be satisfied simultaneously (Figure 6E).

The first condition is that major axes of the Raman LDA space and those of the proteome csLE space correspond (Figure 6E). Consequently, cellular Raman spectra under a condition accompanying the expression of a condition-specific SCG must be significantly different from those under conditions with the expression of other condition-specific SCGs in a manner distinguishable by LDA. Mathematically, this condition is related to the $m\times m$ orthogonal matrix $\mathrm{\Theta }$ that appears in the equation in Figure 6E. For the distributions of the proteome components to be similar in the low-dimensional subspaces of ${\mathrm{\Omega }}_{\mathrm{LE}}$ and ${\mathrm{\Omega }}_{\mathrm{B}}$, $\mathrm{\Theta }$ must be close to the identity matrix with small off-diagonal elements (Figure 6E). We verified this first condition with the data (Appendix 1—figure 9; see ‘Evaluating similarity between orthogonal matrix $\mathrm{\Theta }$ and identity matrix’ in Materials and methods for details).

The second condition relates to the proportionality of the $n$-dimensional vectors ${\mathit{𝒃}}_0$ and ${\mathit{𝒃}}_0^{\mathrm{est}}$ in Figure 6E. This proportionality relation can be transformed into another relation that $d_i$ is proportional to $g_i:=\Vert {\mathit{p}}_i{\Vert }_1/\Vert {\mathit{p}}_i{\Vert }_2$, where ${\parallel {\mathit{𝒑}}_i\parallel }_1$ and ${\parallel {\mathit{𝒑}}_i\parallel }_2$ are the $L^1$ and $L^2$ norms of the expression-level $m$-dimensional vector of protein $i$ across conditions (Figures 4A and 6E, Table 1).

$g_i$ can be interpreted as the expression generality score. When $g_i$ is large, the protein $i$ is expressed generally across conditions; when $g_i$ is small, this is expressed only under specific conditions (Appendix 1—figure 8) (see ‘Interpretation of $L^1$ norm/$L^2$ norm ratio of an expression vector as a quantitative measure of expression generality’ in Materials and methods). Therefore, the proportionality between $d_i$ and $g_i$ indicates that the proteins with high stoichiometry conservation centrality must be expressed nonspecifically to conditions. We also verified this condition with the data, confirming that it is indeed satisfied (Figure 7A and Appendix 1—figure 9).

Figure 7

Download asset Open asset

The spread of the points from the proportionality diagonal line of the E. coli proteome data in Figure 7A was found related to the growth rate under the condition where each protein is expressed the most (see Section 2.2 in Appendix for a detailed analysis on the origin of the deviation). Consequently, one can envisage a growth-rate-dependent expression pattern of each protein on the basis of its relative position in this $g_i$-$d_i$ plot (Figure 7B and C). For example, both BamB and YqjD are expressed nonspecifically to the conditions with nearly identical expression generality scores. However, BamB is expressed at higher levels under fast growth conditions, whereas YqjP is expressed at higher levels under slow growth conditions due to their relative positions to the proportionality line. A similar growth rate dependence is observed for PaaE and DgoA, but with more prominent condition specificity because these proteins are characterized by their low expression generality scores. These growth-rate-dependent deviation patterns might hint at a new growth law that governs the total relative expression changes of the proteome components (see Section 2.2 in Appendix for detailed discussion).

We also examined the generality of the aforementioned two conditions using the Raman and proteome data of E. coli strains with different genotypes (BW25113, MG1655, and NCM3722) under two culture conditions (Schmidt et al., 2016) and the Raman and transcriptome data of S. pombe under 10 culture conditions (Kobayashi-Kirschvink et al., 2018). Applying csLE to the omics data, we again found similar omics structures between ${\mathrm{\Omega }}_{\mathrm{LE}}$ and ${\mathrm{\Omega }}_{\mathrm{B}}$ when visualized using the first few major axes, with homeostatic cores at the centers and condition-specific SCGs at the vertices (Appendix 1—figures 11 and 12).

Proportionality between stoichiometry conservation centrality and expression generality score was also confirmed in both additional datasets (Appendix 1—figure 7). We further used publicly available quantitative proteome data of M. tuberculosis, M. bovis, and S. cerevisiae (Schubert et al., 2015; Lahtvee et al., 2017) to examine this relation and confirmed that the proportionality universally holds (Appendix 1—figures 7 and 13). Almost no deviation from the proportionality line existed in the S. cerevisiae proteome data measured for the cells in different media but cultured in chemostats with an identical dilution rate (thus, identical growth rate), which is consistent with the result of E. coli in which the deviations were related to the growth rate differences.

Cells were washed three times with 0.9% aqueous solution of NaCl, and 5 µL of the suspension was placed on a synthetic quartz slide glass (Toshin Riko Co., Ltd.) and dried. Raman spectra of cells were measured with a Raman microscope (Appendix 1—figure 2), where a custom-built Raman system (STR-Raman, AIRIX) was integrated into a microscope (Ti-E, Nikon). Excitation light was generated by a $532\mathrm{nm}$ continuous-wave diode-pumped solid-state laser (Gem 532, Laser Quantum). We altered the first version of this Raman microscope (Kobayashi-Kirschvink et al., 2018), and light from the laser oscillator was transmitted by mirrors in this research. A $100\times$ and $\mathrm{NA}=0.9$ air objective lens (MPLN100X, Olympus) was used. Raman scattered light was collected by an optical fiber and transmitted to a spectrometer (Acton SP2300i, Princeton Instruments). Dispersed light by a $300\mathrm{gr}/\mathrm{mm}$ grating was projected onto an image sensor of an sCMOS camera (OrcaFlash 4.0 v2, Hamamatsu Photonics). The sCMOS camera was water-cooled at 15°C to reduce dark noise. The exposure time for each cell was $10\mathrm{s}$. Randomly selected 15 cells were measured per condition per replicate. Raman spectrum of background was measured for each cell with $10\mathrm{s}$ exposure in an area close to a targeted cell where neither cells nor NaCl crystals existed.

In our setup, the laser power at the sample stage was $21\mathrm{mW}$. The measurement system and processes were controlled using Micro-Manager 1.4 (Edelstein et al., 2014) and a plugin we made.

Readout noise of sCMOS image sensors is pixel-dependent. A noise reduction filter developed in Kobayashi-Kirschvink et al., 2018, on the basis of Huang et al., 2013, was applied to measured spectral images by using 10,000 blank images obtained with the same sCMOS sensor with exposure time of $10\mathrm{s}$. See Kobayashi-Kirschvink et al., 2018 for details.

After noise reduction with the filter, pixel counts were summed up along the direction perpendicular to wavenumber. A background spectrum was subtracted from a cellular Raman spectrum. A pixel region corresponding to the range from $632{\mathrm{cm}}^{-1}$ to $1862{\mathrm{cm}}^{-1}$ was cropped. The cropped spectrum was smoothed with a Savitzky-Golay filter (Savitzky and Golay, 1964). To minimize the effect of laser excitation variations, each spectrum was normalized by subtracting the average and dividing it by the standard deviation.

We wrote scripts and analyzed data using MATLAB (R2019a and R2023b), except for Brunner-Munzel test, for which we used R (version 4.0.3) (see ‘Centrality-evolutionary conservation correlation’ in Materials and methods).

Related to Figure 2 in the main text, we first performed LDA against the Raman data. LDA is a linear classifier; it finds the most discriminatory bases by maximizing the ratio of the between-class variance to the within-class variance and reduces the dimensions of the data to $m-1$, where $m$ is the number of classes (Huang et al., 2004; De Bie et al., 2005; Goodacre et al., 1998). In the case of our main data, classes are culture conditions. In the verification step of the correspondence between the LDA Raman and omics data, we conducted LOOCV. In LOOCV, one condition is used as test data and the remaining conditions are used as training data. This is repeated by changing the condition to exclude.

The details of the data analyses are provided in the sections below.

In the previous section, we identified SCGs by setting a threshold of cosine similarity for extracting protein pairs. We next removed the threshold and considered the ‘distance’ with respect to cosine similarity for all the protein pairs to capture the global proteome structure that includes SCGs.

The cosine similarity for all the $\left(\genfrac{}{}{0pt}{}{n}{2}\right)$ pairs of proteins can be summarized in one matrix as

where ($i,j$) component represents cosine similarity between proteins $i$ and $j$. Assuming that this matrix is an adjacency matrix in graph theory and network theory, the entire proteomes are considered as a weighted undirected complete graph (with loops), where nodes correspond to protein types and any protein pair is connected by an undirected edge. Each edge is weighted by the cosine similarity between the two protein species at both ends. Note that all the diagonal elements of $A$ are one, which represents that each node has a loop with weight of one. These were introduced just for simplicity.

To ask whether the SCGs identified in the previous section have any unique features in this network, we evaluated the degree to which each node is central in the network structure. In graph theory, ‘centrality’ is known as an index to measure how ‘important’ or ‘influential’ each node is. In particular, we employed a measure called ‘degree centrality’ (for weighted graphs) (Nieminen, 1973; Segarra and Ribeiro, 2015). Degree centrality, which is also called ‘degree’, simply measures ‘influence’ of a node on a network on the basis of links with its direct neighborhood. One can obtain a degree centrality value by calculating the sum of the weights of all the edges connected to each node (see also the definition of the degree matrix $D$ in Equation 14 below). We note that in our graph, degree centrality vector $A{\mathrm{𝟏}}_n(=D{\mathrm{𝟏}}_n)$, where ${\mathrm{𝟏}}_n$ is an $n$-dimensional column vector of which all elements are one, is equal to the eigenvector corresponding to the largest eigenvalue of a ‘normalized’ adjacency matrix ${(D^{-1}A)}^{\top }=A{D}^{-1}$ up to multiplication by a constant. From this perspective, the centrality index we adopted measures ‘influence’ of a node in a recursive manner depending on ‘influence’ of its neighboring nodes. A well-known example of this centrality indicator is Google’s PageRank (Brin and Page, 1998) used for ranking web pages on the World Wide Web. It can be regarded as a variant of ‘eigenvector centrality’ (the eigenvector corresponding to the largest eigenvalue of the adjacency matrix $A$) (Bonacich, 1972; Segarra and Ribeiro, 2015). As explained in the main text, the protein species in the homeostatic core (the largest SCG) had high centrality scores, while those in the other condition-specific SCGs had low centrality scores (Figure 5A).

We directly observed the global stoichiometry conservation structure of this proteome graph using Laplacian eigenmaps (Figures 5K, L, 6D). In general, a graph can be uniquely specified not only by the adjacency matrix $A$, but also by the Laplacian matrix $L$ defined as

where $D=\left(d_{ij}\right)$ is the degree matrix with the components of

where ${\mathrm{𝟏}}_n$ is an $n$-dimensional column vector of which all elements are one. The $(i,i)$-element of $D$ represents the sum of the weights of all the edges connected to node $i$. In our case, it represents the sum of cosine similarity values between protein $i$ and the other proteins. To see the entire proteome graph structure, we specifically employed the normalized Laplacian,

We remark that there are two types of often-used normalized Laplacian matrices, $L_{\mathrm{rw}}$ and $L_{\mathrm{sym}}=D^{-1/2}L{D}^{-1/2}=I-D^{-1/2}A{D}^{-1/2}$, in the field of machine learning (von Luxburg, 2007), and our mathematical analysis can provide a clear interpretation to each of them in the context of the Raman-proteome linear correspondence as described in Appendix 2.1.5.

There exist $m-1$ nontrivial eigenvalues of $L_{\mathrm{rw}}$ that are greater than zero and less than one. We write these $m-1$ eigenvalues as ${\lambda }_{\mathrm{LE1}},\mathrm{\dots },{\lambda }_{\mathrm{LE}(m-1)}$ from the smallest and the corresponding eigenvectors as ${\mathit{𝒗}}_{\mathrm{rw},1},\mathrm{\dots },{\mathit{𝒗}}_{\mathrm{rw},(m-1)}$. Additionally, we denote the eigenvector corresponding to the eigenvalue zero as ${\mathit{𝒗}}_{\mathrm{rw},0}$. Using these eigenvectors, one can construct a matrix $\begin{array}{l}{\tilde{V}}_{\mathrm{r}\mathrm{w}}=\left[\begin{array}{cccc}{\mathit{v}}_{\mathrm{r}\mathrm{w},0}& {\mathit{v}}_{\mathrm{r}\mathrm{w},1}& \cdots & {\mathit{v}}_{\mathrm{r}\mathrm{w},(m-1)}\end{array}\right]\end{array}$ and visualize a proteome, assigning protein $j$ with a coordinate specified by the elements after the second column in the $j$-th row of ${\tilde{V}}_{\mathrm{r}\mathrm{w}}$, i.e., by the $j$-th row of $\left[\begin{array}{ccc}{\mathit{v}}_{\mathrm{r}\mathrm{w},1}& \cdots & {\mathit{v}}_{\mathrm{r}\mathrm{w},(m-1)}\end{array}\right]$. The csLE structure we illustrate in the main figures was produced by using these eigenvectors. For example, the csLE1-csLE2 figure in the main text (Figure 6D) is a scatterplot between ${\mathit{𝒗}}_{\mathrm{rw},1}$ and ${\mathit{𝒗}}_{\mathrm{rw},2}$. Note that the closer to one the cosine similarity of a protein pair is (the more similar their expression patterns are), the ‘closer’ the two protein species are placed (see Section 2.1.5 in Appendix for details).

This method of obtaining low-dimensional representation of data using eigenvectors of a graph Laplacian is known as Laplacian eigenmaps (LE) (Belkin and Niyogi, 2001; Belkin and Niyogi, 2003). Thus, what we explained above is the LE of a graph with edges weighted with cosine similarity of expression patterns of nodes (protein species). It differs from the original and common usages of LE in that the graph we considered is a complete graph (with loops) and that the weight of edges (pairwise similarity of nodes) is cosine similarity. It has made all the mathematical formulations linear, which allowed us to biologically interpret the results with mathematically rigorous analyses. We also remark that our graph representation of proteome does not rely on existing knowledge on the underlying interaction and regulatory networks of proteins and is based only on final expression levels of the proteins. Therefore, the results are robust against the uncertainty of underlying molecular detail.

As we see in Appendix 2.1.5, an orthogonal matrix $\mathrm{\Theta }$ that appears in the relation connecting the two types of proteome structure must be close to an identity matrix to guarantee the structural similarity. To evaluate to what extent $\mathrm{\Theta }$ is close to an identity matrix, we generated many random orthogonal matrices (Mezzadri, 2006) and compared $\mathrm{\Theta }$ and the identity matrix with them.

We first multiplied each orthogonal matrix by itself in the sense of Hadamard product (element-wise product). Then, we regarded the resultant matrix as a scatterplot (Appendix 1—figure 9B) and calculated its Pearson correlation coefficient, assuming that $(i,j)$ element was the frequency of ‘data points’ at the coordinate $(i,j)$. The obtained Pearson correlation coefficient can be regarded as a measure of closeness to the identity matrix. In the case of the identity matrix, the correlation coefficient takes the maximum value, one, because non-zero values are concentrated on the diagonal part.

We calculated the square of each matrix in the sense of Hadamard product for two reasons. First, since all elements of the resultant matrix are non-negative, one can ensure that the number of ‘points’ is non-negative at any coordinate. Note that it is not necessarily an integer here. Second, the sum of all the elements of the resultant matrix is necessarily $m$. Thus, the total number of ‘points’ is equally $m$ for any $m\times m$ orthogonal matrices compared.

In addition to this method, we also evaluated the closeness of $\mathrm{\Theta }$ to the identity matrix (i) by comparing the magnitudes of off-diagonal elements among $\mathrm{\Theta }$, the identity matrix, and random orthogonal matrices, and (ii) by comparing the magnitudes of elements of leading principal submatrices among $\mathrm{\Theta }$, the identity matrix, and random orthogonal matrices. In (i), from a part consisting of $(m-1)$- and $-(m-1)$-diagonals ($(1,m)$ and $(m,1)$ elements) to the whole matrix, we expand step by step the area to consider by including $i$- and $-i$-diagonals ($m-1\ge i\ge 0$, the final step is inclusion of the main diagonal), and calculated the sum of the square of the elements in the area at each step. In (ii), from the smallest leading principal submatrix ($(1,1)$ element) to the whole matrix, we expand the area to consider step by step and calculate the sum of the square of the elements in the area at each step. See also schematic diagrams in the figures in Appendix (e.g. Appendix 1—figure 9D and E).

In Appendix 2.1.5, we will also see that even if $\mathrm{\Theta }$ is close to the identity matrix, there is another condition which must be met to guarantee the similarity of the two types of proteome structure. By considering the mathematics behind the condition, we will reveal that the two indices, stoichiometry conservation centrality (degree) $d_j={\sum }_i\cos {\theta }_{{\mathit{𝒑}}_i{\mathit{𝒑}}_j}$ and expression generality score $g_j={\Vert {\mathit{p}}_j\Vert }_1/{\Vert {\mathit{p}}_j\Vert }_2$, must be mutually proportional. Here, we explain why $g_j$ is a quantitative measure of the generality (or constancy) of expression levels.

First, we note that the ratio ${\Vert {\mathit{p}}_j\Vert }_1/{\Vert {\mathit{p}}_j\Vert }_2$ is independent of the magnitude of the expression vector ${\mathit{𝒑}}_j$. In other words, normalization does not affect the ratio:

On the basis of this, we only consider normalized expression vectors ${\mathit{p}}_j/{\Vert {\mathit{p}}_j\Vert }_2$ without loss of generality.

By definition, $L^2$ norm of a normalized expression vector (the denominator of the most left-hand side of Equation 17) equals one. Thus, the ratio we are considering equals the $L^1$ norm of the normalized expression vector:

Here, we write

where ${\tilde{p}}_{1j},{\tilde{p}}_{2j},\dots ,{\tilde{p}}_{mj}\ge 0$. Then,

Note that ${\sum }_{i=1}^m{\left({\tilde{p}}_{ij}\right)}^2=1$ holds because of normalization. Therefore, any normalized expression vector corresponds to a point on the first orthant division of the unit $(m-1)$-sphere $\sum _{i=1}^m{\left(x_i\right)}^2=1(x_1,\dots ,x_m\ge 0)$.

Next, we consider a hyperplane ${\sum }_{i=1}^m{x}_i=k$ which passes through the point in Equation 19. Since all the coefficients of this hyperplane are equal, all the $m$ intercepts are also equal. The intercept value is $k={\sum }_{i=1}^m{\tilde{p}}_{ij}$, thus equals the ratio in Equation 20. In other words, the ratio from Equation 20 appears as an intercept of the hyperplane passing through the point corresponding to the normalized vector ${\mathit{p}}_j/{\Vert {\mathit{p}}_j\Vert }_2$ with all the coefficients equal to one.

By simple calculation, one can see that the two surfaces ${\sum }_{i=1}^m{\left(x_i\right)}^2=1(x_1,\mathrm{\dots },x_m\ge 0)$ and $\sum _{i=1}^m{x}_i=k$ intersect when $1\le k\le \sqrt{m}$. (In other words, ${\Vert {\mathit{p}}_j\Vert }_2\le {\Vert {\mathit{p}}_j\Vert }_1\le \sqrt{m}{\Vert {\mathit{p}}_j\Vert }_2$ holds.) The intercept value $k={\sum }_{i=1}^m{\tilde{p}}_{ij}$ takes the maximum $k=\sqrt{m}$ when the normalized expression vector points to the ‘center’ of the first orthant division of the unit $(m-1)$-sphere, i.e., when

This means that the expression level is even and constant across the conditions. When this evenness of expression level breaks, the intercept value $k$ decreases, and it attains the minimum $k=1$ when the normalized expression vector overlaps with an axis, i.e., when

which corresponds to a completely ‘condition-specific expression pattern’ (μ is the condition’s index).

See Appendix 1—figure 8A and B for a graphical explanation of the argument for the two- and three-dimensional cases.

As mentioned in the main text, we confirmed that PCA could find a proteome structure (Appendix 1—figure 6H) similar to the csLE structure (Figure 5L and Figure 6D). Since cosine similarity of expression vectors is inner product of the $L^2$-normalized expression vectors, we also performed $L^2$ normalization of proteome data before applying PCA in this analysis. In other words, we applied PCA to a normalized proteome data $\left[{\mathit{p}}_1/{\Vert {\mathit{p}}_1\Vert }_2\cdots {\mathit{p}}_n/{\Vert {\mathit{p}}_n\Vert }_2\right]$ (see Appendix 2.1.5).

We remark that, despite the structural similarity between the PCA structure and the csLE structure, csLE has an advantage over PCA in that the relative proximity of positions reflects the strength of stoichiometry conservation between each element. In addition, as shown in the main text and Section 2 in Appendix, csLE of omics data has a direct quantitative connection to cellular Raman spectra, which is not the case for PCA.

See Materials and methods section in the main text.

To clarify what is nontrivial in the correspondence between LDA Raman and csLE proteome, here we derive rigorous mathematical relations for the correspondence through linear algebraic calculation.