Model preprocessing.

We first ensured the model was consistent by identifying blocked reactions—whose maximum and minimum fluxes are null under open exchange conditions—with flux variability analysis, using the find_blocked_reactions function from the COBRApy package. We also removed any boundary metabolites—usually added to balance exchange reactions—prior to running any reconstruction, gapfill, or analysis method. This model must be feasible in steady-state conditions and must be capable of allowing flux through the biomass pseudo-reaction.

Transcriptomics data as a proxy of enzyme activity.

In this pipeline, we focused specifically on transcriptomics data that is mappable with the model’s gene associations, although some methods allow integration of other data types.

Similarly to previous approaches [19, 36, 44], we used transcriptomics data as a proxy for enzyme presence and flux activity from which we can calculate reaction activity scores (RAS) to serve as inputs for CSMR algorithms. These scores should ideally reflect whether a given reaction in the model is likely to be present in the context represented by the transcriptomics data. The work of Richelle et al. details the implications of several ways to infer RASs and provides several thresholding options [44], which we adapted as a part of our work. Out of the parameterization choices highlighted by the authors, we focused on varying the thresholding approach and GPR integration functions.

Using transcriptomics to characterize enzyme activity is not trivial, since the relationship between messenger RNA and protein expression is not fully understood, despite ongoing progress in quantifying both of these biological entities. However, Nusinow et al. have recently quantified the proteome for a subset of the CCLE panel and found a moderately positive correlation (mean Pearson c.c. = 0.48) between mRNA and protein abundance [45]. In this work, we assumed a linear relationship so that RASs were calculated based on gene expression measurements.

Scoring transcript activity from expression measurements.

RNA-Seq technologies typically produce transcript-level measurements represented as proportions of the entire transcriptome. However, we intended, for the RASs used in this work, to obtain a positive or negative value relative to reference thresholds calculated across all samples. To this end, we processed expression measurements into transcript activity scores (TASs) that can better represent this dichotomy. All the thresholding calculations were done over the transcript-level measurements.

We first defined the concept of a global threshold, where the expression of all genes contributes. This is useful in filtering out transcripts whose expression is low or high enough for them to be assigned as inactive or active, respectively, with a high degree of confidence. However, this type of thresholding does not take into account the variability of measurements between transcripts. While originally available for microarray-based transcriptome quantification [46], a generic expression threshold to barcode all cell types is hard to define and apply in RNA-Seq measurements, due to the different conditions in which experiments are performed.

A local thresholding approach can also be considered to mitigate the aforementioned problems. Rather than condensing the entire measurements into a single value, thresholding can be performed on a per gene basis, yielding a value for each gene independently. Similarly to the global thresholding approach, local thresholds can also be used as a reference for fold change calculations.

In this pipeline, both thresholds were calculated by first determining transcript-wise quantiles for various percentages (from 10% to 90%), yielding multiple sets of local thresholds, one for each percentage. To convert the latter to global thresholds, we used the mean value of a local threshold set to obtain a single representative value for the entire expression dataset.

After establishing appropriate thresholds, we then combined them to establish rules that can be used to determine whether a transcript is active and calculate its TAS to better represent that activity level. We implemented an approach based on the work of Richelle et al [44], where transcript activity can be represented in two main states, namely:

• Inactive: transcript expression is inactive with a high degree of confidence—assigned when expression values are lower than a global lower threshold (g_min). The expected TAS values will always be negative.
• Active: transcript expression is active with a high degree of confidence since its value exceeds a global upper threshold (g_max). The expected TAS values are always positive.

This also implies the existence of an intermediate state of uncertainty for cases where transcript expression lies between the two thresholds. In these cases, we distinguished between active and inactive transcripts by comparing the transcript’s expression with its transcript-specific local threshold (l(y)) and the expected TAS value is constrained between -1 and 1, with its sign reflecting whether it is considered active.

We implemented three thresholding strategies based on the work of Richelle et al [44] with some minor changes to accommodate for the expected distribution of TAS values across multiple states. A global thresholding strategy implies a single global threshold to distinguish between active and inactive transcripts. An extension of this strategy, named localT1, includes local thresholding for transcripts that would otherwise be considered as inactive, and assigns TASs based on the ratio between expression and the transcript’s local threshold. Finally, we also included a localT2 strategy defined by the usage of two thresholds and an intermediate state, as defined above.

