
==== Front
STAR Protoc
STAR Protoc
STAR Protocols
2666-1667
Elsevier

S2666-1667(24)00456-8
10.1016/j.xpro.2024.103291
103291
Protocol
Protocol for genome-scale differential flux analysis to interrogate metabolic differences from gene expression data
Beura Satyajit 14
Das Amit Kumar 1
Ghosh Amit amitghosh@iitkgp.ac.in
235∗
1 Department of Bioscience and Biotechnology, Indian Institute of Technology Kharagpur, West Bengal 721302, India
2 School of Energy Science and Engineering, Indian Institute of Technology Kharagpur, West Bengal 721302, India
3 P.K. Sinha Centre for Bioenergy and Renewables, Indian Institute of Technology Kharagpur, West Bengal 721302, India
∗ Corresponding author amitghosh@iitkgp.ac.in
4 Technical contact

5 Lead contact

04 9 2024
20 9 2024
04 9 2024
5 3 103291© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

Deciphering the functional differences between diseased and healthy cells requires understanding the alterations in biochemical flux patterns. We present a genome-scale differential flux analysis (GS-DFA) protocol to elucidate these metabolic disparities by integrating condition-specific gene expression data into the human genome-scale metabolic model (humanGEM). In this protocol, we describe the steps to normalize and integrate data into the humanGEM and analyze differential flux across the biochemical network between diseased and healthy cells.

For complete details on the use and execution of this protocol, please refer to Nanda et al.1

Graphical abstract

Highlights

• Installation guide for necessary software and libraries for GS-DFA

• Instructions for preprocessing and normalizing the gene expression data

• Steps for integrating the normalized gene expression data into the humanGEM

• Detailed procedures for GS-DFA between the diseased cell and control

Publisher’s note: Undertaking any experimental protocol requires adherence to local institutional guidelines for laboratory safety and ethics.

Deciphering the functional differences between diseased and healthy cells requires understanding the alterations in biochemical flux patterns. We present a genome-scale differential flux analysis (GS-DFA) protocol to elucidate these metabolic disparities by integrating condition-specific gene expression data into the human genome-scale metabolic model (humanGEM). In this protocol, we describe the steps to normalize and integrate data into the humanGEM and analyze differential flux across the biochemical network between diseased and healthy cells.

Subject areas

bioinformatics
RNA-seq
systems biology
==== Body
pmcBefore you begin

Metabolic diseases encompass a broad spectrum of disorders that affect the ability of a cell to process nutrients, leading to disruptions in essential biochemical pathways.2 These metabolic alterations pose significant challenges to human health. Exploring the metabolic alterations at a genome-scale level will provide profound insights into the metabolic dysregulations of diseased cells compared to healthy controls. In this context, condition-specific genome-scale metabolic modeling is introduced to capture the altered flux distributions within the biochemical pathways of specific cells or tissues under different conditions.3,4,5,6,7,8 The condition-specific GEMs are reconstructed by integrating gene expression data with human metabolic models. These models offer a powerful framework for understanding the dynamic interplay between gene regulation and metabolic activity within specific biological conditions. Various highly curated and experimentally validated humanGEMs, including Recon3D9 and HumanGEM,10 have been developed as foundational frameworks, which enclose all the metabolic reactions in human. To make a tissue/condition-specific model, it becomes essential to constraint the reactions in the humanGEM based on cell/tissue-specific gene expression data.1,11 Algorithms such as iMAT,12 INIT,13 tINIT14 etc. have empowered the integration of gene expression data into humanGEM models, facilitating the reconstruction of tissue-, cell-, or condition-specific GEMs. The reconstructed condition-specific GEMs can determine the extent of alterations in metabolic reaction fluxes, allowing us to comprehend the activities of various biochemical pathways between diseased and healthy cells. These altered metabolic activities can serve as biomarkers for disease detection and aid in developing therapeutics. In this protocol, we introduce a methodology for GS-DFA, aimed at interrogating the metabolic variances derived from condition-specific gene expression data. It includes instructions for installing essential software, preprocessing and normalizing data, and integrating it into the HumanGEM framework to create condition-specific models. Finally, the protocol outlines meticulous procedures for analyzing genome-scale differential flux across the biochemical network, employing the condition-specific models. This analytical phase bridges RNA-Seq data with cellular insights, elucidating altered metabolic activities across diverse conditions. The insights can offer valuable perspectives into the intricate metabolic disturbances that play a pivotal role in cellular homeostasis.

Installing relevant libraries and software

Timing: 2 h

1. Install R and RStudio for preprocessing of raw RNA-seq data.

Note: Ignore this step if you already have TPM normalized gene expression data

2. Install MATLAB 2024a from MATLAB: https://matlab.mathworks.com/.

Note: We have tested model reconstruction with MATLAB releases from 2022a to 2024a.

3. Install Anaconda from Anaconda: https://www.anaconda.com.

Note: Launch Jupyter Notebook/Spider in the Anaconda Navigator to perform the model analysis.

4. Install COBRA Toolbox and COBRApy as per instructions from GitHub: https://opencobra.github.io/.

5. Install Git Bash from Git: https://git-scm.com/downloads.

6. Install Gurobi Optimization Solver from Gurobi: https://www.gurobi.com.

Note: Use your academic email ID to get it for free

7. Install RAVEN 2.4.0 and HumanGEM 1.4.1 from GitHub: https://github.com/SysBioChalmers.

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Deposited data	
	
