Proteomic indicators of oxidation and hydration state in colorectal cancer

Data sources

This section describes the data sources and additional data processing steps applied in this study. An attempt was made to locate all currently available proteomic studies for clinical tissue on CRC including, among others, those listed in the “Tissue” and “Tissue subproteomes” sections of the review paper by De Wit et al. (2013) and in Supporting Table 3 (“Clinical Samples”) of the review paper by Martínez-Aguilar et al. (2013). To make the comparisons more robust, only datasets with at least 30 proteins in each of the up- and down-regulated groups were considered; however, all datasets from a given study were included if at least one of the datasets met this criterion. The reference keys for the selected studies shown below and in Table 1 are derived from the names of the authors and year of publication.

In comparisons between groups of up- and down-expressed proteins, the convention in this study is to consider proteins with higher expression in normal tissue or less-advanced cancer stages as a “normal” group (group 1), with number of proteins n₁, while proteins with higher expression in cancer or more-advanced cancer stages are categorized as a “cancer” group (group 2), with number of proteins n₂. For example, in the dataset of Uzozie et al. (2014) comparing normal mucosa and adenoma, the proteins up-expressed in adenoma are assigned to group 2, while in the adenoma—carcinoma dataset of Knol et al. (2014), the proteins with higher expression in adenoma are assigned to group 1 (see Table 1).

Names or IDs of genes or proteins given in the sources were searched in UniProt (The UniProt Consortium, 2015). The corresponding UniProt IDs are provided in the *.csv data files in Dataset S1. Amino acid sequences of human proteins were taken from the UniProt reference proteome (files UP000005640_9606.fasta.gz containing canonical, manually reviewed sequences, and UP000005640_9606_additional.fasta.gz containing isoforms and unreviewed sequences, dated 2016-04-13, downloaded from ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/reference_proteomes/Eukaryota/). Entire sequences were used; i.e., signal peptides and propeptides were not removed when calculating the amino acid compositions. However, amino acid compositions were calculated for particular isoforms, if these were identified in the sources. Files human.aa.csv and human_additional.aa.csv in Dataset S1 contain the amino acid compositions of the proteins calculated from the UniProt reference proteome. In a few cases, amino acid compositions of unreviewed or obsolete sequences in UniProt, not available in the reference proteome, were individually compiled; these are contained in file human2.aa.csv in Dataset S1.

