Daifeng Wang Lab

NetREm is a software package that utilizes network-constrained regularization for biological applications and other network-based learning tasks. In biology, traditional regression methods can struggle with correlated predictors, particularly transcription factors (TFs) that regulate target genes (TGs) in gene regulatory networks (GRNs). NetREm incorporates information from prior biological networks to improve predictions and identify complex relationships among predictors (e.g. TF-TF coordination: direct/indirect interactions among TFs). This approach can highlight important nodes and edges in the network, reveal novel regularized embeddings for genes, provide insights into underlying biological processes, identify subnetworks of predictors that group together to influence the response variable, and improve model accuracy and biological/clinical significance of the models. NetREm can incorporate multiple types of network data, including Protein-Protein Interaction (PPI) networks, gene co-expression networks, and metabolic networks. In summary, network-constrained regularization may bolster the construction of more accurate and interpretable models that incorporate prior knowledge of the network structure among predictors.

Pipeline image of NetREm

The minimum requirement is a computer with 8 GB of RAM and 32 GB of storage. For large prior graph networks, 32 GB of RAM is recommended.

The software uses Python 3.10. After downloading the NetREm Github code, conda/Anaconda users can use the following steps to install:

1. In the Anaconda navigator prompt, create a virtual environment of Python 3.10 by running:
   conda create -n NetREm python=3.10
2. Activate the environment:
   conda activate NetREm
3. Make sure to change the current directory to the NetREm folder.
4. Install the packages and dependencies (math, matplotlib, networkx, numpy, typing, os, pandas, plotly.express, random, scipy, scikit-optimize, sklearn, sys, tqdm, warnings):
   pip install -r requirements.txt

Please note that if you encounter import errors from files or functions in the code folder (such as Netrem_model_builder.py), add an empty file named init.py to the code folder, and add the "code." prefix to all imports from the "code" folder. For example, import Netrem_model_builder as nm ➡️ import code.Netrem_model_builder as nm.

NetREm fits a Network-constrained Lasso regression machine learning model with user-provided weights for the prior network. Here, netrem is the main function with the following usage:

netrem(
        edge_list,
        beta_net = 1,
        alpha_lasso = 0.01,
        default_edge_weight = 0.01,
        edge_vals_for_d = True,
        w_transform_for_d = "none",
        degree_threshold = 0.5,
        gene_expression_nodes = [],
        overlapped_nodes_only = False,
        y_intercept = False,
        view_network = False,
        model_type = "Lasso",
        ...
)

We input an edge list of the prior graph network (constrains the model via network-based regularization) and a beta_net ($\beta_{net} \geq 0$, which scales the network-based regularization penalty). The user may specify the alpha_lasso ($\alpha_{lasso} \geq 0$) manually for the lasso regularization on the overall model (if model_type = Lasso) or NetREm may select an optimal $\alpha_{lasso}$ based on cross-validation (CV) on the training data (if model_type = LasssoCV). Then, netrem builds an estimator object from the class Netrem that can then take in input $X$ and $y$ data: transforms them to $\tilde{X}$ and $\tilde{y}$, respectively, and use them to fit a Lasso regression model with a regularization value of $\alpha_{lasso}$. Ultimately, the trained NetREm machine learning model is more reflective of an underlying network structure among predictors and may be more biologically meaningful and interpretable. Nonetheless, NetREm could be applied in various contexts where a network structure is present among the predictors. Input networks are typically weighted and undirected. We provide details, functions, and help with converting directed networks to undirected networks (of similarity values among nodes) here.

• A NetREm network-regularized linear model estimator object from the NetREmModel class.

We assume that our $X$ and $y$ data correspond to $M$ samples and $N$ predictors. For biological applications, the $X$ data would typically be gene expression data (bulk or single-cell) for the $N$ candidate predictors (usually Transcription Factors (TFs)) for the $M$ samples. Then, the $y$ values would correspond to the gene expression values for the target gene (TG) $y$ for those same $M$ samples.

We can retrieve our model coefficients and other attributes by calling these outputs:

• predict($X$)

We can use our model to predict values for our response variable $y$.

Returns: Numpy array of $\hat{y}$ predicted values for $y$.

• test_mse($X$, $y$)

We can evaluate our model performance capabilities on data like testing data using the Mean Squared Error (MSE) as our metric.