Raw code	GitHub	https://github.com/itsamit/GS-DFA	
RNA-seq data	NCBI	https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE153315	
Transcript length	MTD	https://ngdc.cncb.ac.cn/mtd/mRNA.php	
HumanGEM 1.14.0	Zenodo	https://zenodo.org/records/10303455	
	
Software and algorithms	
	
R 4.3.3	RStudio	https://posit.co/download/rstudio-desktop/	
RStudio	RStudio	https://posit.co/download/rstudio-desktop/	
Anaconda	Anaconda	https://www.anaconda.com/download	
MATLAB 2024a	MATLAB	https://matlab.mathworks.com/	
COBRA Toolbox	Heirendt et al.15	https://opencobra.github.io/	
COBRApy	Ebrahim et al.16	https://opencobra.github.io/	
RAVEN 2.4.0	SysBio.se	https://github.com/SysBioChalmers	
	
Other	
	
Gurobi 10.0.1 (Solver)	Gurobi	https://www.gurobi.com	

Step-by-step method details

Convert gene expression to TPM values

Timing: 1 h

TPM is the most accepted normalization method for RNA sequencing data. The first step is to convert raw counts to normalized values.1. Download the gene expression data from GEO.

Note: The raw read counts of RNA-seq data are available for download from the NCBI database. One can go to NCBI home page, select GEO Datasets and search for a specific study by the GEO ID or using Keyword. In this protocol, we have considered the gene expression data (GEO: GSE153315) from 30 samples, including metformin-responder T2D patients, metformin non-responder T2D patients, and a healthy control group.

Searching the GEO dataset in NCBI

Downloading the GEO dataset

Raw data set

2. Download the Transcript lengths from MTD database.

Note: Transcript lengths vary depending on which sets of exons undergo alternative splicing. It is important to note that normalizing read counts agnostic of specific transcript isoforms can change the TPM values (which is a function of length) and consequently downstream calculations. To obtain the tissue-specific gene length, use the Mammalian Transcriptomic Database (MTD), navigate to the download section, and select the desired cell line. The downloaded file provides the information about the starting and ending nucleotide positions of a specific gene in the genome. By subtracting the start position of the gene from the end position, the gene length can be determined. As in our case, we were interested in transcripts from blood cell, we obtained the transcript lengths for the specific tissue from MTD: https://ngdc.cncb.ac.cn/mtd/.17

Downloading the gene length from MTD

CRITICAL: Save both the RNA sequencing data and Transcript lengths data file in a in the same folder (TPM_normalization) where the data preprocessing code (id_change_RNA_sequencing_data.R, id_change_Transcript_lengths_data.R and Final_input_file_for_TPM.R) and normalization code TPM.R is present.

3. Convert gene symbols to ensemble IDs in RNA sequencing data file.

Note: Execute the scripts id_change_RNA_sequencing_data.R to convert gene symbols to ensemble IDs in RNA sequencing data file.

setwd("C:/Users/DELL/Desktop/Tpm_Norm") # Path

Data1 <- read.csv("GSE153315_counts.csv") # Import the file

library(hgu95av2.db) # Import hgu95av2.db

keytypes(hgu95av2.db) # Converting Gene Symbol to ensembl IDs

ensembl.id <- mapIds(hgu95av2.db, keys= Data1$Gene, column='ENSEMBL', keytype='SYMBOL', multiVals='first')

ensembl.id <- as.matrix(ensembl.id)

write.csv(ensembl.id, "GSE153315_counts_ENSEMBL.csv")

Data2 <- read.csv("GSE153315_counts_ENSEMBL.csv") # Post processing of the file

colnames(Data2)[1]='Gene'

merged_df <- merge(Data2, Data1, by = "Gene", all = TRUE)

print(merged_df)

df_cleaned <- na.omit(merged_df)

print(df_cleaned)

df_without_first_column <- df_cleaned[, -1]

print(df_without_first_column)

colnames(df_without_first_column)[1]='Gene_Symbol'

write.csv(df_without_first_column, "GSE153315_counts_ENSEMBL_Final.csv")

4. Convert gene symbols to ensemble IDs in Transcript lengths data file.

Note: Execute the scripts id_change_Transcript_lengths_data.R to convert gene symbols to ensemble IDs in Transcript lengths data file.

setwd("C:/Users/DELL/Desktop/Tpm_Norm") # Path

Data1 <- read.csv("SRX105932_gene_length.csv") # Import the file

library(hgu95av2.db) # Import hgu95av2.db

keytypes(hgu95av2.db)

ensembl.id <- mapIds(hgu95av2.db, keys= Data1$Gene_Symbol, column='ENSEMBL', keytype='SYMBOL', multiVals='first')

ensembl.id <- as.matrix(ensembl.id)

write.csv(ensembl.id, "SRX105932_gene_length_ENSEMBL.csv")

Data2 <- read.csv("SRX105932_gene_length_ENSEMBL.csv")

colnames(Data2)[1]='Gene_Symbol'

merged_df <- merge(Data2, Data1, by = "Gene_Symbol", all = TRUE)

print(merged_df)

df_cleaned <- na.omit(merged_df)

print(df_cleaned)

df_without_first_column <- df_cleaned[, -1]

print(df_without_first_column)

colnames(df_without_first_column)[1]='Gene_Symbol'

write.csv(df_without_first_column, "SRX105932_gene_length_ENSEMBL_Final.csv")

5. Combine the gene expression data and there corresponding transcript length in a single file.