TASs were then generated by calculating the ratio between measured expression values and an appropriate threshold which is chosen according to the state in which the transcript is assigned. A detailed description of the formulae used in each state for the three employed strategies can be found on Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Functions used to convert transcript expression values into transcript activity scores, assuming x as a vector of expression levels for each transcript.

The “Expression value” column represents the condition that values in x must meet for the corresponding reference threshold t used to calculate a ratio with the formula . Finally, the range of TAS values for each condition is detailed on the last column.

https://doi.org/10.1371/journal.pcbi.1009294.t001

Inferring reaction activity from transcript scores.

The TASs from the aforementioned strategies are then converted to RASs using the gene-protein-reaction (GPR) rules provided with the model. GPR rules are Boolean expressions that describe, for a given reaction, which combinations of transcripts are involved with the synthesis of one or more enzymes capable of catalyzing it. These rules are often expressed or can be converted into disjunctive normal form, where multiple conjunctions (expressions with the AND operator) denoting the various enzymes or isoforms involved are bound by a disjunction (OR operator).

RASs must be presented as continuous scores and, thus, the Boolean operators in GPR rules must be replaced with numerical values. We replaced AND operators with a minimum function—an enzyme’s activity is limited by the lowest expressed transcript/subunit—while OR operations could be replaced with either sum or maximum functions. When using sum, we assumed the reaction activity correlates with the combined activity of all enzymes and isoforms catalyzing it, while the maximum function equates reaction activity with the highest expressed enzyme’s score.

Normalizing inputs for context-specific reconstruction algorithms.

This step includes conversion of RASs into inputs accepted by the different CSMR algorithms, given their diverse nature, and we have implemented two alternatives to perform this conversion in our routines. Although our pipeline is generic, we chose the FASTCORE [34] and tINIT [47] algorithms for context-specific model reconstruction.

Methods such as tINIT, where scores mirror the reactions’ states as present or absent, can take RAS as input without any further processing. On the other hand, algorithms such as FASTCORE require a set of core reactions as input. In this case, a further threshold must be applied for these to be obtained. In our work, we emphasized the division between positive and negative scores to represent activity, and as such, core reactions are those with a RAS above 0.

Algorithm output and post-processing.

With the inputs appropriately adapted, the output of each algorithm is always a binary vector r of size n (equal to the number of reactions in the template model), indicating reactions’ presence. Indeed, this vector includes Boolean flags indicating whether each reaction should be kept or removed in the context-specific model.

The models generated by the CSMR algorithms were then checked for consistency with expected phenotypes. For each of these models, we knocked out (set lower and upper bounds to 0) reactions flagged for removal, before performing any simulation or analysis. We first checked whether the model is capable of allowing non-zero flux through the biomass reaction, to ensure lethality can be tested. When growth medium formulations are available, we can additionally ensure that the model is feasible and capable of growth if the compounds present in growth media are the only ones allowed to be consumed. This was achieved by constraining exchange reactions that do not involve medium metabolites to only allow positive flux values, thus only allowing medium metabolites to be consumed by the model.

Elementary flux mode-based gap filling approach.

Gap filling was performed using a novel EFMGapfill approach, which was implemented in the troppo Python package. This algorithm leverages efficient elementary flux mode (EFM) enumeration algorithms to find minimal sets of active fluxes required for feasibility under a specific condition. Assuming a stoichiometric matrix S of m metabolites and n fluxes, the flux vector v and an identically sized vector y, and a set K of reactions available to fill gaps, the LP formulation employed in EFMGapfill can be defined as follows: (1) (2) (LP1) (3) (4) (5) (6)

In the formulation represented in LP1, constraint 1 defines the steady-state constraint, similarly to other CB approaches, such as FBA. Constraints 2 and 3 associate the binary variables in y to the fluxes in the vector v. In this expanded solution space, the variables in y will hold a value of 1 if their associated flux in the vector v is greater than 1. Otherwise, both variables are set to 0. These variables and indicator constraints are then used to discretize fluxes into active and inactive states. The objective function is dependent on the set of reactions K available for the algorithm to add as a gap filling solution, although the objective is always to minimize the sum of a subset of the vector y whose indices are contained in K.