Reported gene names were converted to UniProt IDs using the UniProt mapping tool (http://www.uniprot.org/mapping), and IPI accession numbers were converted to UniProt IDs using the DAVID conversion tool (https://david.ncifcrf.gov/content.jsp?file=conversion.html). For proteins with no automatically generated matches, manual searches in UniProt of the protein descriptions, where available, were performed. Proteins with missing or duplicated identifiers, or those that could not be matched to a UniProt ID, were omitted from the comparisons here. Therefore, the numbers of proteins actually used in the comparisons (listed in Table 1) may be different from the numbers of proteins reported by the authors and summarized below.

WTK+08: Watanabe et al. (2008) used 2-nitrobenzenesulfenyl labeling and MS/MS analysis to identify 128 proteins with differential expression in paired CRC and normal tissue specimens from 12 patients. The list of proteins used in this study was generated by combining the lists of up- and down-regulated proteins from Table 1 and Supplementary Data 1 of Watanabe et al. (2008) with the Swiss-Prot and UniProt accession numbers from their Supplementary Data 2.

AKP+10: Albrethsen et al. (2010) used nano-LC-MS/MS to characterize proteins from the nuclear matrix fraction in samples from 2 patients each with adenoma (ADE), chromosomal instability CRC (CIN+) and microsatellite instability CRC (MIN+). Cluster analysis was used to classify proteins with differential expression between ADE and CIN+, MIN+, or in both subtypes of carcinoma (CRC). Here, gene names from Supplementary Tables 5–7 of Albrethsen et al. (2010) were converted to UniProt IDs using the UniProt mapping tool.

JKMF10: Jimenez et al. (2010) compiled a list of candidate serum biomarkers from a meta-analysis of the literature. In the meta-analysis, 99 up- or down-expressed proteins were identified in at least 2 studies. The list of UniProt IDs used in this study was taken from Table 4 of Jimenez et al. (2010).

XZC+10: Xie et al. (2010) used a gel-enhanced LC-MS method to analyze proteins in pooled tissue samples from 13 stage I and 24 stage II CRC patients and pooled normal colonic tissues from the same patients. Here, IPI accession numbers from Supplemental Table 4 of Xie et al. (2010) were converted to UniProt IDs using the DAVID conversion tool.

ZYS+10: Zhang et al. (2010) used acetylation stable isotope labeling and LTQ-FT MS to analyze proteins in pooled microdissected epithelial samples of tumor and normal mucosa from 20 patients, finding 67 and 70 proteins with increased or decreased expression (ratios ≥2 or ≤0.5). Here, IPI accession numbers from Supplemental Table 4 of Zhang et al. (2010) were converted to UniProt IDs using the DAVID conversion tool.

BPV+11: Besson et al. (2011) analyzed microdissected cancer and normal tissues from 28 patients (4 adenoma samples and 24 CRC samples at different stages) using iTRAQ labeling and MALDI-TOF/TOF MS to identify 555 proteins with differential expression between adenoma and stage I, II, III, IV CRC. Here, gene names from supplemental Table 9 of Besson et al. (2011) were converted to UniProt IDs using the UniProt mapping tool.

JCF+11: Jankova et al. (2011) analyzed paired samples from 16 patients using iTRAQ-MS to identify 118 proteins with >1.3-fold differential expression between CRC tumors and adjacent normal mucosa. The protein list used in this study was taken from Supplementary Table 2 of Jankova et al. (2011).

MRK+11: Mikula et al. (2011) used iTRAQ labeling with LC-MS/MS to identify a total of 1,061 proteins with differential expression (fold change ≥1.5 and false discovery rate ≤0.01) between pooled samples of 4 normal colon (NC), 12 tubular or tubulo-villous adenoma (AD) and 5 adenocarcinoma (AC) tissues. The list of proteins used in this study was taken from Table S8 of Mikula et al. (2011).

KKL+12: Kim et al. (2012) used difference in-gel electrophoresis (DIGE) and cleavable isotope-coded affinity tag (cICAT) labeling followed by mass spectrometry to identify 175 proteins with more than 2-fold abundance ratios between microdissected and pooled tumor tissues from stage-IV CRC patients with good outcomes (survived more than five years; 3 patients) and poor outcomes (died within 25 months; 3 patients). The protein list used in this study was made by filtering the cICAT data from Supplementary Table 5 of Kim et al. (2012) with an abundance ratio cutoff of >2 or <0.5, giving 147 proteins. IPI accession numbers were converted to UniProt IDs using the DAVID conversion tool.

KYK+12: Kang et al. (2012) used mTRAQ and cICAT analysis of pooled microsatellite stable (MSS-type) CRC tissues and pooled matched normal tissues from 3 patients to identify 1,009 and 478 proteins in cancer tissue with increased or decreased expression by higher than 2-fold, respectively. Here, the list of proteins from Supplementary Table 4 of Kang et al. (2012) was filtered to include proteins with expression ratio >2 or <0.5 in both mTRAQ and cICAT analyses, leaving 175 up-expressed and 248 down-expressed proteins in CRC. Gene names were converted to UniProt IDs using the UniProt mapping tool.

WOD+12: Wiśniewski et al. (2012) used LC-MS/MS to analyze proteins in microdissected samples of formalin-fixed paraffin-embedded (FFPE) tissue from 8 patients; at P < 0.01, 762 proteins had differential expression between normal mucosa and primary tumors. The list of proteins used in this study was taken from Supplementary Table 4 of Wiśniewski et al. (2012).

YLZ+12: Yao et al. (2012) analyzed the conditioned media of paired stage I or IIA CRC and normal tissues from 9 patients using lectin affinity capture for glycoprotein (secreted protein) enrichment by nano LC-MS/MS to identify 68 up-regulated and 55 down-regulated differentially expressed proteins. IPI accession numbers listed in Supplementary Table 2 of Yao et al. (2012) were converted to UniProt IDs using the DAVID conversion tool.

MCZ+13: Mu et al. (2013) used laser capture microdissection (LCM) to separate stromal cells from 8 colon adenocarcinoma and 8 non-neoplastic tissue samples, which were pooled and analyzed by iTRAQ to identify 70 differentially expressed proteins. Here, gi numbers listed in Table 1 of Mu et al. (2013) were converted to UniProt IDs using the UniProt mapping tool; FASTA sequences of 31 proteins not found in UniProt were downloaded from NCBI and amino acid compositions were added to human2.aa.csv.

KWA+14: Knol et al. (2014) used differential biochemical extraction to isolate the chromatin-binding fraction in frozen samples of colon adenomas (3 patients) and carcinomas (5 patients), and LC-MS/MS was used for protein identification and label-free quantification. The results were combined with a database search to generate a list of 106 proteins with nuclear annotations and at least a three-fold expression difference. Here, gene names from Table 2 of Knol et al. (2014) were converted to UniProt IDs.

UNS+14: Uzozie et al. (2014) analyzed 30 samples of colorectal adenomas and paired normal mucosa using iTRAQ labeling, OFFGEL electrophoresis and LC-MS/MS. 111 proteins with expression fold changes (log₂) at least ±0.5 and statistical significance threshold q < 0.02 that were also quantified in cell-line experiments were classified as “epithelial cell signature proteins”. UniProt IDs were taken from Table III of Uzozie et al. (2014).

WKP+14: de Wit et al. (2014) analyzed the secretome of paired CRC and normal tissue from 4 patients, adopting a five-fold enrichment cutoff for identification of candidate biomarkers. Here, the list of proteins from Supplementary Table 1 of de Wit et al. (2014) was filtered to include those with at least five-fold greater or lower abundance in CRC samples and p < 0.05. Two proteins listed as “Unmapped by Ingenuity” were removed, and gene names were converted to UniProt IDs using the UniProt mapping tool.

STK+15: Sethi et al. (2015) analyzed the membrane-enriched proteome from tumor and adjacent normal tissues from 8 patients using label-free nano-LC-MS/MS to identify 184 proteins with a fold change > 1.5 and p-value < 0.05. Here, protein identifiers from Supporting Table 2 of Sethi et al. (2015) were used to find the corresponding UniProt IDs.

WDO+15: Wiśniewski et al. (2015) analyzed 8 matched formalin-fixed and paraffin-embedded (FFPE) samples of normal tissue (N) and adenocarcinoma (C) and 16 nonmatched adenoma samples (A) using LC-MS to identify 2300 (N/A), 1780 (A/C) and 2161 (N/C) up- or down-regulated proteins at p < 0.05. The list of proteins used in this study includes only those marked as having a significant change in SI Table 3 of Wiśniewski et al. (2015).

LPL+16: Li et al. (2016) used iTRAQ and 2D LC-MS/MS to analyze pooled samples of stroma purified by laser capture microdissection (LCM) from 5 cases of non-neoplastic colonic mucosa (NC), 8 of adenomatous colon polyps (AD), 5 of colon carcinoma in situ (CIS) and 9 of invasive colonic carcinoma (ICC). A total of 222 differentially expressed proteins between NC and other stages were identified. Here, gene symbols from Supplementary Table S3 of Li et al. (2016) were converted to UniProt IDs using the UniProt mapping tool.

PHL+16: Peng et al. (2016) used iTRAQ 2D LC-MS/MS to analyze pooled samples from 5 cases of normal colonic mucosa (NC), 8 of adenoma (AD), 5 of carcinoma in situ (CIS) and 9 of invasive colorectal cancer (ICC). A total of 326 proteins with differential expression between two successive stages (and, for CIS and ICC, also differentially expressed with respect to NC) were detected. The list of proteins used in this study was generated by converting the gene names in Supplementary Table 4 of Peng et al. (2016) to UniProt IDs using the UniProt mapping tool.

Basis I

To formulate a thermodynamic description of a chemically reacting system, an important choice must be made regarding the basis species used to describe the system. The basis species, like thermodynamic components, are a minimum number of chemical formula units that can be linearly combined to generate the composition of any chemical species in the system of interest. Stated differently, any species can be formed by combining the components, but components can not be used to form other components (VanBriesen & Rittmann, 1999). Within these constraints, any specific choice of a basis is theoretically permissible. In making the choice of components, convenience (Gibbs, 1875), ease of interpretation and relationship with measurable variables, as well as availability of thermodynamic data (e.g., Helgeson, 1970), and kinetic favorability (May & Murray, 2001) are other useful considerations. Once the basis species are chosen, the stoichiometric coefficients in the formation reaction for any chemical species are algebraically determined.

Following previous studies (e.g., Dick, 2008), the basis species initially chosen here are CO₂, H₂O, NH₃, H₂S and O₂ (Basis I). The reaction representing the overall formation from these basis species of a protein having the formula C_cH_hN_nO_oS_s is (R1)${\displaystyle c{\mathrm{CO}}_2+n{\mathrm{NH}}_3+s{\mathrm{H}}_{2}S+n_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{H}}_{2}O+n_{{\mathrm{O}}_2}{\mathrm{O}}_2\to {\mathrm{C}}_c{\mathrm{H}}_h{\mathrm{N}}_n{\mathrm{O}}_o{\mathrm{S}}_s}$ where $n_{{\mathrm{H}}_2\mathrm{O}}=\left(h-3n-2s\right)∕2$ and $n_{{\mathrm{O}}_2}=\left(o-2c-n_{{\mathrm{H}}_2\mathrm{O}}\right)∕2$. Dividing n_H2O by the length of the protein gives the water demand per residue (${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$), which is used here because proteins in the comparisons generally have different sequence lengths.

These or similar sets of inorganic species (such as H₂ instead of O₂) are often used in studying reaction energetics in geobiochemistry (e.g., Shock & Canovas, 2010). However, as seen in Figs. 1A and 1B, there is a high correlation between Z_C of protein molecules and ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ in the reactions to form the proteins from Basis I (note that the choice of basis species affects only ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ and not Z_C). Because of this stoichiometric interdependence, changes in either redox or hydration potential, while holding the chemical potentials of the remaining basis species constant, have correlated effects on the energetics of chemical transformations (see ‘Comparison with inorganic basis species’ below). A different set of basis species can be chosen that reduces this correlation and affords a more informative description of the compositional changes in proteomic transformations.

Basis II

In this exploratory study, we restrict attention to at most two variables, with the implication that the others are held constant. In a subcellular setting, assuming that the chemical potentials of CO₂, NH₃ and H₂S do not change during a proteomic transformation, as implied by varying the chemical potentials of O₂ and H₂O in Basis I, may be less appropriate than assuming constant chemical potentials of more complex metabolites. In thermodynamic models for systems of proteins, constant chemical activities of chemical components having the compositions of amino acids might be a reasonable provision.

Although 1140 3-way combinations can be made of the 20 common proteinogenic amino acids, only 324 of the combinations contain cysteine and/or methionine (one of these is required to provide sulfur), and of these only 300, when combined with O₂ and H₂O, are compositionally nondegenerate. The slope, intercept and R² of the linear least-squares fits between Z_C and ${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$ using each possible basis containing O₂, H₂O and three amino acids are listed in file AAbasis.csv in Dataset S1. Many of these combinations have lower R² and lower slopes than found for Basis I (Figs. 1A and 1B), indicating a decreased correlation. From those with a lower correlation, but not the lowest, the basis including cysteine (Cys), glutamic acid (Glu), glutamine (Gln), O₂ and H₂O (Basis II) has been selected for use in this study. The scatterplots and fits between Z_C and ${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$ using Basis II are shown in Figs. 1C and 1D.

A secondary consideration in choosing this basis instead of others with even lower R² is the centrality of glutamine and glutamic acid in many metabolic pathways (e.g., DeBerardinis & Cheng, 2010). Accordingly, these amino acids may be kinetically more reactive than others in pathways of protein synthesis and degradation. The presence of side chains derived from cysteine and glutamic acid in the abundant glutathione molecule (GSH), associated with redox homeostasis, is also suggestive of a central metabolic requirement for these amino acids. Again, it must be stressed that the current provisional choice of basis species is neither uniquely determined nor necessarily optimal for a thermodynamic description of any particular system. More experience with thermodynamic modeling and better biochemical intuition will likely provide reasons to refine these calculations using a different basis, perhaps including metabolites other than amino acids.

A general formation reaction using Basis II is (R2)${\displaystyle n_{\mathrm{Cys}}{\mathrm{C}}_3{\mathrm{H}}_7{\mathrm{NO}}_2\mathrm{S}+n_{\mathrm{Glu}}{\mathrm{C}}_5{\mathrm{H}}_9{\mathrm{NO}}_4+n_{\mathrm{Gln}}{\mathrm{C}}_5{\mathrm{H}}_{10}{\mathrm{N}}_2{\mathrm{O}}_3+n_{{\mathrm{H}}_2\mathrm{O}}{\mathrm{H}}_2\mathrm{O}+n_{{\mathrm{O}}_2}{\mathrm{O}}_2\to {\mathrm{C}}_c{\mathrm{H}}_h{\mathrm{N}}_n{\mathrm{O}}_o{\mathrm{S}}_s}$

where the reaction coefficients (n_Cys, n_Glu, n_Gln, n_H2O and n_O2) can be obtained by solving (2)${\displaystyle \left[\begin{array}{ccccc}3\hfill & 5\hfill & 5\hfill & 0\hfill & 0\hfill \\ 7\hfill & 9\hfill & 10\hfill & 2\hfill & 0\hfill \\ 1\hfill & 1\hfill & 2\hfill & 0\hfill & 0\hfill \\ 2\hfill & 4\hfill & 3\hfill & 1\hfill & 2\hfill \\ 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\times \left[\begin{array}{c}{n}_{\mathrm{Cys}}\hfill \\ n_{\mathrm{Glu}}\hfill \\ n_{\mathrm{Gln}}\hfill \\ n_{{\mathrm{H}}_2\mathrm{O}}\hfill \\ n_{{\mathrm{O}}_2}\hfill \end{array}\right]=\left[\begin{array}{c}c\hfill \\ h\hfill \\ n\hfill \\ o\hfill \\ s\hfill \end{array}\right].}$ Although the definition of basis species requires that they are themselves compositionally nondegenerate, the matrix equation emphasizes the interdependence of the stoichiometric reaction coefficients. A consequence of this multiple dependence is that single variables such as n_O2 and n_H2O are not simple variables, but are influenced by both the chemical composition of the protein and the choice of basis species used to describe the system.

The combination of molecules shown in Reaction (R2) does not represent the actual mechanism of synthesis of the proteins. Instead, reactions such as this account for mass-conservation requirements and permit the subsequent generation of thermodynamic models for the potential for formation of different proteins as a function of system parameters (i.e., chemical potentials of O₂ and H₂O).

As an example of a specific calculation, consider the following reaction for MUC1, a chromatin-binding protein that is highly up-expressed in CRC cells (Knol et al., 2014). (R3)${\displaystyle 7{\mathrm{C}}_3{\mathrm{H}}_7{\mathrm{NO}}_2\mathrm{S}+535.6{\mathrm{C}}_5{\mathrm{H}}_9{\mathrm{NO}}_4+515.2{\mathrm{C}}_5{\mathrm{H}}_{10}{\mathrm{N}}_2{\mathrm{O}}_3{\to \mathrm{C}}_{5275}{\mathrm{H}}_{8231}{\mathrm{N}}_{1573}{\mathrm{O}}_{1762}{\mathrm{S}}_7+895.2{\mathrm{H}}_2\mathrm{O}+522.4{\mathrm{O}}_2.}$

As with the other reactions shown above, this reaction is not a mechanism, but represents the stoichiometric requirements for the formation from the basis species of one mole of the protein. Water is released in Reaction (R3), so the water demand (n_H2O) is negative. The length of this protein is 1,255 amino acid residues, giving the water demand per residue, ${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}=-895.2∕1,255=-0.71$. The average oxidation state of carbon (Z_C) in MUC1, which does not depend on the choice of basis species, is 0.005 (Eq. 1). The value of Z_C indicates that MUC1 is a relatively highly oxidized protein, while its ${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$ places it near the median water demand for up-expressed proteins in cancer in this dataset (see Fig. 2A below).

Thermodynamic calculations

Standard molal thermodynamic properties of the amino acids and unfolded proteins estimated using amino acid group additivity were calculated as described by Dick, LaRowe & Helgeson (2006), taking account of updated values for the methionine sidechain group (LaRowe & Dick, 2012). In this study, the Gibbs energies of hypothetically non-ionized proteins were used, and calculations were carried out at 37 °C and 1 bar. The temperature dependence of standard Gibbs energies was calculated using the revised Helgeson–Kirkham–Flowers (HKF) equations of state (Helgeson, Kirkham & Flowers, 1981; Tanger IV & Helgeson, 1988). Thermodynamic properties for O₂ (gas) were calculated using data from Wagman et al. (1982) and the Maier–Kelley heat capacity function (Kelley, 1960). Properties of liquid H₂O were calculated using data and extrapolations coded in Fortran subroutines from the SUPCRT92 package (Johnson, Oelkers & Helgeson, 1992), as provided in the CHNOSZ package (Dick, 2008).

Chemical affinities of reactions were calculated using activities of amino acids in the basis equal to 10⁻⁴, and activities of proteins equal to 1/(protein length) (i.e., unit activity of amino acid residues). Continuing with the example of Reaction (R3), an estimate of the standard Gibbs energy ($\mathrm{\Delta }{G}_f^o$) of the aqueous protein molecule (Dick, LaRowe & Helgeson, 2006; LaRowe & Dick, 2012) at 37 °C is −40, 974 kcal/mol; combined with the standard Gibbs energies of the basis species, this give a standard Gibbs energy of reaction ($\mathrm{\Delta }{G}_r^o$) equal to 66,889 kcal/mol. At loga_H2O = 0 and logf_O2 = − 65, with activities of the amino acid basis species equal to 10⁻⁴ , the overall Gibbs energy (ΔG_r) is 24,701 kcal/mol. The negative of this value is the chemical affinity (A) of the reaction. The per-residue chemical affinity for formation of protein MUC1 in the stated conditions is −19.7 kcal/mol. (This calculation can be reproduced using the function reaction() in file plot.R in Dataset S1.)

In a given system, proteins with higher (more positive) chemical affinity are relatively energetically stabilized, and theoretically have a higher propensity to be formed. Therefore, the differences in affinities reflect not only the amino acid compositions of the protein molecules but also the potential for local environmental conditions to influence the relative abundances of proteins.

Weighted rank difference

The contours on relative stability diagrams for the groups of differentially expressed proteins (see Fig. 6 below) depict the weighted rank differences of chemical affinities of formation of proteins. To illustrate this calculation, consider a hypothetical system composed of 3 proteins with higher expression in cancer (C) and 4 with higher expression in normal samples (down-expressed in cancer, i.e., having higher expression in a healthy state) (H). Suppose that under one set of conditions (i.e., specified loga_H2O and logf_O2), the per-residue affinities of the proteins give the following ranking in ascending order (I): ${\displaystyle \begin{array}{ccccccc}\mathrm{C}\hfill & \mathrm{C}\hfill & \mathrm{C}\hfill & \mathrm{H}\hfill & \mathrm{H}\hfill & \mathrm{H}\hfill & \mathrm{H}\hfill \\ 1\hfill & 2\hfill & 3\hfill & 4\hfill & 5\hfill & 6\hfill & 7\hfill \end{array}}$

This gives as the sum of ranks for up-expressed (C) proteins ∑r_C = 6, and for down-expressed (H) proteins ∑r_H = 22. The difference in sum of ranks is Δr_C−H = − 16; the negative value is associated with a higher rank sum for the down-expressed proteins, indicating that these as a group are more stable than the down-expressed proteins. In a second set of conditions, we might have (II): ${\displaystyle \begin{array}{ccccccc}\mathrm{H}\hfill & \mathrm{H}\hfill & \mathrm{H}\hfill & \mathrm{H}\hfill & \mathrm{C}\hfill & \mathrm{C}\hfill & \mathrm{C}\hfill \\ 1\hfill & 2\hfill & 3\hfill & 4\hfill & 5\hfill & 6\hfill & 7\hfill \end{array}}$ Here, the difference of rank sums is Δr_C−H = 18 − 10 = 8.

For systems where the numbers of proteins in the two groups are equal, the maximum possible differences in rank sums would have equal absolute values, but that is not the case in this and other systems having unequal numbers of up- and down-expressed proteins. To characterize these datasets, the weighted rank-sum difference can be calculated using (3)${\displaystyle \mathrm{\Delta }\overline{r}=2\left(\frac{n_{\mathrm{H}}}{n}\sum r_{\mathrm{C}}-\frac{n_{\mathrm{C}}}{n}\sum r_{\mathrm{H}}\right)}$ where n_H, n_C and n are the numbers of down-expressed, up-expressed, and total proteins in the comparison. In the example here, we have n_H∕n = 4∕7 and n_C∕n = 3∕7. Equation (3) then gives $\mathrm{\Delta }\overline{r}=-12$ and $\mathrm{\Delta }\overline{r}=12$, respectively, for conditions (I) and (II) above, showing equal weighted rank-sum differences for the two extreme rankings.

We can also consider a situation where the ranks of the proteins are evenly distributed: ${\displaystyle \begin{array}{ccccccc}\mathrm{H}\hfill & \mathrm{C}\hfill & \mathrm{H}\hfill & \mathrm{C}\hfill & \mathrm{H}\hfill & \mathrm{C}\hfill & \mathrm{H}\hfill \\ 1\hfill & 2\hfill & 3\hfill & 4\hfill & 5\hfill & 6\hfill & 7\hfill \end{array}}$ Here the absolute difference of rank sums is Δr_C−H = 12 − 16 = − 4, but the weighted rank-sum difference is $\mathrm{\Delta }\overline{r}=0$. The zero value for an even distribution and the opposite values for the two extremes demonstrate the applicability of this weighting scheme.

Software and data availability

All statistical and thermodynamic calculations were performed using R (R Core Team, 2016). Thermodynamic calculations were carried out using R package CHNOSZ (Dick, 2008). Effect sizes (see below) were calculated using R package orddom (Rogmann, 2013). Figures were generated using CHNOSZ and graphical functions available in R together with the R package colorspace (Ihaka et al., 2015) for constructing an HCL-based color palette (Zeileis, Hornik & Murrell, 2009). With the mentioned packages installed, Table 1 and the figures in this paper can be reproduced using the code (plot.R) and data files (*.csv) in Dataset S1.

Compositional comparisons of human proteins

Comparisons of proteome composition in terms of average oxidation state of carbon (Z_C) and water demand per residue (${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$) are presented in Fig. 2 and Table 1. Figure 2 shows scatterplots of individual protein compositions for proteomes in three representative studies. Each of these exhibits a strongly differential trend in Z_C or ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ that can be visually identified. In Fig. 2A, chromatin-binding proteins highly expressed in carcinoma (Knol et al., 2014) as a group exhibit a lower Z_C than those found to be more abundant in adenoma. In Fig. 2B, proteins relatively highly expressed in epithelial cells in adenoma (Uzozie et al., 2014) tend to have higher ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ than the proteins more highly expressed in paired normal tissues. Differentially expressed proteins between adenoma and normal tissue identified in a recent deep-proteome analysis (Wiśniewski et al., 2015) are compared in Fig. 2C, showing that proteins up-expressed in adenoma are relatively oxidized (i.e., have higher Z_C).

In order to quantify these differences, Table 1 shows the numbers of proteins in each comparison (n₁ for normal or less advanced cancer stage; n₂ for tumor or more advanced cancer stage), differences of means (MD), common language effect size as percentages (ES), and p-values calculated using the Wilcoxon rank sum test. This non-parametric test is suitable for data which may not be normally distributed. For a given experiment, the common language effect size, or probability of superiority, describes the probability that Z_C or ${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$ of a protein is higher in the cancer group than in the normal group. That is, percent values of the ES greater than (or less than) 50 indicate a greater proportion of pairwise higher (or lower) Z_C or ${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$ of proteins in the n₂ compared to n₁ groups. The ES and p-value are used here to allow for a subjective assessment of the compositional differences. ES values ≥60 or ≤40 and p-values < 0.05 are highlighted in the table. The corresponding mean differences are underlined for p < 0.05, or bolded if ES is also ≥60 or ≤40. These cutoffs highlight datasets with the largest and most significant differences in Z_C and ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$. Mean and median values of Z_C and ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ are given in file summary.csv in Dataset S1.

Counting the underlined and bolded MD values in Table 1, the number of datasets with a significant difference in Z_C (18) is greater than those with a significant difference in ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ (10). Of the 10 unique studies yielding at least one dataset with a significant difference in Z_C in a comparison with normal tissue, 8 exhibit a higher mean value in adenoma or carcinoma compared to normal tissue. One of the other studies (Besson et al., 2011) has datasets with higher mean Z_C in proteins up-expressed in adenoma, but lower mean Z_C in proteins up-expressed in stage II and III carcinoma, compared to normal tissue. A second study, which analyzed proteins in stromal cells (Li et al., 2016), shows a significantly lower Z_C in adenoma compared to normal tissue.

Most of the studies analyzed proteins in whole or microdissected tissue, but two datasets in studies from the same laboratory represent the nuclear matrix or chromatin-binding fraction (Albrethsen et al., 2010; Knol et al., 2014). These two datasets give lower mean Z_C of proteins more highly expressed in carcinoma than adenoma. One other dataset has a lower mean Z_C of proteins up-expressed in carcinoma compared to adenoma (Wiśniewski et al., 2015), and one has a higher mean value (Mikula et al., 2011).

The datasets with a significant difference in ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ all show higher mean values for proteins in adenoma (5 datasets) or carcinoma (3 datasets) compared to normal tissue, up- expressed compared to down-expressed serum biomarker candidates (Jimenez et al., 2010), and secreted proteins of tumor tissue compared to normal tissue (de Wit et al., 2014). Interestingly, none of the datasets with a significant difference in ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ corresponds to a carcinoma/adenoma comparison.

Natural variability inherent in the heterogeneity of tumors, as well as differences in experimental design and technical analysis, may underlie the opposite trends in Z_C between some datasets that compare the same stages of cancer. However, there is a preponderance of datasets with higher values of Z_C and ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ for the proteins up-expressed in adenoma or carcinoma compared to normal tissue.

Compositional comparisons of microbial proteins

Summary data on microbial populations from four studies were selected for comparison here. First, in a study of 16S RNA of fecal microbiota, Wang et al. (2012) reported genera that are significantly increased or decreased in CRC compared to healthy patients. In order to compare the chemical compositions of the microbial populations, single species with sequenced genomes were chosen to represent each of these genera (see Table 2). Where possible, the species selected are those mentioned by Wang et al. (2012) as being significantly altered, or are species reported in other studies to be present in healthy or cancer states (see Table 2).

In the second study considered (Zeller et al., 2014), changes in the metagenomic abundance of fecal microbiota associated with CRC were analyzed for their potential as a biosignature for cancer detection. The species shown in Fig. 1A of Zeller et al. (2014) with a log odds ratio greater than 0.15 were selected for comparison, and are listed in Table 3. Zeller et al. (2014) found a strong enrichment of Fusobacterium in cancer, consistent with previous reports (Kostic et al., 2012; Castellarin et al., 2012). In a third study, Candela et al. (2014) reported the findings of a network analysis that identified 5 microbial “co-abundance groups” at the genus level. As before, single representative species were selected in this study, and are listed in Table 2. Except for the presence of Fusobacterium, the co-abundance groups show little genus-level overlap with community profiles derived from the previous two studies.

Finally, Table 4 lists the “best aligned strain” from Supplementary Dataset 5 of Feng et al. (2015) for all species shown there with negative enrichment in cancer, and for selected species with positive enrichment in cancer. Although every uniquely named strain given by Feng et al. (2015) was used in the comparisons below (n = 44; see Fig. 3D below), for clarity only the up-enriched species that appear in the calculated stability diagram (see Fig. 4D below) are listed in Table 4 and labeled in Fig. 3D. File microbes.csv in Dataset S1 contains the complete list of Bioproject IDs and calculated Z_C and ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$ for all the microbial species considered here.

For each of the microbial species listed in Tables 2–4, a mean protein composition was calculated by combining amino acid sequences of all proteins downloaded from the NCBI genome page associated with the Bioproject IDs shown in the Tables (see file microbial.aa.csv in Dataset S1). This method does not account for actual protein abundances in organisms, and excludes any post-translational modifications. Calculation of the mean amino acid composition of proteins in this way is not an exact representation of the cellular protein composition, but provides a starting point for identifying environmental signals in protein composition. Mean amino acid compositions or amino acid frequencies deduced from microbial genomes, calculated without weighting for actual protein abundance, have been used in many studies making evolutionary and/or environmental comparisons (e.g., Tekaia & Yeramian, 2006; Zeldovich, Berezovsky & Shakhnovich, 2007; Brbić et al., 2015). In the future, more refined calculations may be possible by using genome-wide estimates of protein expression levels based on codon usage patterns (e.g., Moura, Savageau & Alves, 2013; Brbić et al., 2015).

The water demand per residue (${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$) vs. oxidation state of carbon (Z_C) in the mean amino acid compositions of proteins from all of the microbial species considered here are plotted in Figs. 1B and 1D, and for individual datasets in Fig. 3. The groups of proteins in the microbes enriched in cancer patients have somewhat lower Z_C than those enriched in healthy donors in the same study. The dataset from Feng et al. (2015) (Fig. 3D) shows a more complex distribution, where the microbes with a relative enrichment in healthy individuals form two clusters at high and low Z_C. The Fusobacterium species identified in the studies of Zeller et al. (2014), Candela et al. (2014) and Feng et al. (2015) have the lowest Z_C of any microbial species considered here. The mean human protein composition is also plotted in Fig. 3, revealing a higher Z_C than any of the mean microbial proteins except for Actinomyces viscosus and Bifidobacterium animalis, identified in the study of Feng et al. (2015) (Fig. 3D). The tendency for microbial organisms to be composed of more reduced biomolecules than the host may reflect the relatively reducing conditions in the gut.

Thermodynamic descriptions: background

Going beyond compositional comparisons, thermodynamic descriptions can account for stoichiometric and energetic constraints and provide a richer interpretation of proteomic data in the context of tumor microenvironments.

By combining both stoichiometric and energetic variables, a thermodynamic description of proteomic data reveals possible biochemical constraints that may arise within cells and in tumor microenvironments. To give an example of how relative stabilities of up- and down-expressed proteins in a proteomic dataset can be calculated as a function of chemical potentials, consider Reaction (R3) above written for the formation of one mole of MUC1. In order to compare proteins of different lengths, the formula of the protein is written per residue. The corresponding reaction is then (R4)${\displaystyle 0.006{\mathrm{C}}_3{\mathrm{H}}_7{\mathrm{NO}}_2\mathrm{S}+0.427{\mathrm{C}}_5{\mathrm{H}}_9{\mathrm{NO}}_4+0.411{\mathrm{C}}_5{\mathrm{H}}_{10}{\mathrm{N}}_2{\mathrm{O}}_3{\to \mathrm{C}}_{4.203}{\mathrm{H}}_{8.557}{\mathrm{N}}_{1.253}{\mathrm{O}}_{2.403}{\mathrm{S}}_{0.006}+0.714{\mathrm{H}}_2\mathrm{O}+0.416{\mathrm{O}}_2.}$

An expression for the chemical affinity (Kondepudi & Prigogine, 1998; Helgeson et al., 2009) of Reaction (R4) is (4)${\displaystyle \mathit{A}=2.303RTlog\left(K∕Q\right)}$ where 2.303 is shorthand for the natural logarithm of 10, R is the gas constant, T is temperature, log represents the common (decimal) logarithm, and the activity quotient Q is given by (5)${\displaystyle logQ=log{a}_{{\mathrm{C}}_{4.203}{\mathrm{H}}_{8.557}{\mathrm{N}}_{1.253}{\mathrm{O}}_{2.403}{\mathrm{S}}_{0.006}}+0.714log{a}_{{\mathrm{H}}_2\mathrm{O}}+0.416log{f}_{{\mathrm{O}}_2}-0.006log{a}_{{\mathrm{C}}_3{\mathrm{H}}_7{\mathrm{NO}}_2\mathrm{S}}-0.427log{a}_{{\mathrm{C}}_5{\mathrm{H}}_9{\mathrm{NO}}_4}-0.411log{a}_{{\mathrm{C}}_5{\mathrm{H}}_{10}{\mathrm{N}}_2{\mathrm{O}}_3}.}$

The equilibrium constant is given by $-2.303RTlogK=\mathrm{\Delta }{G}_r^o$, where $\mathrm{\Delta }{G}_r^o$ is the standard Gibbs energy of the reaction. As noted above, the standard Gibbs energies of species used to calculate $\mathrm{\Delta }{G}_r^o$ at T = 37 °C are generated using amino acid group additivity for the proteins and published values for standard thermodynamic properties of the basis species in the reaction.

To compare the potential for formation of metastable molecules, the per-residue formulas of the proteins are assigned equal activities (1). Then, Eqs. (4)–(5) show that the affinity, and hence relative potential for formation of different proteins, is a function of not only their amino acid composition (which determines the chemical formulas and standard Gibbs energies of the proteins used here), but also system parameters including temperature and the chemical potentials of the components. In this study, the chemical activities of the amino acid basis species are provisionally set to constant values (10⁻⁴), while logf_O2 and loga_H2O are considered to be adjustable parameters that are used as exploratory variables. The ranges of these variables shown on the diagrams are selected in order to encompass the stability boundaries between groups of proteins differentially enriched in cancer and normal samples.

There are combinations of chemical activities of basis species in Eq. (5) where the per-residue formation reactions of two proteins have an equal affinity, indicating equal chemical stability of the proteins. Other combinations of chemical activities of basis species give the result that one protein-residue formula has a higher affinity than the other(s), indicating greater stability of this protein. This concept provides an approach for constructing stability diagrams, which may be called the “maximum affinity method”, that can be used to reproduce published equilibrium and metastable equilibrium diagrams for many inorganic and organic systems as shown by examples in the CHNOSZ package (Dick, 2008) and is used below for microbial proteins. Because of the greater numbers of individual proteins in human proteomic datasets, a new method based on the difference in weighted sums of ranks of affinities is used here to compare the relative stabilities of groups of up- and down-expressed proteins in cancer.

Relative stability fields for microbial proteins

Stability diagrams are shown in Figs. 4A–4D for the four sets of microbial proteins described above. The first diagram, representing significantly changed genera detected in fecal 16S rRNA (Wang et al., 2012; first part of Table 2), shows maximal stability fields for proteins from 5 species relatively enriched in healthy patients, and 3 species enriched in CRC patients. The other 4 proteins in the system are less stable than the others within the range of logf_O2 and loga_H2O shown and do not appear on the diagram. The relative positions of the stability fields in Fig. 4A are roughly aligned with the values of Z_C and ${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$ of the proteins; note for example the high-logf_O2 positions of the fields for the relatively high-Z_C Escherichia coli and Alistipes indistinctus, and the high-loga_H2O position of the field for the high-${\overline{n}}_{{\mathrm{H}}_2\mathrm{O}}$ Peptostreptococcus stomatis. Except for E. coli, the proteins from the species enriched in CRC occupy the lower logf_O2 (reducing) and higher loga_H2O zones of this diagram.

In thermodynamic calculations for proteins from bacteria detected in fecal metagenomes (Zeller et al., 2014; Table 3), the mean protein compositions of 3 of 6 healthy-enriched microbes and 4 of 7 cancer-enriched microbes exhibit maximal relative stability fields (Fig. 4B). Here, the cancer-associated proteins occupy the more reducing (Fusobacterium nucleatum subsp. vincentii and subsp. animalis) or more oxidizing (Clostridium hylemonae, Porphyromonas asaccharolytica) regions, while the proteins from bacteria more abundant in healthy individuals are relatively stable at moderate oxidation–reduction conditions.

For the bacterial species representing microbial co-abundance groups (Candela et al., 2014; second part of Table 2), all of the 5 mean protein compositions are present on the diagram (Figs. 4C). Here, the proteins from cancer-enriched bacteria are more stable at reducing conditions and those from healthy-enriched microbes are stabilized by oxidizing conditions.

A stability diagram for proteins of bacteria identified in a second metagenomic study (Feng et al., 2015) shows a similar result (Fig. 3D) for the 11 mean protein compositions with highest stability at some point the diagram. These patterns in relative stability again reflect the differences in Z_C of the proteins, although in this case, a greater proportion of proteins (33 out of the 44 included in the calculations) are not found to have maximal stability fields. The resulting stability diagram is therefore a more limited portrayal of the available data.

Figure 4E is a composite representation of the calculations, in which higher cumulative counts of maximal stability of proteins from bacteria enriched in normal and cancer samples in the four studies are represented by deeper blue and red shading, respectively. According to this diagram, the chemical conditions predicted to be most favorable for the formation of proteins in many bacteria enriched in CRC are characterized by low logf_O2. Proteins from bacteria that are abundant in healthy patients tend to be stabilized by moderate values of logf_O2. Despite the differences in experimental design and microbial identification between studies, the thermodynamic calculations reveal a shared pattern of relative stabilities among the four datasets considered here.

Relative stability fields for human proteins

Diagrams like those shown above that portray the maximally stable protein compositions are inadequate for analysis of larger datasets such as those generated in proteomic studies. It is apparent in Fig. 5 that only three different proteins up-expressed in cancer, from the 106 proteins in the KWA+14 dataset (chromatin-binding proteins in carcinoma/adenoma), are maximally stable across a range of logf_O2. However, visual inspection reveals a differential sensitivity to oxygen fugacity in the whole dataset, with lower logf_O2 providing relatively higher potential for the formation of many of the up-expressed proteins in carcinoma samples. How can these responses be quantified in order to explore the data in multiple dimensions, including both loga_H2O and logf_O2?

In Fig. 5B, the difference in mean values of chemical affinity per residue of carcinoma and adenoma-associated proteins appears as a straight line as a function of logf_O2. This linear behavior would translate to evenly spaced isostability (taken as constant mean affinity difference) contours on a logf_O2–loga_H2O diagram. The weighted rank difference of affinities (see Methods), shown by the curved line Fig. 5B, is a summary function that is more informative of changing chemical conditions. The variable slope is greatest near the zone of convergence for affinities of individual proteins (Fig. 5A), corresponding to the transition zone between groups of proteins. The resulting two-dimensional stability diagrams shown below have curved and diversely spaced isostability (taken as constant weighted rank difference of affinity) contours.

The diagrams in Fig. 6 portray weighted rank differences of chemical affinities of formation between groups of up- and down-expressed proteins reported for proteomic experiments. These combined depictions of stoichiometric and energetic differences constitute a theoretical prediction of the relative chemical (not conformational) stabilities of the proteins.

The slopes of the equal-stability lines and the positions of the stability fields reflect the magnitude and sign of differences in Z_C and ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$. Figures 6A–6C show results for datasets that are dominated by differences in ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$; the nearly horizontal lines show that relative stabilities are accordingly more sensitive to loga_H2O than logf_O2. The second row depicts relative stabilities in the three datasets from Mikula et al. (2011), which have large changes in, sequentially, ${\bar{n}}_{{\mathrm{H}}_2\mathrm{O}}$, Z_C, then both of these (Table 1). Accordingly, the equal-stability lines for these datasets are closer to horizontal, closer to vertical, or have a more diagonal trend (Figs. 6D–6F).

The last row shows results for datasets that are characterized by large changes in Z_C; the relative stabilities depend strongly on logf_O2. According to Fig. 6G, higher oxygen fugacity increases the relative potential for the formation of proteins up-expressed in cancer (dataset of Jankova et al., 2011). However, in a dataset for up- and down-expressed chromatin-binding proteins in carcinoma (Knol et al., 2014), lower logf_O2 is predicted to promote formation of the proteins up-expressed in carcinoma. This is the opposite trend to that found for most of the other datasets with significant differences in Z_C. These opposing trends might be attributed to different biochemical constraints acting at the subcellular and cellular or tissue levels during carcinogenesis.

The full set of diagrams for all datasets listed in Table 1 is provided in Fig. S1. It is notable that for the datasets where the relative stabilities are strongly a function of loga_H2O (sub-horizontal lines), the equal-stability lines are within a few log units of 0 (unit activity). Equal-stability lines that are diagonal often cross unit activity of H₂O at a moderate value of logf_O2, near −65 to −60 (see Fig. S1). This could be indicative of a tendency for these proteomic transformations to be partially buffered by other redox reactions in the cell, and/or by liquid-like H₂O with close to unit activity.

Effective values of oxidation–reduction potential (Eh) can be calculated by considering the water dissociation reaction, i.e., (R5)${\displaystyle {\mathrm{H}}_2\mathrm{O}\rightleftharpoons \frac{1}{2}{\mathrm{O}}_2+2{\mathrm{H}}^{+}+2e^{-}.}$ If one assumes that loga_H2O = 0 (unit water activity, as in an infinitely dilute solution), this reaction can be used to interconvert logf_O2, pH and pe (or, in conjunction with the Nernst equation, Eh) (e.g., Garrels & Christ, 1965, p. 176; Anderson, 2005, p. 363). However, in the approach utilized here for assessing the relative stabilities of proteins in a subcellular context, no such assumptions are made on the operational value of loga_H2O. Instead, it is used as an indicator of the internal state of the system, and is not necessarily buffered by an aqueous solution. Consequently, the effective Eh is considered to be a function of variable logf_O2 and loga_H2O, as shown in Fig. 6I for pH = 7.4 and T = 37 °C. This comparison gives some perspective on operationally reasonable ranges of logf_O2 and loga_H2O.

The subcellular reduction potential monitored by the reduced glutathione (GSH)/oxidized glutathione disulfide (GSSG) couple ranges from ca. −260 mV for proliferating cells to ca. −170 mV for apoptotic cells (Schafer & Buettner, 2001), lying toward the middle part of the range of conditions shown in Fig. 6. A physiologically plausible Eh value of −0.2 V, corresponding to logf_O2 = − 62.8 at unit activity of H₂O, is close to the stability transitions for many of the datasets considered here (see also Fig. S1).

Comparison with inorganic basis species

Figures made using Basis I (inorganic basis species, e.g., Reaction (R1)) are provided in the Supplemental Information (human proteins: Fig. S2; microbial proteins: Fig. S3). The stability boundaries in loga_H2O–logf_O2 diagrams constructed using Basis I cluster around a common, positive slope, in contrast with the greater diversity of slopes appearing on the corresponding diagrams constructed using Basis II (Fig. S1).

As noted above, all mathematically possible choices for the basis species of a system are thermodynamically valid, but it appears that Basis II affords a greater convenience for interpretation. That is, compared to Basis I, Basis II yields a greater degree of separation of the effects of changing chemical potentials of H₂O and O₂ under the assumption that the activities of the remaining basis species (inorganic species in Basis I, or amino acids in Basis II) are held constant. However, it is also notable that two of the diagrams constructed using Basis I (Fig. S2), unlike the others, have nearly horizontal equal-stability lines, showing that increasing activity of H₂O at constant activity of CO₂, NH₃, H₂S and fugacity of O₂ gives an energetic advantage to the formation of potential up-expressed serum biomarkers (dataset JKMF10; Jimenez et al., 2010) and proteins up-expressed in an “epithelial cell signature” for adenoma (dataset UNS+14; Uzozie et al., 2014). These datasets are also found to be among those having significantly differential water demand using Basis II (Table 1; Fig. S1). Based on the similar results for these datasets using different choices of chemical components, it can be suggested that the compositions of the differentially expressed proteins in these datasets are especially indicative of changes in hydration potential.