Note: Execute the scripts Final_input_file_for_TPM.R and save the ‘TPM_Input_File.csv’ file.

setwd("C:/Users/DELL/Desktop/Tpm_Norm") # Path

Data1 <- read.csv("SRX105932_gene_length_ENSEMBL_Final.csv")

df_without_first_column_D1 <- Data1[, -1]

Data2 <- read.csv("GSE153315_counts_ENSEMBL_Final.csv")

df_without_first_column_D2 <- Data2[, -1]

merged_df <- merge(df_without_first_column_D2, df_without_first_column_D1, by = "Gene_Symbol", all = TRUE)

print(merged_df)

df_cleaned <- na.omit(merged_df)

write.csv(df_cleaned, "TPM_Input_File.csv")

Output of raw data processing:

6. Perform TPM Normalization (R script).

library(dplyr) # Load necessary library

data <- read.csv("TPM_Input_File.csv") # Load your dataset

gene_names <- data[, 2] # Extract gene names from the first column

expression_values <- data[, 3:32] # Extract expression values from columns 2-31

gene_lengths <- data[, ncol(data)] # Gene lengths are in the last column

gene_lengths_kbp <- gene_lengths / 1000

Read_per_kelobase <- sweep(expression_values, 1, gene_lengths_kbp, "/")

Normalization_factor <- colSums(Read_per_kelobase) / 1000000

TPM <- sweep(Read_per_kelobase, 2, Normalization_factor, "/")

TPM <- cbind(genes = gene_names, TPM) # Add gene names as the first column to the TPM data frame

print(TPM) # Print the resulting data frame with TPM values

write.csv(TPM,"TPM_File.csv",row.names = FALSE) # Save the TPM file

Output of TPM normalization

Note: To verify the correctness of the normalization data, sum all the gene expression values of a sample. If the total is 1,000,000, then the normalization data is correct.

Reconstruction of condition-specific genome-scale metabolic model

Timing: 7 h

Condition-specific HumanGEM refers to an in silico model that represents all the metabolic reactions occurring in specific tissues under certain conditions. These condition/tissue-Specific models focus on the unique metabolic activities of particular tissues (like liver, muscle, or brain) or specific physiological or pathological conditions (such as fasting, exercise, or disease states). In this section, we added the method and script to preprocess tissue-specific gene expression data and incorporate it into the HumanGEM to reconstruct condition-specific models. This helps us better understand how metabolism varies across different tissues and conditions.7. Perform input data (TPM Normalized) preprocessing before integrating into the model.

Note: The HumanGEM model encompasses approximately 3000 metabolic genes. However, the extracted RNA-seq (TPM Normalized) data comprises a vast number of genes, in our study it is 12000. Consequently, prioritizing the integration of the metabolic genes into the model is imperative. To extract the gene expression data specifically for metabolic genes, save the TPM normalized file, the HumanGEM gene list file, and the input data processing script in the 'Model_reconstruction' folder, then execute the script.

import pandas as pd # import the library

df1 = pd.read_csv('Human-GEM_genes.csv') # read csv data

df2 = pd.read_csv('TPM_normalized_data.csv') # read csv data

inner_join = pd.merge(df1, df2, on ='genes', how ='inner') # merge the data files

column_index_to_drop = 1 # drop Unnamed column

df = inner_join.drop(inner_join.columns[column_index_to_drop], axis=1)

sample_groups = [['s1', 's3', 's7', 's8', 's9', 's12', 's13', 's14', 's18', 's20'], ['s2', 's4', 's5', 's6', 's10', 's11', 's15', 's16', 's17', 's19'], ['s21', 's22', 's23', 's24', 's25', 's26', 's27', 's28', 's29', 's30']] # Grouping samples

average_values = [] # Calculate average values for each group

for group in sample_groups:

 selected_samples = df[group]

 average_values.append(selected_samples.mean(axis=1))

new_df = pd.DataFrame(average_values).T # Create a new DataFrame to save the average values

new_df.columns = ['non-responder_T2D', 'responder_T2D', 'control']

df1=df.genes # Save the file

new_df1 = pd.DataFrame(df1)

joined_df = pd.concat([new_df1, new_df], axis=1)

joined_df.set_index(joined_df.columns[0], inplace=True)

joined_df.to_csv('TPM_normalized_data_final.csv')

Output of processed TPM normalization:

8. Integration of gene expression data into the HumanGEM model.a. Download MATLAB, Git Bash, and Gurobi software, and proceed to install them on your system.

b. Download the COBRA Toolbox, RAVEN, and HumanGEM model and add it to the MATLAB path by clicking “add with subfolders” and save it.

c. Further, go to “Bourse the folder” section of the MATLAB console and add the necessary files like: TPM normalized data (processed data), essential metabolic tasks, HAM media, and model integration script.

Tools added to MATLAB path:

d. Then, run the script i.e., Model_Reconstruction_Code_WithCuration.mat to get the condition-specific model.Note: The successful completion of the data integration process will generate condition-specific models. For our case it is control, responder_T2D, and non_responder_T2D model.

i. Import TPM normalized data and HumanGEM model for data integration into the model.initCobraToolbox();

changeCobraSolver('gurobi','all');

setRavenSolver('gurobi');

data = readtable('TPM_normalized_data_final.csv');

data_struct.tissues = data.Properties.VariableNames(2:end)'; % sample (tissue) names

data_struct.genes = data.genes; % gene names

data_struct.levels = table2array(data(:, 2:end)); % gene TPM values

data_struct.threshold = 1;