We adapted the input K according to a Boolean vector r (a set of reaction indices). K will have all reactions from the template model not included in r. The vector r is typically the output of a previously determined CSMR reconstruction. Furthermore, we also defined and constrained an objective reaction u (usually the biomass pseudo-reaction) to always carry non-zero flux, representing a phenotype that is expected to be maintained upon tailoring the model to the subset of reactions in r.

We identified two possible gap filling scenarios that can be accomplished using this approach. In the first, we did not assume constraints on external metabolite exchanges and thus, we also excluded these reactions from K. The resulting solution from our gap filling approach is the smallest set of intracellular reactions not found in r that should be included so that the context-specific model is capable of carrying flux through u.

An alternative scenario may arise where the set of reactions r must not be manipulated, but the model still requires gap filling to predict growth. The growth medium, rather than the enzyme content of this model must be the target for manipulation. In this case, all intracellular reactions not in r must be constrained and exchange reactions must be split into forward and reverse reactions carrying flux in opposite directions. To find the minimal set of extracellular metabolites required for the model to carry flux through u, the set K must be defined as the set of reverse exchange reactions in the model.

Gene essentiality.

Gene essentiality screens, such as those performed with CRISPR, provide a directly quantifiable measurement of the impact of gene deletions on cell viability, which can be modeled on CBMs through metabolic tasks [19]. The biomass objective function, included in most human models, groups most of these tasks’ demands by aggregating the necessary components for cell division and maintenance. Given the computational demand of checking multiple gene knockouts for each task and each omics sample, we will focus on predicting lethal gene knockouts using the biomass objective function as a measure of cell growth.

The CBM workflow used to predict essential genes used GPRs to determine the set of reactions to exclude given a knocked-out gene g. To this end, we first obtained a mapping ω(g, r), which evaluates the GPR expression of reaction r with every gene marked as active, except for g. To apply the gene knockout, we must first determine the set K = {r|r ∈ R, ∀¬ω(g, r)}, which identifies the reactions that are disabled upon deletion of g; then, we set the lower and upper bounds of each reaction in K to 0. Adding these constraints to the model, the simulation can be run using FBA, yielding predicted growth rates for each gene deletion.

Finally, flux distributions resulting from gene knockouts can be evaluated. It is useful to always compare predicted mutant growth rates with wild-type levels. We considered several growth rate thresholds based on the wild-type value to represent viability, but we assumed that values of the mutant biomass flux below 0.1% of the wild-type biomass flux were considered a knockout lethal. Additionally, infeasible solutions are considered as non-viable. We then discretized each gene knockout’s result as essential or non-essential and compared them with the experimental screening.

We used Matthews’ correlation coefficient (MCC) to assess the predictive ability of our models. The multiclass definition of MCC as implemented in the scikit-learn package is presented on Eq 7, assuming a generic classifier to predict K classes, t as a vector with the amount of true positives and p the vector with the amount of predictions each class k, while c is the total amount of true positive samples for all classes and s is the number of samples. (7)

Predicted fluxes.

Alternatively, flux distributions obtained from the model using an appropriate phenotype prediction method can be directly compared with experimentally measured fluxes, obtained from techniques such as isotope labeling coupled with metabolic flux analysis. In this work, we employed parsimonious enzyme usage flux balance analysis (pFBA) to predict phenotypes using our context-specific reconstructions. We have chosen this method since it requires no prior knowledge and reduces the admissible solution space of FBA by assuming cells not only attempt to achieve the predefined cell objective, but also minimise the overall sum of metabolic fluxes to do so.

A key limitation in using CB models to predict flux values is the lack of reliable measurements for substrate uptake fluxes. This directly influences the predicted growth rate and intracellular fluxes. A more reliable comparison can be made by discretizing flux values into three classes: forward active, if the flux is positive, reverse active if it is negative (flux is active, but carried in the reverse direction), or null when there is no flux. Although less precise, this discards the usage of experimentally measured external metabolite consumption rates. The model’s predictive ability can then be ascertained by using metrics suitable for multiclass predictive models such as Matthews’ correlation coefficient or weighted F1 scores.