Returns: Numeric value corresponding to the Mean Square Error (MSE).

$$ MSE = \frac{1}{m} \sum_{i=1}^m (y_i - \hat{y_i})^2 $$

=======

The goal is to build a machine learning model to predict the gene expression levels of our target gene (TG) $y$ based on the gene expression levels of $N = 5$ Transcription Factors (TFs) [TF₁, $TF_{2}$, $TF_{3}$, $TF_{4}$, $TF_{5}$] in a particular cell-type. Assume the gene expression values for each TF are [X₁, $X_{2}$, $X_{3}$, $X_{4}$, $X_{5}$], respectively. We generate $10,000$ random cells (rows: samples) of data where the Pearson correlations ($r$) between the gene expression of each TF ($X$) with gene expression of TG $y$ as corrVals: [cor(TF₁, $y$) = 0.9, cor(TF₂, $y$) = 0.5, cor(TF₃, $y$) = 0.4, cor(TF₄, $y$) = -0.3, cor(TF₅, $y$) = -0.8]. We note that the sparsity_factor_perc is 40%, so roughly 4,000 out of the 10,000 cells for each variable will be ≈0; we incorporate these dropoout, so that the underlying data mimics single-cell expression data. Thus, $M$ is $7,000$ cells since that is what is used for training the NetREm model.

The dimensions of $X$ are therefore 10,000 rows by 5 columns (predictors). More details about our generate_dummy_data function (and additional parameters we can adjust for) are in Dummy_Data_Demo_Example.pdf. Our NetREm estimator also incorporates a constraint of an undirected weighted prior graph network of biological relationships among only 5 TFs based on a weighted Protein-Protein Interaction (PPI) network ([TF₁, $TF_{2}$, $TF_{3}$, $TF_{4}$, $TF_{5}$]), where higher edge weights $w$ indicate stronger biological direct and/or indirect pairwise interactions among TFs at the protein-level. These linked TFs may be involved in several shared functions at the molecular level, such as coordinating to co-regulate shared target genes (TGs).

Please note that the code for this demo example is demo_toy.py in the demo folder.

from DemoDataBuilderXandY import generate_dummy_data
from Netrem_model_builder import netrem
import PriorGraphNetwork as graph
import error_metrics as em 
import essential_functions as ef
import netrem_evaluation_functions as nm_eval

num_samples = 10000 # number of samples (e.g. single-cells)
corrs_list =  [0.9, 0.5, 0.4, -0.3, -0.8] # the individual correlations (r) of each X variable with y
train_p = 70 # the % of original data we use for M
dummy_data = generate_dummy_data(corrVals = corrs_list, # the # of elements in corrVals is the # of predictors (X)
                                 num_samples_M = num_samples, # the number of samples M is: (num_samples_M)*(train_data_percent)
                                 sparsity_factor_perc = 40, # approximately what % of data for each variable is sparse (0)
                                 train_data_percent = train_p) # the remainder out of 10,000 will be kept for testing. If 100, then ALL data is used for training and testing.

The Python console or Jupyter notebook will print out the following:

same_train_test_data = False
Please note that since we hold out 30.0% of our 100000 samples for testing, we have:
X_train = 7,000 rows (samples) and 5 columns (N = 5 predictors) for training.
X_test = 3,000 rows (samples) and 5 columns (N = 5 predictors) for testing.
y_train = 7,000 corresponding rows (samples) for training.
y_test = 3,000 corresponding rows (samples) for testing.

The $X$ data should be in the form of a Pandas dataframe as below:

X_df = dummy_data.X_df
X_df.head(10) # please view the first 10 rows

y_df = dummy_data.y_df
y_df.head(10) # please view the first 10 rows

# M = 7,000 samples (cells) for training data (used to train and fit NetREm model)
X_train = dummy_data.view_X_train_df()
y_train = dummy_data.view_y_train_df()

# 3,000 samples (cells) for testing data
X_test = dummy_data.view_X_test_df()
y_test = dummy_data.view_y_test_df()

Our generated single-cell gene expression data (X, y) looks like this:

X_train.corr() # pairwise correlations among the training samples

X_test.corr() # pairwise correlations among the testing samples

dummy_data.combined_correlations_df # breakdown of the correlations in training, testing, and overall data