load('Human-GEM.mat');

ihuman = addBoundaryMets(ihuman);

ii. Check the essential metabolic functions of model before integration.essentialTasks = parseTaskList('metabolicTasks_Essential_xlsx.xlsx');

checkTasks(ihuman, [], true, false, false, essentialTasks);

refModel = ihuman; % the reference model from which the GEM will be extracted

celltype = []; % used if tissues are subdivided into cell type, which is not the case here

hpaData = []; % data structure containing protein abundance information (not used here)

arrayData = data_struct; % data structure with gene (RNA) abundance information

metabolomicsData = []; % list of metabolite names if metabolomic data is available

removeGenes = true; % (default) remove lowly/non-expressed genes from the extracted GEM

taskFile = []; % we already loaded the task file, so this input is not required

useScoresForTasks = true; % (default) use expression data to decide which reactions to keep

printReport = true; % (default) print status/completion report to screen

taskStructure = essentialTasks; % metabolic task structure (used instead "taskFile")

params.TimeLimit = 5000; % additional optimization parameters for the INIT algorithm

paramsFT = []; % additional optimization parameters for the fit-tasks algorithm

iii. Integrate the gene expression data into the HumanGEM model and optimize the model in the HAM media.for i = 2:4

 tissue = char(data.Properties.VariableNames(i));

 disp(tissue);

 model = getINITModel2(refModel, tissue, celltype, hpaData,

 arrayData, metabolomicsData, removeGenes,

 taskFile, useScoresForTasks, printReport,

 taskStructure, params, paramsFT);

 model.id = tissue;

 model = simplifyModel(model);

 model.b = repelem(0,length(model.b))';

 media=readtable('Trial_HAM_HumanGem1.14_csv.csv');

 metList=table2array(media(:,8));

 lb=table2array(media(:,9));

 exc=model.rxns(findExcRxns(model));

 mediacomps=exc(ismember(exc,findRxnsFromMets(model,metList)));

 model=changeRxnBounds(model,exc,0,'l'); %Constrain by HAM media

 model=changeRxnBounds(model,mediacomps,lb,'l');

 writeCbModel(model, 'fileName' , tissue , 'format','mat');

 pause(900);

end

Successful execution of the script will print the final model statistics:

Note: For model optimization, either Gurobi or CPLEX can be used. However, we have chosen Gurobi for this protocol.

Note: The integration of transcriptomics data into the HumanGEM model using the tINIT algorithm is not the only option. Other tools, such as iMAT and GIMME, can also be used.

Model characteristics

Timing: 30 min

The primary analysis for the condition-specific models (Control, Responder_T2D, and Non_responder_T2D GEMs) involves assessing its characteristics. This includes evaluating the number of activated genes, metabolites, and reactions within the respective models.9. Save the models and the script in a dedicated folder for the evaluation of model characteristics, and proceed to execute the script. A successful execution of the code will generate the model characteristics plot (Figure 1).a. Import the models and print the number of active reactions, metabolites and genes from each model.import sys # import the libraries

import os

import cobra

modelCON=cobra.io.load_matlab_model('control.mat') # import the models

modelT2D_R=cobra.io.load_matlab_model('responder_T2D.mat')

modelT2D_NR=cobra.io.load_matlab_model('non_responder_T2D.mat')

# Characteristics of control model

print(len(modelCON.genes), len(modelCON.reactions), len(modelCON.metabolites))

# Characteristics of responder_T2D model

print(len(modelT2D_R.genes), len(modelT2D_R.reactions), len(modelT2D_R.metabolites))

# Characteristics of non_responder_T2D model

print(len(modelT2D_NR.genes), len(modelT2D_NR.reactions), len(modelT2D_NR.metabolites))

b. Use the number of active reactions, metabolites and genes to plot model characteristics.import numpy as np

import matplotlib.pyplot as plt

# Data for three samples, each containing counts of genes, reactions, and metabolites

samples = ['Control', 'Responder_T2D', 'Non_responder_T2D']

genes = [1803, 1788, 1753] # Gene counts for each sample

reactions = [6232, 6054, 6071] # Reaction counts for each sample

metabolites = [4460, 4405, 4430] # Metabolite counts for each sample

bar_width = 0.25 # Set the width of the bars

r1 = np.arange(len(samples)) # Set the position of the bars on the x-axis

r2 = [x + bar_width for x in r1]

r3 = [x + bar_width for x in r2]

# Plotting the grouped bar chart

plt.bar(r1, genes, color='b', width=bar_width, edgecolor='grey', label='Genes')

plt.bar(r2, reactions, color='g', width=bar_width, edgecolor='grey', label='Reactions')

plt.bar(r3, metabolites, color='r', width=bar_width, edgecolor='grey', label='Metabolites')

plt.xlabel('Samples', fontweight='bold') # Add xticks on the middle of the group bars

plt.xticks([r + bar_width for r in range(len(samples))], samples)

plt.legend() # Adding a legend and title

plt.title('Model Characteristics')

plt.savefig('bar_plot.png') # Save the plot in the working directory

plt.show() # Show plot

Output

Note: In this protocol, we have considered the condition-specific transcriptomics data, such as control, T2D Metformin responder, and T2D Metformin non-responder, to reconstruct the condition-specific model. The reconstructed models were named Control, Responder_T2D, and Non_responder_T2D, respectively. The reconstructed models, underwent evaluation of model characteristics, revealing the number of genes, active reactions, and metabolites in each model (Figure 1). Moreover, these models were also used to study the differential flux profiles among the three different conditions in the subsequent sections.