Parameter importance.

We assessed the effect of the model reconstruction parameters on each validation task by fitting a linear regression model with parameters as inputs (independent variables) and performance metrics (MCC) as the output (dependent) variable. First, we built a feature matrix with one-hot encoded parameters for all reconstructed models. We then fit the linear regression using this feature matrix to predict the MCC value. This model’s coefficients for each parameter are then used to infer its impact on the predictive performance.

Gene essentiality predictions.

The results summarized on Fig 2 show that global thresholds have a greater impact on gene essentiality predictions, since the localT1 strategy, which places a greater emphasis on local thresholding leads to worse gene essentiality predictions. Additionally, the average of all models reconstructed using the global thresholding strategy is close to that found in localT2 models. Our best models obtained a MCC value of 0.44 while the best models from the work of Robinson as his colleagues was lower than 0.4, as represented on Fig 2A.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Overview of the influence of parameterization on the models’ performance when predicting essential genes as determined by their MCC.

For B and C panels, the essentiality threshold selected was -0.75. A: MCC value distribution for each CERES score threshold (horizontal axis) and algorithm combination (coloured box and whiskers). B: MCC value distribution for each thresholding strategy (horizontal axis) and algorithm combination (coloured box and whiskers). C: Linear coefficients for each individual parameter value on a regression model aimed at predicting MCC values. Each parameter variable was one-hot encoded as multiple binary variables.

https://doi.org/10.1371/journal.pcbi.1009294.g002

Despite this similarity when comparing the average of all models for each strategy, the best predictions were achieved using the localT2 strategy, which is a clear indicator that a combination of both thresholding approaches are useful to estimate RASs. Although the performance achieved using the localT2 strategy could be attributed to the usage of two (rather than one) global thresholds, we observed that the g_min parameter when using the localT2 strategy seems to have little impact on the models’ predictive power. This supports the claim that when both local and global strategies are combined, they provide a positive effect on the TAS calculation strategy, since localT1 does not show good results, and the ones from the global ones are improved when coupled with the localT2 strategy. The complete results are provided on S1 File.

We were also able to infer some of the properties associated with the dataset, where g_max and local thresholds at the 25th percentile seem to have the most positive effect on predictive ability in all models. Although the CERES score threshold representing the median essential gene knockout was set at -1, our models show slightly increased predictive power at -0.75.

Aside from data preprocessing related parameters, we have found the best parameter combinations are to use the tINIT algorithm in conjunction with the maximum function as replacement for the AND operator. We also observed that refining the model with EFMGapfill to allow growth using only the defined growth media metabolites as substrate did not result in better gene essentiality predictions.

Flux activity predictions.

We also performed a similar assessment on the ability of our models to correctly predict reaction activity/inactivity and directionality for the MCF7 cell line under three growth medium compositions with associated fluxomics measurements. Note that MCC is used given that the predicted variable is discrete. We did not address a quantitative prediction, as the preliminary results showed poor results in directly predicting flux values.

For each parameter combination, we generated three corresponding predictions using the growth medium as an additional flux constraint and calculated the MCC between the measured and predicted flux activities. In Fig 3A, we can see that most parameters affect flux and essentiality predictions similarly. The best performing strategies are still those based on global and local thresholding with 2 states. In both algorithms, we also observed that constraining nutrient uptake to the metabolites that could be matched with the growth medium led to higher predictive power.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Overview of the influence of parameterization on the models’ performance when predicting flux activity and using MCC as the evaluation metric.

Top (A): Linear coefficients for each individual parameter value on a regression model aimed at predicting MCC values for each parameter combination. Each parameter variable was one-hot encoded as multiple binary variables. Bottom (B): Relationship between average MCC value and standard deviation for each group of 3 simulations (conditions) that make up a single parameter combination. Different colors represent different algorithms and/or baseline comparison models.

https://doi.org/10.1371/journal.pcbi.1009294.g003