We note that NetREm will make this input TF-TF Protein-Protein Interaction (PPI) Network (PPIN) fully connected by adding a default edge weight ŋ = 0.01 for missing edges. This is for numerical stability and to also help propel the discovery of novel TF-TF coordination links.

# prior network edge_list (missing edges or edges with no edge weight will be added with the default_edge_list so the network is fully-connected):
edge_list = [["TF1", "TF2", 0.8], ["TF4", "TF5", 0.95], ["TF1", "TF3"], ["TF1", "TF4"], ["TF1", "TF5"], 
             ["TF2", "TF3"], ["TF2", "TF4"], ["TF2", "TF5"], ["TF3", "TF4"], ["TF3", "TF5"]]

beta_network_val = 1 
alpha_lasso_val = 0.1 # the default, so we do not need to specify it. 
# by default, model_type is Lasso, so alpha_lasso_val will be specified for the alpha_lasso parameter. 
# However, if we specify model_type = Lasso, so our alpha_lasso is determined by cross-validation on training data).

# Building the network regularized regression model: 
# By default, edges are constructed between all of the nodes; nodes with a missing edge are assigned the default_edge_weight. 
netrem_demo = netrem(edge_list = edge_list, 
                     beta_net = beta_network_val,
                     alpha_lasso = alpha_lasso_val, # the default is 0.1. 
                     model_type = "Lasso",
                     view_network = True)

# Fitting the NetREm model on training data: X_train and y_train:
netrem_demo.fit(X_train, y_train)

To view and extract the predicted model coefficients c for the predictors:

netrem_demo.model_coef_df

In the context of gene regulation (in biology), we predict that predictors with negative NetREm coefficients for target gene (TG) $y$ may be repressors (their activity focuses on reducing expression of $y$) and those with positive coefficients for $y$ may be activators.

To view the TG-specific TF-TF coordination (of direct and/or indirect interactions among TFs) that NetREm learned for this target gene $y$, in our given cell-type, we can view the B_interaction_df.

netrem_demo.B_interaction_df

Based on the NetREm model, we predict:

• TF predictors with positive $c^{*} > 0$ coefficient signs are activators of TG: linked to increased gene expression levels $y$ of the TG.
• TF predictors with negative $c^{*} < 0$ coefficient signs are repressors of TG: linked to decreased gene expression levels $y$ of the TG.
• TF predictors with negative $c^{*} = 0$ coefficient signs are not selected as final TFs for the TG based on the input data.

Our TG-specific TF-TF coordination network (i.e. coordination network among predictors for this $y$ variable) has coordination scores $B$ that are between -100 and 100. That is, -100 ≤ $B$ ≤ 100.

These tend to reflect the grouped variable selection property of network-regularized regression models that tend to select TF predictors that also group together in biological pathways in the original Protein-Protein Interaction Network (PPIN).

We predict positive scores 0 < $B$ ≤ 100 between predictors with the same coefficient sign (e.g. $TF_1$ and $TF_2$ have a score of 71.466718, which is higher than their links with $TF3$ since TF3 has weak links with them in the input TF-TF PPIN). This implies potential cooperativity (e.g. cobinding, pioneer/settler models) among them. They have a shared goal and $B &gt; 0$:

• TF1, TF2, TF3 have a common goal of co-activating TG expression levels $y$ ➡️ increased $y$
• TF4 and TF5 have a common goal of co-repressing TG expression levels $y$ ➡️ decreased $y$

We predict negative scores -100 ≤ $B$ < 0 between predictors with different coefficients sign (e.g. $TF_1$ and $TF_5$ have a score of -97.923796). This implies potential antagonistic activity (e.g. competition for binding sites, sequestration) among them: activator-repressor antagonism. For instance, these TFs may compete to regulate the $TG$ through biological mechanisms that may be investigated further through experiments. They have opposing goals that may partly cancel each other out so $B &lt; 0$:

• TF1, TF2, TF3 each have antagonistic relations with TF4 and TF5.

We can test the performance of our data on testing data $X_{test}$ (3,000 samples), to understand better the generalizability of our NetREm model on new, unseen, data.

pred_y_test = netrem_demo.predict(X_test) # predicted values for y_test
mse_test = netrem_demo.test_mse(X_test, y_test)