Note: In addition to gene, active reaction, and metabolite information, one can also print media components, exchange reactions, and model growth rate etc. Refer to the COBRApy documentation (openCOBRA: https://cobrapy.readthedocs.io/en/latest/) for scripts to execute various characteristics of the reconstructed model.

Figure 1 Model characteristics of healthy and two different types of diseased cells, blue represents genes, green represents reactions, and red represents metabolites present in different models

Perform flux-sampling using ACHR

Timing: 3 h

Flux sampling is a computational technique crucial in analyzing the flux profile of genome-scale metabolic models as well as the community scale metabolic models.18 By randomly sampling the solution space, flux sampling generates sets of feasible flux distributions that adhere to model constraints and objectives.1,18 Overall, flux sampling is a powerful tool for gaining insights into genome-scale metabolic networks, facilitating their manipulation for biotechnological and biomedical purposes. In this section, we have added the script to perform flux sampling using ACHR method.10. Save the models and the script for flux sampling in a dedicated folder, and proceed to execute the script. A successful execution of the code will generate flux sampling files.

import cobra # import the libraries

from cobra.sampling import ACHRSampler

modelCON=cobra.io.load_matlab_model('control.mat') # import the models

modelT2D_R=cobra.io.load_matlab_model('responder_T2D.mat')

modelT2D_NR=cobra.io.load_matlab_model('non_responder_T2D.mat')

achr_modelCON = ACHRSampler(modelCON, thinning=100) # performed flux sampling

achr_modelT2D_R= ACHRSampler(modelT2D_R,thinning=100)

achr_modelT2D_NR = ACHRSampler(modelT2D_NR, thinning=100)

samples_modelCON=achr_modelCON.sample(10000)

samples_modelT2D_R=achr_modelT2D_R.sample(10000)

samples_modelT2D_NR=achr_modelT2D_NR.sample(10000)

samples_modelCON.to_csv("modelCON.csv") # save the FS files

samples_modelT2D_R.to_csv("modelT2D_R.csv")

samples_modelT2D_NR.to_csv("modelT2D_NR.csv")

Note: Our methodology employed Artificial Coordinate Hit and Run (ACHR) sampler for sampling fluxes of a community model.19 To ensure sparsity and uncorrelatedness within the fluxes, a thinning factor of 100 was applied, facilitating comprehensive coverage of the flux space. The optimization process guiding ACHR flux sampling was conducted for each model, generating 10,000 sample points using the sampling object. The flux sampling script has been coded in Python.

Note: The “thinning factor” in flux sampling refers to a parameter used to reduce autocorrelation in the generated samples by selecting only one value in every nth (in this protocol, n = 100) samples. This helps to ensure that the retained samples are more representative of the true distribution by reducing the dependency between consecutive samples. Increasing the thinning factor (n > 100) can more effectively reduce autocorrelation, potentially leading to better results in terms of sample independence. However, an increase in the thinning factor also escalates computational costs.

Perform differential flux analysis

Timing: 30 min

In this section, the sample flux calculated from the condition-specific models have been used to determine the differential flux through metabolic reactions associated with a biochemical pathway between healthy and diseased cell/tissue. A generalized figure (Figure 2) has been added to illustrate how the differential flux analysis data can be mapped.11. Save the flux sampling files and the script for differential flux analysis in a dedicated folder, and proceed to execute the script. A successful execution of the code will generate differential flux files.a. Import libraries and flux sampling data files.import cobra # import the libraries

import matplotlib.pyplot as plt

import scipy as sp

import numpy as np

import pandas as pd

from sklearn import decomposition

from sklearn import datasets

from sklearn.preprocessing import scale

from sklearn.linear_model import LinearRegression

import statsmodels

from scipy.stats import sem, t

from scipy import mean

from statsmodels.sandbox.stats.multicomp import multipletests

import seaborn as sns

from scipy.stats import hypergeom

samples_HC = pd.read_csv("modelCON.csv") # import FS files

samples_T2D_R = pd.read_csv("modelT2D_R.csv")

samples_T2D_NR = pd.read_csv("modelT2D_NR.csv")

b. Perform differential flux analysis between control and T2D responder.def kstest(samplesUI,samplesI,file_name): # differential flux analysis

 rxns1=set(samplesUI.columns)

 rxns2=set(samplesI.columns)

 rxn_c=rxns1.intersection(rxns2)

 pvals=[]

 rxnid=[]

 fc=[]

 for rxn in rxn_c:

 data1=samplesUI[rxn].round(decimals=4)

 data2=samplesI[rxn].round(decimals=4)

 data1=data1.sample(n=1000)

 data2=data2.sample(n=1000)

 if((data1.std()!=0 and data1.mean()!=0) or (data2.std()!=0 and data2.mean()!=0)):

 kstat,pval=sp.stats.ks_2samp(data1,data2)

 foldc=(data1.mean()-data2.mean())/abs(data1.mean()+data2.mean())

 pvals.append(pval)

 rxnid.append(rxn)

 fc.append(foldc)

 data_mwu=pd.DataFrame({'Reaction':rxnid,'Pvalue':pvals})

 data_mwu=data_mwu.set_index('Reaction')

 plt.hist(data_mwu['Pvalue'],100)

 plt.xlabel('P-value')

 plt.ylabel('Frequency')

 plt.title('P-value distribution')

reject,padj,_,_=statsmodels.stats.multitest.multipletests(data_mwu['Pvalue'], alpha=0.05, method='fdr_bh', is_sorted=False, returnsorted=False)

 data_mwu['Padj']=padj

 data_mwu['Reject']=reject

 data_mwu['FC']=fc

 data_sigFC=data_mwu.loc[(abs(data_mwu['FC'])>0.82) & (data_mwu['Padj']<0.05),:]

 file=file_name+"_data_sigFC.csv"

 tile=file_name+"_mwu.csv"

 data_sigFC.to_csv(file)

 data_mwu.to_csv(tile)

kstest(samples_HC,samples_T2D_R,'T2D_R')

Note: For calculating the differential flux analysis between control and T2D non responder replace “samples_T2D_R″ with “samples_T2D_NR” in Kstest (Last line of the script).

Note: In differential flux analysis, the flux change (FC) cutoff value can be adjusted to identify significantly impacted metabolic reactions.

Figure 2 A schematic diagram illustrates the differential flux analysis

The differential flux map displays the flux changes (FC) through metabolic reactions, with green arrows indicating increased flux, red arrows indicating decreased flux, and black arrows showing no change. The width of the arrow indicates the flux change value.

Generate potential pathways altered between control and experimental conditions

Timing: 1 h

Hypergeometric enrichment analysis can performed to identify subsystems within the HumanGEM model that are overrepresented in the set of altered reactions. Precisely, the augmented list of reactions obtained from above (differential flux analysis) can employed to conduct pathway enrichment analysis (Figure 3).12. Save the differential flux file, subsystem file (Pathway information) and the script for pathway enrichment analysis in a dedicated folder, and proceed to execute the script.a. Import libraries, differential flux analysis data and subsystem file.import cobra # import libraries

import matplotlib.pyplot as plt

import scipy as sp

import numpy as np

import pandas as pd

from sklearn import decomposition

from sklearn import datasets

from sklearn.preprocessing import scale

from sklearn.linear_model import LinearRegression

import statsmodels

from scipy.stats import sem, t

from scipy import mean

from statsmodels.sandbox.stats.multicomp import multipletests

import seaborn as sns

from scipy.stats import hypergeom

%matplotlib inline

rxnlist=pd.read_csv('T2D_R_data_sigFC.csv') # Import differential flux analysis data

subs=pd.read_csv('HumanGem_Subsytems.csv') # Import subsystem file

b. Perform pathway enrichment analysis.dataset=pd.DataFrame() # Pathway enrichment analysis

for path in subs['Var2'].unique():

 reaction_set=subs.loc[subs['Var2']==path,'Var1']

 rxn=reaction_set.reset_index(drop=True)

 df_temp=pd.DataFrame({path:rxn})

 dataset=pd.concat([dataset,df_temp],axis=1)

listSize=len(rxnlist)

listrxnSize=[]

setSize=[]

rxnSize=13025

for col in dataset.columns:

 df=pd.DataFrame({'Reaction':dataset[col]})

 out=df.merge(rxnlist)

 listrxnSize.append(len(out))

 setSize.append(len(dataset[col].dropna()))

hyperdata=pd.DataFrame({'Pathways':dataset.columns,'ListReactions':listrxnSize,'SetSize':setSize})

hits=hyperdata['ListReactions']

pool=hyperdata['SetSize']

allrxns=hyperdata['SetSize'].sum()

targetrxns=hyperdata['ListReactions'].sum()

pvalList=[]

for h,p in zip(hits,pool):

rv=hypergeom(allrxns-p,p,targetrxns)

pval=rv.pmf(h)

pvalList.append(pval)

hyperdata['P-value']=pvalList

reject,padj,_,_=statsmodels.stats.multitest.multipletests(hyperdata['P-value'], alpha=0.05, method='fdr_bh', is_sorted=False, returnsorted=False)

hyperdata['P-valueadj']=padj

hyperdata['Reject']=reject

df = hyperdata.sort_values(by='ListReactions',ascending=False)

df.to_csv("T2D_R_Pathway_enrichment.csv")

c. Generate a plot using pathway enrichment data file.hyperdata_sig=hyperdata[(hyperdata['Reject']) & (hyperdata['ListReactions']!=0)]

hyperdata_sorted=hyperdata_sig.sort_values(by='P-valueadj',ascending=False)

plt.figure(figsize=(14,4))

sc=plt.scatter(hyperdata_sorted['P-valueadj'],np.arange(0,len(hyperdata_sorted['Pathways'])),s=hyperdata_sorted['ListReactions'],color=(0.7,0.2,0.9,0.7))

plt.xlabel('Adjusted p-value')

plt.yticks(np.arange(0,len(hyperdata_sorted['Pathways'])),labels=hyperdata_sorted['Pathways'])

handles, labels = sc.legend_elements(prop="sizes", alpha=0.8)

plt.legend(handles, labels, bbox_to_anchor=(1.6,1.02),loc='upper right',title="Reactions")

# plt.grid(axis='y')

plt.tight_layout()

plt.savefig('T2D_R_pic_sigFC.png',dpi=600)

Output

Note: A subsystem file containing pathway-level information is available within the HumanGEM model file. For pathway enrichment analysis, this subsystem file should be extracted from the HumanGEM file and saved in a dedicated folder, along with the python script and the differential flux file.

Figure 3 GS-DFA highlights the affected pathways in diseased cell enriched with altered metabolic reactions

(A) Affected metabolic pathways in metformin non-responder T2D patients.

(B) Affected metabolic pathways in metformin responder T2D patients. The presented reactions have a fold change (FC) cutoff of 10 and an adjusted p-value less than 0.05.

Expected outcomes

The condition-specific models for both disease and a healthy cell can be reconstructed by integrating gene expression data into the HumanGEM. Furthermore, the Flux sampling can be performed using the reconstructed model, followed by differential flux analysis and pathway enrichment analysis. Differential flux analysis can pinpoint alterations in metabolic flux through specific biochemical reactions within diseased cells compared to healthy ones. Further, the differential flux profile can be mapped onto biochemical pathways through pathway enrichment analysis. This analysis can unveil disruptions in cellular metabolism within diseased cells compared to healthy controls (Figure 4). The GS-DFA offers a powerful approach to dissect the complex metabolic network of human cell and elucidate the aberrations associated with disease. This approach can also pave the way for developing targeted therapeutic interventions aimed at restoring normal cell function and improve metabolism of the patient.Figure 4 GS-DFA reveals the correlation between the metabolic pathways in diseased cells and their impact on the patients

Green arrow indicates the upregulated metabolic pathways whereas the red arrow indicated the downregulated pathways.

Limitations

Data quality

The reliability and quality of gene expression data can vary significantly depending on experimental conditions, techniques, and sample sources. Integration of heterogeneous data sources can introduce noise and biases, potentially impacting the accuracy of the metabolic model predictions.

Dynamic nature of metabolism

Metabolic networks are highly dynamic and can adapt to various environmental stimuli and cellular states. Gene expression data provide a snapshot of cellular activity at a specific moment, but they may not capture the full spectrum of metabolic dynamics, particularly for fast-changing processes or transient responses.

Missing metabolic reactions and pathways

Even comprehensive metabolic models like HumanGEM may lack certain reactions or pathways due to gaps in current knowledge or incomplete annotations. This can limit the model’s ability to accurately simulate cellular metabolism, especially in specialized or poorly characterized cell types or conditions.

Troubleshooting

Problem 1

Preprocessing of the RNA-seq data (Related to Step 3).

Potential solution

The first step in GS-DFA is TPM normalization, so preprocess the RNA-Seq data is crucial. To normalize the data the file should consist of raw read counts and a uniform gene ID. Given the variance in gene annotation across different RNA-Seq datasets, standardizing it to the ENSEMBL ID format becomes necessary, a task facilitated by the widely utilized hgu95av2.db library. This format alignment is essential as the HumanGEM model operates on the ENSEMBL ID format. Notably, hgu95av2.db encompasses a comprehensive collection of 27 gene ID formats. Therefore, before selecting the data, matching the gene IDs is imperative.

Problem 2

If preprocessed gene expression data (normalized) is available in NCBI instead of raw read count data. (Related to Step 6).

Potential solution

In some cases, instead of raw count data, preprocessed gene expression data (normalized) is available in the NCBI repository. In such instances, the processed data can be directly used for integrating transcriptomics data into the HumanGEM model.

Problem 3

Working directory. (Related to all the major steps i.e., Steps 3–12).

Potential solution

Since we’re using three distinct programming languages i.e., R, Python, MATLAB, for the successful completion of the project, it’s essential to set the working directory before executing the code.

Problem 4

Preparing the MATLAB environment for integrating the gene expression data into the HumanGEM model. (Related to Step 8).

Potential solution

First, install MATLAB, Gurobi, and Git Bash on your system, then proceed to follow the instructions below.

For the easy execution of the script, we have enclosed the HumanGEM, TPM normalized data, essential metabolic tasks, HAM media, and model integration script with this manuscript. The COBRA Toolbox, RAVEN, and HumanGEM model need to be downloaded and added in the MATLAB path by clicking “add with subfolders” and save it. Further, go to “Browse the folder” section and add the necessary files like: TPM normalized data, essential metabolic tasks, HAM media, and model integration script to MATLAB. Then, run the script to get the condition-specific models.

Problem 5

Preparing media file. (Related to Step 8).

Potential solution

In this protocol, we have used HAM media. However, the media components will vary depending on the cell/tissue and experiment. In such cases, generate a separate media file by following the format of the HAM media file included in this manuscript. Then, add it to the working directory before running the script.

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Prof. Amit Ghosh (amitghosh@iitkgp.ac.in).

Technical contact

For technical information, including software and code, contact Satyajit Beura (satyajit_beura@kgpian.iitkgp.ac.in).

Materials availability

This is a computational biology pipeline, so no biological materials or reagents were used.

Data and code availability

The authors declare that all other relevant data are available within the article and its supplementary data files. The associated models and codes are available in GitHub (https://github.com/itsamit/GS-DFA) and Zenodo (https://doi.org/10.5281/zenodo.13208796).

Acknowledgments

The authors are thankful to the Department of Science and Technology-GoI (grant no. CRG/2020/002080 ), the Department of Biotechnology-GoI (grant no. BT/PR37958/GET/119/297/2020 ), and the Scheme for Promotion of Academic and Research Collaboration (SPARC), MHRD-GoI (grant no. SPARC/2019-2020/P1991/SL ). S.B. thanks the support from the Prime Minister’s Research Fellows (PMRF) Scheme for the fellowship.

Author contributions

S.B., A.K.D., and A.G. conceived the study. S.B. preprocessed the input data and constructed and analyzed the models. S.B. wrote the original draft and generated the figures. S.B. and A.G. wrote the final manuscript. A.G. supervised the study. All authors read and approved the manuscript.

Declaration of interests

The authors declare no competing interests.
==== Refs
References

1 Nanda P. Ghosh A. Genome Scale-Differential Flux Analysis reveals deregulation of lung cell metabolism on SARS-CoV-2 infection PLoS Comput. Biol. 17 2021 e1008860 10.1371/JOURNAL.PCBI.1008860
2 Clemente-Suárez V.J. Martín-Rodríguez A. Redondo-Flórez L. López-Mora C. Yáñez-Sepúlveda R. Tornero-Aguilera J.F. New Insights and Potential Therapeutic Interventions in Metabolic Diseases Int. J. Mol. Sci. 24 2023 10672 10.3390/IJMS241310672 37445852
3 Beura S. Kundu P. Das A.K. Ghosh A. Metagenome-scale community metabolic modelling for understanding the role of gut microbiota in human health Comput. Biol. Med. 149 2022 105997 10.1016/J.COMPBIOMED.2022.105997
4 Beura S. Kundu P. Das A.K. Ghosh A. Sinha P.K. Genome-scale community modelling elucidates the metabolic interaction in Indian type-2 diabetic gut microbiota Sci. Rep. 14 2024 17259 10.1038/s41598-024-63718-0
5 Kundu P. Beura S. Mondal S. Das A.K. Ghosh A. Machine learning for the advancement of genome-scale metabolic modeling Biotechnol. Adv. 74 2024 108400 10.1016/J.BIOTECHADV.2024.108400
6 Chandrasekaran S. Price N.D. Probabilistic integrative modeling of genome-scale metabolic and regulatory networks in Escherichia coli and Mycobacterium tuberculosis Proc. Natl. Acad. Sci. USA 12 2010 17845 17850 10.1073/pnas.1005139107
7 Duarte N.C. Becker S.A. Jamshidi N. Thiele I. Mo M.L. Vo T.D. Srivas R. Palsson B.Ø. Global reconstruction of the human metabolic network based on genomic and bibliomic data Proc. Natl. Acad. Sci. USA 104 2007 1777 1782 17267599
8 Folger O. Jerby L. Frezza C. Gottlieb E. Ruppin E. Shlomi T. Predicting selective drug targets in cancer through metabolic networks Mol. Syst. Biol. 7 2011 501 10.1038/MSB.2011.35 21694718
9 Brunk E. Sahoo S. Zielinski D.C. Altunkaya A. Dräger A. Mih N. Gatto F. Nilsson A. Preciat Gonzalez G.A. Aurich M.K. Recon3D enables a three-dimensional view of gene variation in human metabolism Nat. Biotechnol. 36 2018 272 281 10.1038/nbt.4072 29457794
10 Robinson J.L. Kocabaş P. Wang H. Cholley P.E. Cook D. Nilsson A. Anton M. Ferreira R. Domenzain I. Billa V. An atlas of human metabolism Sci. Signal. 13 2020 eaaz1482 10.1126/SCISIGNAL.AAZ1482
11 Islam M.M. Goertzen A. Singh P.K. Saha R. Exploring the metabolic landscape of pancreatic ductal adenocarcinoma cells using genome-scale metabolic modeling iScience 25 2022 104483 10.1016/J.ISCI.2022.104483
12 Zur H. Ruppin E. Shlomi T. iMAT: an integrative metabolic analysis tool Bioinformatics 26 2010 3140 3142 10.1093/BIOINFORMATICS/BTQ602 21081510
13 Agren R. Bordel S. Mardinoglu A. Pornputtapong N. Nookaew I. Nielsen J. Reconstruction of genome-scale active metabolic networks for 69 human cell types and 16 cancer types using INIT PLoS Comput. Biol. 8 2012 e1002518 10.1371/JOURNAL.PCBI.1002518
14 Agren R. Mardinoglu A. Asplund A. Kampf C. Uhlen M. Nielsen J. Identification of anticancer drugs for hepatocellular carcinoma through personalized genome-scale metabolic modeling Mol. Syst. Biol. 10 2014 721 10.1002/MSB.145122 24646661
15 Heirendt L. Arreckx S. Pfau T. Mendoza S.N. Richelle A. Heinken A. Haraldsdóttir H.S. Wachowiak J. Keating S.M. Vlasov V. Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0 Nat. Protoc. 14 2019 639 702 10.1038/S41596-018-0098-2 30787451
16 Ebrahim A. Lerman J.A. Palsson B.O. Hyduke D.R. COBRApy: COnstraints-Based Reconstruction and Analysis for Python BMC Syst. Biol. 7 2013 74 10.1186/1752-0509-7-74/FIGURES/2 23927696
17 Sheng X. Wu J. Sun Q. Li X. Xian F. Sun M. Fang W. Chen M. Yu J. Xiao J. MTD: a mammalian transcriptomic database to explore gene expression and regulation Brief. Bioinform. 18 2017 28 36 10.1093/BIB/BBV117 26822098
18 Herrmann H.A. Dyson B.C. Vass L. Johnson G.N. Schwartz J.M. Flux sampling is a powerful tool to study metabolism under changing environmental conditions NPJ Syst. Biol. Appl. 5 2019 32 10.1038/s41540-019-0109-0 31482008
19 Bordel S. Agren R. Nielsen J. Sampling the Solution Space in Genome-Scale Metabolic Networks Reveals Transcriptional Regulation in Key Enzymes PLoS Comput. Biol. 6 2010 e1000859 10.1371/JOURNAL.PCBI.1000859