The relationship between average MCC and its standard deviation is shown in Fig 3B, where we can firstly observe that the flux predictions are, in general, of poor quality, with MCC values near 0. FASTCORE reconstructed models were able to reach the highest correlation with the experimental fluxomics data, although they are more sensitive to parameterization. INIT, on the other hand, showed higher average MCC values across all parameter combinations with less dispersion and yielded models that rank closer when evaluated with this metric. We also compared our models with a baseline MCF7 cell line model featured in the work of Robinson et al. [19], which ranks significantly lower than our best FASTCORE and INIT models. The data used can be seen on the S2 File.

FASTCORE appears as the more consistent tool to extract context-appropriate sets of reactions from a template model. Due to its low computational demand, these reconstructions can be repeated with alternative parameters or to sample large amounts of models. tINIT, on the other hand, shows consistently lower sensitivity to the parameters selection, as see on Figs 2A, 2B and 3B, leading to poorer performing models when compared to the best ones from FASTCORE.

Exploring metabolic variability in breast cancer.

We used the models and their respective predicted fluxes to explore the metabolic heterogeneity among various breast cancer cell lines. To do so, we retrieved molecular subtype annotations from the DepMap repository and used PCA to project these flux distributions using reduced features and obtain relevant information on the flux patterns present in different breast cancer subtypes. These results are summarized on Fig 4.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Cell line models reconstructed for breast cancer cell lines and projected in lower dimensions through the usage of PCA.

The left figure shows the first PC against the second, while the right figure displays the second PC against the third, in the horizontal and vertical axes, respectively.

https://doi.org/10.1371/journal.pcbi.1009294.g004

The subset of models belonging to breast cancer were extracted and the 200 fluxes (filtered by ANOVA, before the PCA) with most statistically significant differences among subtypes were used to decompose the dataset. We included all reconstructed models, regardless of their parameters, and reduced the latent space to 3 principal components (PCs) capable of explaining 33% of the observed variance. The loadings of each principal component were obtained, with each flux being summarised in their corresponding pathways by calculating the average of the absolute ratios between the weights of each flux and the maximum observed values in the PC loadings.

Firstly, we can conclude that the decomposition of predicted fluxes can adequately distinguish between major molecular subtypes of breast cancer. The second PC (PC2) marks a good distinction between basal and luminal cell lines and correlates with reported prognosis and aggressiveness [53], with luminal BCs with better prognosis assigned to positive values as opposed to basal BCs which appear in this PC as negative values.

Predicted metabolic fluxes as relevant features.

The lack of fluxomics data for the whole set of cell lines featured in DepMap does not allow to carry out a large-scale systematic comparison of the pFBA flux distribution predictions with experimental data. However, we set out to assess whether or not these predicted fluxes could be useful in predicting several clinical features associated with each sample. To do so, we established a supervised classification task, where the disease’s primary location would be predicted using various datasets as inputs, namely, (1) standardized expression values (TPM from RNASeq) using the entire gene set, as well as only those genes that can be integrated in the metabolic model, (2) TASs generated for each sample, (3) predicted fluxes (using pFBA over reconstructed models) and (4) the reaction presence (binary) from the CSMR algorithm outputs. Our results, using random forests as the machine learning model to address these tasks, are summarized on Fig 5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Distribution of average Matthews’ correlation coefficient values for cross-validated classifiers trained with various datasets.

https://doi.org/10.1371/journal.pcbi.1009294.g005

Classifiers trained with standardized transcriptomics data showed good relative performance (MCC mean = 0.570, sd = 0.038) with the subset corresponding to metabolic genes only slightly outperforming it. Processing these data and generating TASs slightly increased predictive capabilities (MCC mean = 0.594, sd = 0.030), which further justifies applying our preprocessing workflows before analysing and integrating omics data. However, the outputs of CSMR algorithms, namely, the presence or absence of each reaction resulted in models capable of predicting a cancer cell line’s primary site with an average MCC of 0.525 (sd = 0.031). Furthermore, pFBA simulations resulted in even worse classifiers that could only reach an average MCC of 0.298 (sd = 0.042).

Overall, our results show that context-specific model reconstruction and flux balance analysis approaches are not yet consistent enough for accurate quantitative flux predictions, as predicted metabolic fluxes by themselves did not appear to be relevant features for complex classification tasks.