print(f"The testing Mean Square Error (MSE) is {mse_test}")

The testing Mean Square Error (MSE) is 0.21834495195611514

We can analyze more metrics about our NetREm model results as below, where:

• low NMSE values are great
• high values for SNR and PSNR are ideal:

nmse_test  = em.nmse(y_test.values.flatten(), pred_y_test)
print(f"The testing Normalized Mean Square Error (NMSE) is {nmse_test}")

The testing Normalized Mean Square Error (NMSE) is 0.22013278683483584

snr_test  = em.snr(y_test.values.flatten(), pred_y_test)
print(f"The testing Signal-to-Noise Ratio (SNR) is {snr_test}")

The testing Signal-to-Noise Ratio (SNR) is 6.573152683008163

psnr_test  = em.psnr(y_test.values.flatten(), pred_y_test)
print(f"The testing Peak Signal-to-Noise Ratio (PSNR) is {psnr_test}")

The testing Peak Signal-to-Noise Ratio (PSNR) is 17.629710652383

These are additional metrics:

netrem_demo.final_corr_vs_coef_df

In the context of gene regulation, our results may thereby be interpreted in the following way:

Nonetheless, NetREm can be applied to solve a suite of regression problems where there is an underlying connection among the predictors (that group together in meaningful subnetworks to influence the outcome $y$) and their correlation with one another may be utilized jointly for the predictive task rather than discarded.

We also provide a suite of evaluation functions and explanations of more advanced functionalities in our User Guide. For instance, we provide more details on our netrem and netremcv functions in terms of predicting gene regulation in myelinating Schwann cells here.

netrem_demo.combined_df

netrem_demo.TF_coord_scores_pairwise_df

organize_predictor_interaction_network(netrem_demo)

We also note that since w_transform_for_d = "none", then to calculate the degree of each TF predictor node, we calculate the sum of all edges $w_{ij}$ for $j = 1$ to $N$ where $j \neq i$. Then, the degree $d$ of $TF_i$ is denoted as $d_i$, and it can be expressed as:

$$ d_i = \sum_{\substack{j=1 \ j \neq i}}^N w_{ij} $$

Instead of self-loops, we have:

$$ w_{ii} = \frac{d_i}{N-1} $$

We preprocess the input TF-TF PPI network to make it fully-connected for all pairwise TF-TF edges (default_edge_weight of 0.01) and diagonals reflecting averaged connectivity values.

Please note that the manuscript associated with NetREm is available as a pre-print on biorxiv: Link to Paper

Citation:

• Saniya Khullar, Xiang Huang, Raghu Ramesh, John Svaren, Daifeng Wang, NetREm: Network Regression Embeddings reveal cell-type transcription factor coordination for gene regulation, bioRxiv 2023.10.25.563769; doi: https://doi.org/10.1101/2023.10.25.563769

[1]: Caiyan Li, Hongzhe Li, Network-constrained regularization and variable selection for analysis of genomic data, Bioinformatics, Volume 24, Issue 9, May 2008, Pages 1175–1182, https://doi.org/10.1093/bioinformatics/btn081

[2]: Wilkinson, M.D., Dumontier, M., Aalbersberg, I.J. et al. The FAIR Guiding Principles for scientific data management and stewardship. Sci Data 3, 160018 (2016). https://doi.org/10.1038/sdata.2016.18

[3]: Zhang, Y., Akutsu, T., & Ching, W. K. (2016). Incorporating network topology and phylogenetic constraints for protein-protein interaction prediction. BMC bioinformatics, 17(1), 1-15. https://doi.org/10.1186/s12859-016-1310-4

[4]: Jia, C., Zhang, Y., Chen, K., & Zhang, S. (2018). A novel feature extraction method with improved consistency for identifying cell cycle regulated genes. Bioinformatics, 34(5), 896-903. https://doi.org/10.1093/bioinformatics/btx657

[5]: Lu, Y., Chen, X., & Hu, Z. (2017). Recognition of protein/gene names from text using an ensemble of classifiers and effective abbreviation resolution. BMC bioinformatics, 18(1), 1-11. https://doi.org/10.1186/s12859-017-1515-1

[6]: Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., ... & Vanderplas, J. (2011). Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12(Oct), 2825-2830.