
==== Front
Comput Struct Biotechnol J
Comput Struct Biotechnol J
Computational and Structural Biotechnology Journal
2001-0370
Research Network of Computational and Structural Biotechnology

S2001-0370(24)00270-8
10.1016/j.csbj.2024.08.010
Research Article
NJGCG: A node-based joint Gaussian copula graphical model for gene networks inference across multiple states
Huang Yun ab1
Huang Sen c1
Zhang Xiao-Fei d
Ou-Yang Le leouyang@szu.edu.cn
c⁎
Liu Chen liuchen1985@fjmu.edu.cn
efg⁎⁎
a Department of Geriatrics, The First Affiliated Hospital of Fujian Medical University, Fuzhou 350005, China
b Clinical Research Center for Geriatric Hypertension Disease of Fujian province, The First Affiliated Hospital of Fujian Medical University, Fuzhou 350005, China
c Guangdong Key Laboratory of Intelligent Information Processing, College of Electronics and Information Engineering, Shenzhen University, Shenzhen, China
d School of Mathematics and Statistics & Hubei Key Laboratory of Mathematical Sciences, Central China Normal University, Wuhan, China
e Department of Oncology, Molecular Oncology Research Institute, The First Affiliated Hospital of Fujian Medical University, Fuzhou 350005, China
f Department of Oncology, National Regional Medical Center, Binhai Campus of The First Affiliated Hospital, Fujian Medical University, Fuzhou 350212, China
g Fujian Key Laboratory of Precision Medicine for Cancer, The First Affiliated Hospital of Fujian Medical University, Fuzhou 350005, China
⁎ Corresponding author. leouyang@szu.edu.cn
⁎⁎ Corresponding author at: Department of Oncology, Molecular Oncology Research Institute, The First Affiliated Hospital of Fujian Medical University, Fuzhou 350005, China. liuchen1985@fjmu.edu.cn
1 These authors contributed equally to this work.

22 8 2024
12 2024
22 8 2024
23 31993210
14 4 2024
5 8 2024
11 8 2024
© 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/).
Inferring the interactions between genes is essential for understanding the mechanisms underlying biological processes. Gene networks will change along with the change of environment and state. The accumulation of gene expression data from multiple states makes it possible to estimate the gene networks in various states based on computational methods. However, most existing gene network inference methods focus on estimating a gene network from a single state, ignoring the similarities between networks in different but related states. Moreover, in addition to individual edges, similarities and differences between different networks may also be driven by hub genes. But existing network inference methods rarely consider hub genes, which affects the accuracy of network estimation. In this paper, we propose a novel node-based joint Gaussian copula graphical (NJGCG) model to infer multiple gene networks from gene expression data containing heterogeneous samples jointly. Our model can handle various gene expression data with missing values. Furthermore, a tree-structured group lasso penalty is designed to identify the common and specific hub genes in different gene networks. Simulation studies show that our proposed method outperforms other compared methods in all cases. We also apply NJGCG to infer the gene networks for different stages of differentiation in mouse embryonic stem cells and different subtypes of breast cancer, and explore changes in gene networks across different stages of differentiation or different subtypes of breast cancer. The common and specific hub genes in the estimated gene networks are closely related to stem cell differentiation processes and heterogeneity within breast cancers.

Graphical abstract

Keywords

Gene network
Graphical model
Gene expression
==== Body
pmc1 Introduction

Genes rarely function in isolation. They typically perform specific biological roles through interactions with other genes or biological macromolecules. The intricate regulatory relationships between genes form a complex gene network, serving as the mechanism through which various biological functions are executed [1]. Consequently, constructing gene networks is crucial for understanding the functional organization of living organisms and elucidating the pathogenesis of complex diseases [2], [3], [4].

Traditional methods rely on biological experiments to identify the regulatory relationships between genes. However, these methods are often low-throughput and time-consuming, making it challenging to dynamically and systematically generate large-scale gene networks. Advances in high-throughput sequencing technology have enabled the acquisition of extensive gene expression data, creating new opportunities to use computational methods for inferring gene networks. In recent years, a series of computational network inference methods have been proposed to leverage gene expression data for the inference of gene networks [5], [6], [7], [8], [9], [10], [11], [12], [13].

Existing network inference methods can be broadly categorized into three types: correlation coefficient-based methods [14], [15], mutual information-based methods [6], [8], [16], and Gaussian graphical model-based methods [17]. Correlation coefficient-based methods calculate the correlation coefficient to infer the co-expression relationships between genes, but the correlations inferred by the correlation coefficient encompass both direct and indirect correlations. Mutual information-based methods are capable of inferring the nonlinear associations between genes [6], [8], [16], but they do not distinguish between direct and indirect associations. Methods based on Gaussian graphical models can infer the direct associations between genes by examining their conditional dependencies [18]. However, the sample sizes for gene expression data are often small relative to the number of genes being studied. Moreover, research indicates that real gene networks are usually sparse [19], [17], [20], [21], [22], which necessitates the use of sparse estimation methods for accurate network inference. To address these challenges, various sparse Gaussian graphical models have been proposed [19], [17], [23], [24]. While Gaussian graphical models depend on the assumption of multivariate Gaussianity, real gene expression data may not follow a Gaussian distribution. To address this issue, some semiparametric Gaussian copula models have been introduced to relax the Gaussian assumption [25], [26], [27].

Traditional gene network inference methods have primarily focused on inferring a single network [19], [17], [25], [26], [27]. However, gene networks are dynamic and can change with varying environmental and cellular states. For instance, during the differentiation of mouse embryonic stem cells [28], similar network structures and key differences (see Fig. 1) may exist between the gene networks at different differentiation stages [29]. Combining samples from multiple states to estimate a single gene network can overlook these state-specific differences. Conversely, estimating gene networks separately for each state may ignore the similarities between states, which can impact the accuracy of network inference. Most existing methods for network inference struggle with multi-network estimation. However, methods based on Gaussian graphical models can more naturally extend from single network inference to multi-network inference. Recently, several multi-network joint inference methods based on Gaussian graphical models have been proposed [30], [31], [32], [33], [34], [35], [36], [37], [38].Fig. 1 Motivation and illustration of our model. There are some central genes involved in the processes of cell differentiation. During the four stages of differentiation, gene 3 serves as a common hub gene shared across all stages, while genes 1, 5 and 6 are three specific hub genes unique to certain differentiation stages.

Fig. 1

While existing multi-network joint inference methods, such as JEGN [39], are capable of capturing shared and specific subnetworks between different networks and handling non-Gaussian data, they still have several limitations. Firstly, these methods often overlook the role of hub genes within the networks. Real gene networks typically include highly connected genes that interact with a large number of other genes [40]. The hub genes in the gene networks of different states can either be common hub genes shared across all gene networks (e.g., Gene 3 in Fig. 1), or specific hub genes present only in certain gene networks (e.g., Genes 1 and 6 in Fig. 1). Identifying both common and specific hub genes when estimating multiple gene networks can provide a clearer understanding of how gene networks change across different states. Secondly, existing methods often neglect the issue of missing values in gene expression data. Real gene expression data frequently contains missing values due to technical limitations, indicating that some truly expressed genes may not be detected [41]. For instance, single-cell RNA sequencing (scRNA-seq) data often contains an excessive number of zero counts, particularly for genes with low or moderate expression, due to the presence of dropout events. These missing values can impede accurate gene network inference.

In order to address the aforementioned challenges, this study introduces a novel node-based joint Gaussian copula graphical (NJGCG) model to jointly estimate multiple gene networks based on gene expression data collected from various states. The NJGCG model does not depend on a specific data distribution and is adept at handling missing values, making it suitable for a wide range of gene expression data types. Additionally, by incorporating a tree-structured group Lasso penalty (Fig. 2), our model can identify both common and specific hub genes across different gene networks. To solve the optimization problem, an alternating direction method of multipliers (ADMM) [42], [43], [44] is employed. Simulation studies demonstrate that out NJGCG model outperforms other methods in all tested scenarios. Furthermore, we apply our model to two real datasets to infer the gene networks for different stages of cell differentiation and various breast cancer subtypes. The common and specific hub genes identified in our estimated gene networks have been found to be closely associated with the differentiation processes of stem cells and the heterogeneity within breast cancers.Fig. 2 Tree-structured group Lasso penalty. Functions in different levels are used to capture different aspects of gene networks: The first level produces a sparse network; The second level identifies specific hub genes that are unique to certain networks; The third level captures common hub genes that shared across all states.

Fig. 2

2 Methods

2.1 Non-paranormal distribution

A p-dimensional random variable X=(X1,…,Xp)T follows a non-paranormal distribution X∼NPNp(f,Σ) if there is a set of univariate monotone functions f1,…,fp and f(X)=(f1(X1),…,fp(Xp)) follows a multi-variate normal distribution f(X)∼Np(0,Σ). In this context, the conditional dependency relationships between the p random variables in X can be determined via the inverse of covariance matrix (or precision matrix) Θ=Σ−1, where Θij≠0 indicates random variables Xi and Xj are conditionally dependent given all other variables [25], [26].

2.2 Problem statement and notations

Suppose we have collected the gene expression data of K sets of samples in different states (e.g., different states can be different stages of differentiation, or different cancer subtypes, etc.), which characterize K different gene networks. These gene expression data include a total of n=n1+n2+⋯+nk samples that are independently observed for p common genes. For k-th (k∈{1,2,…,K}) state, there are nk samples: x1(k),x2(k),…,xnk(k)∈Rp. Assuming that the samples under each state are independently and identically distributed and follow a p-dimensional non-paranormal distribution: x1(k),x2(k),…,xnk(k)∼NPNp(f(k),Σ(k)), where Σ(k)∈Rp×p denotes the corresponding covariance matrix under k-th state. Based on this assumption, the conditional dependency relationships between p genes can be determined via the precision matrix Θ(k)=Σ(k)−1. Our goal is to infer K precision matrices (which correspond to K gene networks) from gene expression data collected from K different states.

Before introducing our model, we briefly describe some symbols that will be used. For example, trace(⋅) denotes the trace of a matrix, det(⋅) denotes the determinant of a matrix. ‖A‖1=∑i,j=1p|aij| denotes the ℓ1-norm of matrix A. ‖A‖F=∑i,j=1paij2 denotes the Frobenius norm of matrix A, and 〈A,B〉=trace(ABT) denotes the inner product of two matrices A and B with same size.

2.3 Model formulation

Let X(k)=[x1(k);…;xnk(k)]∈Rnk×p denote the gene expression matrix of k-th state, where x1(k),…,xnk(k)∈Rp are nk samples. Since there may be missing values in the gene expression matrix X(k), let bij(k)=1 if Xij(k) is observed, and bij(k)=0 otherwise.

To utilize the similarities between samples under different states to improve the accuracy of the estimated networks Θ={Θ(1),…,Θ(K)}, we first introduce the following loss function to jointly estimate multiple networks [17]:(1) L(X,Θ)=∑k=1Knk[−log⁡(det⁡(Θ(k)))+trace(S(k)Θ(k))],

where S(k) denotes the empirical sample covariance matrix of k-th state. In order to estimate S(k) from the samples of k-th state (i.e., X(k)), which follow a non-paranormal distribution, we adopt a rank-based estimator named Kendall's tau estimator [25], [26], [39]. Note that there are missing values in the gene expression data. To handle samples with missing values, following previous studies [45], [46], we introduce the following improved Kendall's tau estimator to calculate the empirical estimation of S(k).

We compute the Kendall's tau correlation between genes y and z using the effective independent samples myz(k)=∑i=1nkbiy(k)biz(k) that have values for both genes. Specifically, the Kendall's tau correlation between genes y and z is computed as follows:(2) τˆyz(k)=m¯yz(k)∑i≠i′nkτii′yz(k),

where τii′yz(k)=biy(k)biz(k)bi′y(k)bi′z(k)sign((xiy(k)−xi′y(k))(xiz(k)−xi′z(k))) and m¯yz(k)=1myz(k)(myz(k)−1). Following [47], [48], [26], we adopt the following definition to estimate S(k):(3) Syz(k)={sin⁡(π2τˆyz(k))ify≠z,1ify=z.

To make sure that the sample covariance matrices S(k) are positive semidefinite, following [39], we project the estimated matrices into the cone of positive semidefinite matrices by solving the following optimization problem:(4) Sˆ(k)=arg⁡minR⪰0⁡‖S(k)−R‖∞.

Then replacing the value of S(k) with Sˆ(k).

In practical scenarios, the number of edges in a gene network is typically much smaller than that in a fully connected network, indicating that the gene network should exhibit sparsity [30], [33], [39]. Additionally, prior research has demonstrated the existence of key driver genes within gene networks [49], [50], which regulate a large number of downstream genes. Consequently, there are densely connected genes (referred to as hub genes) in the gene network [40]. Moreover, gene networks under different conditions may share common network structures, implying the potential presence of shared hub genes. Therefore, we assume that gene networks across different conditions contain both common hub genes shared among all networks and specific hub genes unique to individual networks. To simultaneously control the sparseness of the estimated gene networks and capture all types of hub genes (including common and specific hub genes), we define the precision matrix Θ(k) as Θ(k)=Z(k)+V(k)+V(k)T, where Z(k) is used to capture the sparse interactions between genes, and V(k) is used to capture the common hub nodes shared across different states and the specific hub nodes which are unique in certain states.

To make use of the similarities between networks under different states and make the estimated networks sparse, we combine the sparse penalty and tree-structured group Lasso penalty [51], [52] to get the following penalty function:(5) P({Θ(k)})=λ1∑i,j,k|Zij(k)|+λ2[ω1∑i,j,k|Vij(k)|+ω2∑j,k∑i(Vij(k))2+ω3∑j∑i,k(Vij(k))2],Θ(k)=Z(k)+V(k)+V(k)T,k=1,…,K,

where Z(k) is a symmetric matrix and V(k)T is the transpose of V(k) for k=1,…,K. Here V(k) does not need to be symmetrical as V(k)+V(k)T is symmetrical. Based on the above penalty function, we can identify hub genes by group penalizing the columns of V(k). In order to identify both common and specific hub genes, we apply tree-structured group Lasso penalty to {V(k)}k=1,…,K. For j-th column of V(k), a total of Kp parameters are divided into three levels according to the tree-structure group Lasso penalty (Fig. 2).

The first term of the tree-structured group Lasso penalty imposes sparse penalty on V(k) and encourages the sparsity within V(k). The second term imposes group sparsity penalty on columns of V(k) to identify hub nodes that are unique to certain networks. The third term imposes group sparsity penalty on columns of V(k), for k=1,…,K to identify hub nodes that are common to all networks. Here, different parameters λ1, λ2, ω1, ω2 and ω3 are assigned to different penalty functions, where λ1, λ2 and ω1 control the sparsity of the estimate networks, λ2, ω2 and ω3 control the identification of common and specific hub nodes.

By combining the loss function in Eq. (1) and the penalty function in Eq. (5), we obtain the following node-based joint Gaussian copula graphical model (NJGCG) to jointly estimate multiple gene networks which have common and specific hub genes:(6) minΘ,Z,V∑k=1Knk[−log⁡(det⁡(Θ(k)))+trace(S(k)Θ(k))]+λ1∑i,j,k|Zij(k)|+λ2[ω1∑i,j,k|Vij(k)|+ω2∑j,k∑i(Vij(k))2+ω3∑j∑i,k(Vij(k))2],s.t.Θ(k)=Z(k)+V(k)+V(k)T,

where λ1 and λ2 are non-negative tuning parameters controlling the sparsity of the estimated networks and the identification of hub genes, respectively. We will introduce our parameter selection strategy in the last subsection of this section.

2.4 Algorithm

Inspired by [40], we use alternate direction method of multipliers (ADMM) algorithm to solve Eq. (6) [42], [43], [44]. In particular, we rewrite Eq. (6) into Eq. (7) to ensure a convergent solution.(7) f(B)=L(X,Θ)+P({Θ(k)}),

where B=(Θ,V,Z), and Θ=(Θ(1),Θ(2),…,Θ(K))T, V=(V(1),V(2),…,V(K))T, Z=(Z(1),Z(2),…,Z(K))T, K is the number of networks. Then we define following functions:(8) g(B˜)={0ifΘ˜(k)=Z˜(k)+V˜(k)+V˜(k)T,∞otherwise,

where B˜=(Θ˜,V˜,Z˜), and Θ˜=(Θ˜(1),Θ˜(2),…,Θ˜(K))T, V˜=(V˜(1),V˜(2),…,V˜(K))T, Z˜=(Z˜(1),Z˜(2),…,Z˜(K))T. Then the optimization problems become Eq. (9), and the algorithm convergence of this optimization problem still follows the classical solution [43], [44].(9) minB,B˜{f(B)+g(B˜)},s.t.B=B˜.

The augmented Lagrangian term of Eq. (9) is as follows:(10) L(B,B˜,W)=L(X,Θ)+P({Θ(k)})+g(B˜)+ρ2‖B−B˜+W‖F2.

where B and B˜ is primal variables, and W=(WΘ,WZ,WV) is dual variable [43].

Then we introduce the following steps to update parameters:1. B[t+1]←argminBL(B,B˜[t],W[t]),

2. B˜[t+1]←argminB˜L(B[t+1],B˜,W[t]),

3. W[t+1]←W[t]+B[t+1]−B˜[t+1].

Note that the variables: Θ,V,Z can be calculated separately. Θ can be updated base on the loss function L(X,Θ), and V and Z can be updated according to the rules given in Algorithm 1.

Minimizing B˜ is equal to:(11) min⁡{ρ2‖Θ−Θ˜+WΘ‖F2+ρ2‖V−V˜+WV‖F2+ρ2‖Z−Z˜+WZ‖F2},s.t.Θ˜=Z˜+V˜+V˜T.

Let Γ(k) be the p×p Lagrange multiplier matrix for the equality constraint. Then the Lagrangian for Eq. (11) is:L(Θ˜(k),Z˜(k),V˜(k),Γ(k))=ρ2‖Θ˜(k)−(Θ(k)+WΘ(k))‖F2+ρ2‖Z˜(k)−(Z(k)+WZ(k))‖F2+ρ2‖V˜(k)−(V(k)+WV(k))‖F2+tr[Γ(k)T(Θ˜(k)−Z˜(k)−V˜(k)−V˜(k)T)],

and(12) {Θ˜(k)=Θ(k)+WΘ(k)−1ρΓ(k),Z˜(k)=1ρΓ(k)+Z(k)+WZ(k),V˜(k)=1ρ(Γ(k)+Γ(k)T)+V(k)+WV(k),

where(13) Γ(k)=ρ6[(Θ(k)+WΘ(k))−(Z(k)+WZ(k))

(14) −(V(k)+WV(k))−(V(k)+WV(k))T].

Algorithm 1 ADMM Algorithm for Solving NJGCG.

Algorithm 1

It has been proved that the optimization problem in step (a) (iii) admits an analytical solution [53]. In this study, following [54], [55] we use the algorithm contained in the Sparse Learning with Efficient Projections (SLEP) package to solve this problem.

2.5 Selection of tuning parameters

For NJGCG, there are five tuning parameters: ω1, ω2, ω3, λ1 and λ2 that need to be predefined. In this study, to simplify our model and facilitate the use of our model, similar to previous study [55], we set the empirical values of ω1, ω2 and ω3 as (ω1,ω2,ω3)=(0.01,ln⁡p,ln⁡(p(K−1))). The tuning parameters λ1 and λ2 control the network sparsity and the identification of hub nodes, respectively. Inspired by [56], [57], [58], we choose the tuning parameters for NJGCG by using an approximation of the Bayesian Information Criterion (BIC).(16) BIC(Θˆ,Vˆ,Zˆ)=−∑k=1Knk[log⁡(det⁡(Θˆ(k)))−trace(S(k)Θˆ(k))]+∑k=1Klog⁡(nk)⋅[|Zˆ(k)|+ν(k)+c⋅(|Vˆ(k)|−ν(k))],

where ν(k) is the number of estimated hub nodes of k-th network, i.e., ν(k)=∑j=1p1{‖Vˆj‖0>0}, c is a constant between zero and one, and |Zˆ(k)| and |Vˆ(k)| are the cardinalities (the number of unique non-zeros) of Zˆ(k) and Vˆ(k), respectively. We select the values of tuning parameters (λ1,λ2) that minimize the quantity BIC(Θˆ,Vˆ,Zˆ). Note that when the constant c is small, BIC(Θˆ,Vˆ,Zˆ) will favor more hub nodes in Vˆ. In this manuscript, we set c=0.2.

3 Simulation studies

3.1 Data generation

3.1.1 Gaussian data

In simulation studies, we consider both scale-free network and Erdos–Renyi network [59], [60]. We generate K scale-free networks (or K Erdos–Renyi networks), each of them contains p=100 common nodes. In each network, the number of common hub nodes is set to ncomm=2, and the number of specific hub nodes is set to nspec=4. The sample size and the number of conditions are set to nk=200,400,600 and K=4, respectively. The detailed processes of generating simulated data are as follows:Step 1: Generating a scale-free (or a Erdos–Renyi) network with p nodes, and A∈{0,1}p×p denotes the adjacency matrix, where Aij=1 indicate that there is an edge between nodes i and j, and Aij=0 otherwise.

Step 2: For each edge in A, a uniform distribution Unif([−1,−0.5]∪[0.5,1]) is used to generate its weight, i.e.,Aij∼i.i.d.{0ifAij=0,Unif([−1,−0.5]∪[0.5,1])ifAij=1.

Step 3: Randomly selecting nhubs=ncomm+nspec nodes in the network as hub nodes, where ncomm represents the number of common hubs shared by all networks, and nspec represents the number of specific hubs unique to individual networks.

Step 4: For hub nodes, the sparsity of the elements in the corresponding rows and columns of matrix A is reset to 0.3, with the non-zero elements drawn from a uniform distribution Unif([−1,−0.5]∪[0.5,1]).

Step 5: To make the adjacency matrix symmetric, we set A¯=12×(A+AT). To ensure that the adjacency matrix is positive definite, we adjust it to Θ(k)=A¯+(0.1−λmin(A¯))×I, where λmin(A¯) denotes the smallest eigenvalue of A¯, and I is a p×p identity matrix.

Step 6: Generating nk independent samples x1(k),…,xnk(k) from a normal distribution N(0,(Θ(k))−1). Let X(k)=[x1(k);…;xnk(k)]∈Rnk×p denote the sample matrix.

Step 7: Repeating Steps 2-6 to generate K networks and their corresponding sample matrices.

Step 8: To simulate the missing data scenarios, the entries in X(k) are set to zero with a probability δ. Then the sample matrices are used to calculate the sample covariance matrices S(k).

3.1.2 Non-Gaussian data

Similar to [61], we utilize the data generation procedure outlined in [62] to create non-Gaussian data. The networks (scale-free or Erdos–Renyi networks) are generated similarly to those for Gaussian data (following Steps 1–3). Then we employ the ‘‘SeqNet’’  package to generate non-Gaussian gene expression data, utilizing the default settings for the breast cancer dataset.

3.2 Simulation results

We evaluate the performance of NJGCG by comparing it with three state-of-the-art network inference methods, i.e., joint graphical Lasso (GGL) [30], the co-hub node joint graphical lasso (CNJGL) [49] and JEGN [39].

The GGL model is designed for the joint estimation of multiple networks by employing the group Lasso penalty to regulate similarities between individual networks. In contrast, to account for the presence of hub nodes in the network, the CNJGL model incorporates the hub structure of networks to facilitate the joint estimation of multiple networks. While GGL emphasizes the sharing of common edge support across networks, CNJGL focuses on encouraging shared node support. Consequently, the CNJGL model identifies common hub nodes shared across all networks but cannot detect hub nodes unique to specific networks. Comparing the NJGCG model with the CNJGL model will elucidate the effectiveness of our NJGCG in detecting both common and specific hub nodes. Additionally, neither the GGL nor the CNJGL model addresses issues related to non-Gaussian data distributions and missing values. Although JEGN can handle non-Gaussian data and capture both shared and specific subnetworks, it does not account for hub genes and fails to address the issue of missing data. To evaluate the advantages of accounting for missing values, we also present experimental results of our model without accounting for missing values (denoted as NJGCG-N). Specifically, for GGL, CNJGL, JEGN and NJGCG-N, missing values are set to zeros.

Both GGL and CNJGL have two parameters λ1 and λ2, where λ1 controls the sparsity of the estimated network, λ2 of GGL captures the similarities between different networks, and λ2 of CNJGL controls the hub structures in networks. Similar to [30], we reset the tuning parameters of GGL to ω1=λ1+(12)λ2 and ω2=(12)λ2λ1+(12)λ2. For JEGN, there are two parameters λ and α, where λ controls the sparsity of the estimated networks and α controls the balance between shared and specific subnetworks. For NJGCG and NJGCG-N, λ1 and λ2 control the sparsity and hub structure of the estimated networks, respectively. As λ2 can also affect the sparsity of the estimated networks, we set λ2=λrλ1. Then the tuning parameters become λ1 and λr. Here, the ranges for tuning the parameters of the GGL model are set to ω1∈[0.001,10] and ω2∈{0.5,0.7,1.0}. The ranges for tuning the parameters of the CNJGL model are set to λ1∈[0.001,10] and λ2∈{1,5,10}. The ranges for tuning the parameters of the JEGN model are set to λ∈[0.001,10] and α∈{0.1,0.3,0.5}. The ranges for tuning the parameters of the NJGCG and NJGCG-N models are set to λ1∈[0.001,10] and λr∈{2.5,3.0,3.5}.

The performance of various methods is evaluated using the ROC curve. Let ζˆij(k) denote the (i,j)-th entry of an estimator Θˆ(k) and ζij(k) denote the (i,j)-th entry of the ground truth Θ(k), for k=1,…,K, the true positive rate (TPR) and false positive rate (FPR) are defined as follows:TPR=∑k=1K∑i<jI(ζˆij(k)≠0andζij(k)≠0)∑k=1K∑i<jI(ζij(k)≠0),

FPR=∑k=1K∑i<jI(ζˆij(k)≠0andζij(k)=0)∑k=1K∑i<jI(ζij(k)=0),

where I{⋅} is the indicator function. The ROC curve is used to evaluate the performance as a function of the tuning parameters that controls the sparsity level of the estimated networks.

Fig. 3, Fig. 4 show the results of various methods on simulated Gaussian data of scale-free and Erdos–Renyi networks, respectively. Fig. 5, Fig. 6 show the results of various methods on simulated non-Gaussian data of scale-free and Erdos–Renyi networks, respectively. The results of all methods are averaged over 5 repetitions. The columns correspond to different missing rates δ, and the rows correspond to different sample sizes nk. Within each subplot, each colored line corresponds to the results of a method as the value of tuning parameter changes, and each point on the line represents the result for a specific value of the tuning parameter.Fig. 3 Performance of different methods in terms of the ROC curves on simulated Gaussian data of scale-free networks. The columns correspond to different missing data ratio δ, and the rows correspond to different sample sizes nk. In each plot, each colored line corresponds to the result of a method as the tuning parameter is varied.

Fig. 3

Fig. 4 Performance of different methods in terms of the ROC curves on simulated Gaussian data of Erdos–Renyi networks. The columns correspond to different missing data ratio δ, and the rows correspond to different sample sizes nk. In each plot, each colored line corresponds to the result of a method as the tuning parameter is varied.

Fig. 4

Fig. 5 Performance of different methods in terms of the ROC curves on simulated non-Gaussian data of scale-free networks. The columns correspond to different missing data ratio δ, and the rows correspond to different sample sizes nk. In each plot, each colored line corresponds to the result of a method as the tuning parameter is varied.

Fig. 5

Fig. 6 Performance of different methods in terms of the ROC curves on simulated non-Gaussian data of Erdos–Renyi networks. The columns correspond to different missing data ratio δ, and the rows correspond to different sample sizes nk. In each plot, each colored line corresponds to the result of a method as the tuning parameter is varied.

Fig. 6

We can find from Fig. 3, Fig. 4 that our NJGCG model consistently outperforms GGL, CNJGL, and JEGN across all scenarios. This highlights the superior effectiveness of NJGCG in estimating multiple networks that include both common and specific hub nodes. By accounting for missing data, NJGCG outperforms NJGCG-N, suggesting that handling missing values are more effective than simply treating them as zeros. We can also find that as the proportion of missing data increases, the performance advantage of our method becomes more pronounced, demonstrating the necessity of accounting for missing data.

As shown in Fig. 5, Fig. 6, for non-Gaussian data, all methods experience a decline in performance, but our NJGCG model still outperforms all the comparative methods. We found that the performance of GGL and CNJGL declines most noticeably, demonstrating the importance of considering data distribution.

As a graphical model, our NJGCG model has relatively high computational complexity. As the number of nodes p increases, the computation time of our model also increases. To present this more clearly, we varied the number of nodes from 100 to 1000 and measured the corresponding computation times, which are shown in Fig. 7. The experiments were conducted on a Windows PC equipped with an Intel i9-10900K CPU @ 3.70GHz. The results indicate that the computation time of our algorithm grows with the number of nodes. However, since network inference is usually performed offline, our primary focus is on improving the accuracy of network inference. Thus, the trade-off of higher computational complexity for enhanced inference accuracy is justified.Fig. 7 The running time of our NJGCG model for different numbers of nodespon simulated Gaussian data of scale-free networks. The x-axis represents different numbers of p ∈ {100,200,…,1000}. The y-axis represents the running time.

Fig. 7

4 Real data application

4.1 Mouse embyonic stem-cell (mESC) differentiation

During the differentiation of mouse embryonic stem cells (mESC), various regulatory factors govern the entire process, and these factors vary at different stages of differentiation, leading to changes in the gene network of mESC throughout the differentiation processes. Consequently, our objective is to reconstruct the gene networks at different stages of mESC differentiation. To achieve this, we utilize single-cell RNA sequencing (scRNA-seq) data [28] obtained from four time points (days 0, 2, 4, and 7) during the differentiation of mouse embryonic stem cells into epithelial cells. Assuming that the daily data represents the dominant cell state at different stages of differentiation, the dataset quantifies gene expression at different stages. The dataset (downloaded from GEO database with accession GSE65525) comprises 2717 cells and 24175 genes, with the number of cells on the 0th, 2nd, 4th, and 7th days being 933, 303, 683, and 798, respectively. Given our focus on exploring the differentiation processes of mouse embryonic stem cells, following previous studies [63], [64], [41], we concentrate on a list of 84 genes, including essential housekeeping genes, key regulators of mESC differentiation, and differentiation markers. The proportions of zero values in these four gene expression datasets are 42.84%, 82.53%, 88.88% and 52.32% respectively.

We applied our NJGCG model to the aforementioned scRNA-seq data of mouse embryonic stem cells. Considering that the gene expression data does not follow a Gaussian distribution and contains missing values, we employed the improved Kendall's tau estimator (Eq. (3)) to calculate the sample covariance matrix. Subsequently, we utilized our NJGCN model to estimate the gene networks.

To ensure interpretable results, we reset (ω1,ω2,ω3)=(1,ln⁡p,ln⁡(p(K−1))). The BIC criterion, as introduced in Eq. (16), is employed to determine the values of λ1 and λ2, and the selected values are λ1=0.3 and λ2=0.276. Our estimated gene networks are shown in Fig. 8, where genes closer to the center have a higher number of connections.Fig. 8 The gene networks of mESC differentiation estimated by NJGCG. Genes that are closer to the center have more connections. Hub nodes analyzed in the article are highlighted in green. (A) Day 0, (B) Day 2, (C) Day4, (D) Day 7.

Fig. 8

We primarily focus on genes with higher degrees in the estimated networks, as they may play a crucial role in regulating the entire differentiation process [50]. Our NJGCG model can identify both common hub genes and specific hub genes across networks at different stages. In this study, genes with degrees greater than twice the average degree of the network are considered hub genes and are highlighted in green in Fig. 8. Dppa5a, Pou5f1, Zfp42, and Cd9 are identified as hub genes in the first stage, while Dpp5a, Pou5f1, Zfp42, Lama1, Lamc1, and Fn1 are identified as hub genes in the second stage. Pou5f1, Lama1, and Lamc1 are identified as hub genes in the third stage, and Cd9, Lamc1, and Fn1 are identified as hub genes in the fourth stage.

Dppa5a, identified as a hub gene in the first two stages, is known to be related to pluripotency development, essential for the self-renewal of pluripotent embryonic stem cells and the generation of germ cells [65]. The POU transcription factor Pou5f1 (also known as Oct4), identified as a hub gene in the first three stages, is a key regulator in pluripotent stem cells and is critical for the formation of a pluripotent basal cell population in mammalian embryos [66]. Research has shown that Zfp42 (also called Rex1) can regulate the differentiation of pluripotent stem cells along different cell lineages found in early embryos [67]. Cd9 encodes a transmembrane protein that plays a role in many cellular processes, including differentiation, adhesion, and signal transduction [68]. Lama1 encodes the production of laminin α-1, which participates in the formation of embryos and is essential for early embryo growth [69]. Lamc1 encodes the subunit of laminin γ1, and its absence can impact the differentiation of embryonic stem cells [70]. Fn1 encodes endogenous fibronectin, which is essential for the self-renewal of mouse embryonic stem cells [71].

In summary, the results demonstrate the effectiveness of our NJGCG model in identifying important genes that contribute to understanding the functional processes of cell differentiation.

4.2 TCGA breast cancer

Breast cancer is a complex disease characterized by different molecular subtypes, including luminal A, luminal B, HER2-enriched, and basal-like [72]. However, limited research has been conducted to explore the relationships between different breast cancer subtypes based on their gene networks. To gain insights into the characteristics of breast cancer subtypes from the perspective of network rewiring, we aim to infer the gene networks of different cancer subtypes. To achieve this, we obtain RNA sequencing data of breast cancer patients (level 3) from the TCGA database (version: May 6, 2017). This dataset comprises gene expression data from 489 patients with primary breast cancer across 15,718 genes. The subtype information of the patients is obtained from [72], revealing 214 luminal A subtype, 117 luminal B subtype, 55 HER2-enriched subtype, and 90 basal-like subtype patients. In this study, we focus on the breast cancer pathway (hsa05224), which can be obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database [73], resulting in 139 common genes.

The data is preprocessed using the method described in the above section. In order to produce interpretable results, we still set (ω1,ω2,ω3)=(1,ln⁡p,ln⁡(p(K−1))) and use the BIC rule in Eq. (16) to select the values of λ1 and λ2. Finally, the selected values are λ1=0.6 and λ2=0.63. The estimated gene networks for different subtypes are shown as Fig. 9.Fig. 9 The gene networks of different breast cancer subtypes estimated by NJGCG. Genes that are closer to the center have more connections. Hub nodes analyzed in the article are highlighted in green. (A) Luminal A, (B) Luminal B, (C) HER2-enriched, (D) Basal-like.

Fig. 9

We also analyze the biological significance of the identified hub genes in the estimated networks. Here, genes with degrees greater than twice the average degree of the network are considered hub genes and are highlighted in green in Fig. 9.

The hub genes APC and FGF7 are common across all four cancer subtypes. The Adenomatous Polyposis Coli (APC) tumor suppressor is mutated or hypermethylated in up to 70% of sporadic breast cancers, depending on the subtype [74]. Additionally, studies suggest that Linc00460 sponges endogenous miRNA-489-5p to activate FGF7-AKT signaling, thereby promoting breast cancer progression [75]. Epidermal Growth Factor Receptor (EGFR) is identified as a hub gene in Luminal A, Luminal B, and Basal-like subtypes. Over-expression of EGFR has been linked to the development of various cancers, including breast cancer, and it is a well-known feature of molecular subtypes and a useful therapeutic target [72], [76]. AKT3 is identified as a hub gene in Luminal A, Luminal B, and HER2-enriched subtypes. Research has shown that regulating the expression of miRNA can enhance the expression of AKT3 to inhibit breast cancer [77]. Except for Luminal B, SOS1 is identified as a hub gene. Studies have found that SOS1 is part of the EGFR-dependent pathway, and its over-expression promotes cell growth and survival, leading to drug resistance and increasing the activation of the quiescent nuclear factor κB (NFκB) [78]. Fibroblast Growth Factor 23 (FGF23) is only identified as a hub gene in the HER2-enriched subtype. The high interaction of FGF-FGFR in it may be the cause of breast cancer [79]. Overall, the results demonstrate the effectiveness of our NJGCG model in identifying important genes that can aid in understanding the underlying mechanisms of different breast cancer subtypes.

5 Discussion and conclusion

In this study, we introduce a novel node-based joint graphical model to simultaneously infer gene networks across multiple states. The novelty of our model lies in its capability to identify common and specific hub genes among different gene networks. Specifically, by incorporating a row-column overlap norm penalty and a tree-structured group Lasso penalty on the networks, our model facilitates the identification of common and specific hub genes in the estimated gene networks. Additionally, our model is adept at handling non-Gaussian data with missing values. Simulation studies have validated the effectiveness of our model in inferring the underlying gene networks of different states. Furthermore, we have applied our model to infer gene networks at various stages of stem cell differentiation and for different breast cancer subtypes. The hub genes we identified, which are common to all networks or specific to certain networks, offer valuable insights into the differentiation processes of stem cells and the heterogeneity among different cancer subtypes.

In addition to graphical models, methods based on mutual information have been widely utilized for reconstructing gene networks, as they have the ability to quantify the nonlinear relationships between genes [6], [8], [16]. However, mutual information-based methods are limited in their capacity to reconstruct only a single gene network from a single dataset, making it challenging for them to jointly estimate multiple networks. Furthermore, graphical model-based methods are primarily suited for measuring linear relationships. Hence, the development of new methods that can simultaneously explore nonlinear relationships and jointly estimate multiple networks represents a potential future direction.

The primary focus of this study is to identify dependency relationships between genes. In addition to genes, various other biomolecules, such as microRNA and long non-coding RNA, play crucial roles in regulatory processes. Different omics data typically adhere to distinct probability distributions. For instance, gene mutations are commonly modeled by the Bernoulli distribution. The approach we have employed to handle missing values can be adapted to accommodate data following different distributions. In the future, we aim to extend our method to effectively handle multi-omics data.

Supplementary data and code

The source codes of our NJGCG model are available at https://github.com/HuangY2022/NJGCG.

CRediT authorship contribution statement

Yun Huang: Writing – original draft, Methodology, Investigation, Conceptualization. Sen Huang: Writing – original draft, Validation, Methodology, Investigation, Data curation, Conceptualization. Xiao-Fei Zhang: Writing – review & editing, Validation, Methodology. Le Ou-Yang: Writing – review & editing, Supervision, Methodology, Investigation, Conceptualization. Chen Liu: Writing – review & editing, Supervision, Investigation, Conceptualization.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

This work was supported by the 10.13039/501100003392 Natural Science Foundation of Fujian Province , China (Grant No. 2020J01956 , 2022J01214 ), Fujian Traditional Chinese Medicine Science and Technology Project of 10.13039/501100014125 Fujian Provincial Health Commission (No. 2021zyyj68 ), Joint Funds for the Innovation of Science and Technology, Fujian Province (No. 2021Y9102 ).
==== Refs
References

1 Barabási A.-L. Gulbahce N. Loscalzo J. Network medicine: a network-based approach to human disease Nat Rev Genet 12 1 2011 56 68 21164525
2 Marbach D. Costello J.C. Küffner R. Vega N.M. Prill R.J. Camacho D.M. Wisdom of crowds for robust gene network inference Nat Methods 9 8 2012 796 804 22796662
3 Hill S.M. Heiser L.M. Cokelaer T. Unger M. Nesser N.K. Carlin D.E. Inferring causal molecular networks: empirical assessment through a community-based effort Nat Methods 13 4 2016 310 318 26901648
4 You Z.-H. Zhou M. Luo X. Li S. Highly efficient framework for predicting interactions between proteins IEEE Trans Cybern 47 3 2016 731 743 28113829
5 Wang Y. Joshi T. Zhang X.-S. Xu D. Chen L. Inferring gene regulatory networks from multiple microarray datasets Bioinformatics 22 19 2006 2413 2420 16864593
6 Zhang X. Zhao X.-M. He K. Lu L. Cao Y. Liu J. Inferring gene regulatory networks from gene expression data by path consistency algorithm based on conditional mutual information Bioinformatics 28 1 2012 98 104 22088843
7 Wang Z. Xu W. San Lucas F.A. Liu Y. Incorporating prior knowledge into gene network study Bioinformatics 29 20 2013 2633 2640 23956306
8 Zhao J. Zhou Y. Zhang X. Chen L. Part mutual information for quantifying direct associations in networks Proc Natl Acad Sci 113 18 2016 5130 5135 27092000
9 Liu X. Wang Y. Ji H. Aihara K. Chen L. Personalized characterization of diseases using sample-specific networks Nucleic Acids Res 44 22 2016 e164 27596597
10 Şenbabaoğlu Y. Sümer S.O. Sanchez-Vega F. Bemis D. Ciriello G. Schultz N. A multi-method approach for proteomic network inference in 11 human cancers PLoS Comput Biol 12 2 2016 e1004765
11 Zhang Y. Ouyang Z. Zhao H. A statistical framework for data integration through graphical models with application to cancer genomics Ann Appl Stat 11 1 2017 161 30956747
12 Ou-Yang L. Zhang X.-F. Wu M. Li X.-L. Node-based learning of differential networks from multi-platform gene expression data Methods 129 2017 41 49 28579401
13 Ou-Yang L. Yan H. Zhang X.-F. Identifying differential networks based on multi-platform gene expression data Mol BioSyst 13 1 2017 183 192
14 Langfelder P. Horvath S. Wgcna: an r package for weighted correlation network analysis BMC Bioinform 9 1 2008 1 13
15 Zhang J. Huang K. Normalized imqcm: an algorithm for detecting weak quasi-cliques in weighted graph with applications in gene co-expression module discovery in cancers Cancer Inform 13 2014 CIN–S14021
16 Zhang X. Zhao J. Hao J.-K. Zhao X.-M. Chen L. Conditional mutual inclusive information enables accurate quantification of associations in gene regulatory networks Nucleic Acids Res 43 5 2015 e31 25539927
17 Friedman J. Hastie T. Tibshirani R. Sparse inverse covariance estimation with the graphical lasso Biostatistics 9 3 2008 432 441 18079126
18 Koller D. Friedman N. Probabilistic graphical models: principles and techniques 2009 MIT Press
19 Yuan M. Lin Y. Model selection and estimation in the Gaussian graphical model Biometrika 94 1 2007 19 35
20 Zhao S.D. Cai T.T. Li H. Direct estimation of differential networks Biometrika 101 2 2014 253 268 26023240
21 Zhang X.-F. Ou-Yang L. Yang S. Hu X. Yan H. Diffgraph: an r package for identifying gene network rewiring using differential graphical models Bioinformatics 34 9 2018 1571 1573 29309511
22 Zhang X.-F. Ou-Yang L. Yang S. Hu X. Yan H. Diffnetfdr: differential network analysis with false discovery rate control Bioinformatics 35 17 2019 3184 3186 30689728
23 Banerjee O. El Ghaoui L. d'Aspremont A. Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data J Mach Learn Res 9 2008 485 516
24 Dalal O. Rajaratnam B. Sparse Gaussian graphical model estimation via alternating minimization Biometrika 104 2 2017 379 395
25 Liu H. Lafferty J. Wasserman L. The nonparanormal: semiparametric estimation of high dimensional undirected graphs J Mach Learn Res 10 2009 2295 2328
26 Liu H. Han F. Yuan M. Lafferty J. Wasserman L. High-dimensional semiparametric Gaussian copula graphical models Ann Stat 40 4 2012 2293 2326
27 Xue L. Zou H. Regularized rank-based estimation of high-dimensional nonparanormal graphical models Ann Stat 40 5 2012 2541 2571
28 Klein A.M. Mazutis L. Akartuna I. Tallapragada N. Veres A. Li V. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells Cell 161 5 2015 1187 1201 26000487
29 Mukherjee S. Carignano A. Seelig G. Lee S.-I. Identifying progressive gene network perturbation from single-cell rna-seq data 2018 40th annual international conference of the IEEE engineering in medicine and biology society (EMBC) 2018 IEEE 5034 5040
30 Danaher P. Wang P. Witten D.M. The joint graphical lasso for inverse covariance estimation across multiple classes J R Stat Soc, Ser B, Stat Methodol 76 2 2014 373
31 Lee W. Liu Y. Joint estimation of multiple precision matrices with common structures J Mach Learn Res 16 1 2015 1035 1062 26568704
32 Guo J. Levina E. Michailidis G. Zhu J. Joint estimation of multiple graphical models Biometrika 98 1 2011 1 15 23049124
33 Wang B. Singh R. Qi Y. A constrained ℓ 1 minimization approach for estimating multiple sparse Gaussian or nonparanormal graphical models Mach Learn 106 9 2017 1381 1417
34 Ma J. Michailidis G. Joint structural estimation of multiple graphical models J Mach Learn Res 17 1 2016 5777 5824
35 Cai T.T. Li H. Liu W. Xie J. Joint estimation of multiple high-dimensional precision matrices Stat Sin 26 2 2016 445 28316451
36 Huang F. Chen S. Joint learning of multiple sparse matrix Gaussian graphical models IEEE Trans Neural Netw Learn Syst 26 11 2015 2606 2620 25751876
37 Huang F. Chen S. Huang S.-J. Joint estimation of multiple conditional Gaussian graphical models IEEE Trans Neural Netw Learn Syst 29 7 2017 3034 3046 28678717
38 Huang F. Chen S. Learning dynamic conditional Gaussian graphical models IEEE Trans Knowl Data Eng 30 4 2017 703 716
39 Zhang X.-F. Ou-Yang L. Yan T. Hu X.T. Yan H. A joint graphical model for inferring gene networks across multiple subpopulations and data types IEEE Trans Cybern 51 2 2019 1043 1055
40 Tan K.M. London P. Mohan K. Lee S.-I. Fazel M. Witten D. Learning graphical models with hubs J Mach Learn Res 15 2014 3297 3331 25620891
41 Wu N. Yin F. Ou-Yang L. Zhu Z. Xie W. Joint learning of multiple gene networks from single-cell gene expression data Comput Struct Biotechnol J 18 2020 2583 2595 33033579
42 Eckstein J. Bertsekas D.P. On the Douglas—Rachford splitting method and the proximal point algorithm for maximal monotone operators Math Program 55 1 1992 293 318
43 Boyd S. Parikh N. Chu E. Distributed optimization and statistical learning via the alternating direction method of multipliers 2011 Now Publishers Inc.
44 Eckstein J. Yao W. Augmented Lagrangian and alternating direction methods for convex optimization: a tutorial and some illustrative computational results RUTCOR Res Rep 32 3 2012 44
45 Ou-Yang L. Cai D. Zhang X.-F. Yan H. Wdne: an integrative graphical model for inferring differential networks from multi-platform gene expression data with missing values Brief Bioinform 22 6 2021 bbab086
46 Wang H. Fazayeli F. Chatterjee S. Banerjee A. Gaussian copula precision estimation with missing values Artificial intelligence and statistics PMLR 2014 978 986
47 Fang H.-B. Fang K.-T. Kotz S. The meta-elliptical distributions with given marginals J Multivar Anal 82 1 2002 1 16
48 Kruskal W.H. Ordinal measures of association J Am Stat Assoc 53 284 1958 814 861
49 Mohan K. London P. Fazel M. Witten D. Lee S.-I. Node-based learning of multiple Gaussian graphical models J Mach Learn Res 15 1 2014 445 488 25309137
50 Hao D. Ren C. Li C. Revisiting the variation of clustering coefficient of biological networks suggests new modular structure BMC Syst Biol 6 1 2012 1 10 22222070
51 Kim S. Xing E.P. Tree-guided group lasso for multi-task regression with structured sparsity ICML 2010
52 Kim S. Xing E.P. Tree-guided group lasso for multi-response regression with structured sparsity, with an application to eqtl mapping Ann Appl Stat 2012 1095 1117
53 Liu J. Ye J. Moreau-Yosida regularization for grouped tree structure learning Adv Neural Inf Process Syst 23 2010 1459 1467
54 Liu J. Ji S. Ye J. Slep: sparse learning with efficient projections Ariz State Univ 6 491 2009 7
55 Xu T. Ou-Yang L. Yan H. Zhang X.-F. Time-varying differential network analysis for revealing network rewiring over cancer progression IEEE/ACM Trans Comput Biol Bioinform 18 4 2019 1632 1642
56 Yuan M. Lin Y. Model selection and estimation in the Gaussian graphical model Biometrika 94 1 2007 19 35
57 Zou H. Hastie T. Tibshirani R. On the “degrees of freedom” of the lasso Ann Stat 35 5 2007 2173 2192
58 Yuan M. Lin Y. Model selection and estimation in regression with grouped variables J R Stat Soc, Ser B, Stat Methodol 68 1 2006 49 67
59 Barabási A.-L. Albert R. Emergence of scaling in random networks Science 286 5439 1999 509 512 10521342
60 Barabási A.-L. Scale-free networks: a decade and beyond Science 325 5939 2009 412 413 19628854
61 Chen Y. Zhang X.-F. Ou-Yang L. Inferring cancer common and specific gene networks via multi-layer joint graphical model Comput Struct Biotechnol J 21 2023 974 990 36733706
62 Grimes T. Datta S. Seqnet: an r package for generating gene-gene networks and simulating rna-seq data J Stat Softw 98 12 2021
63 Przybyla L.M. Voldman J. Attenuation of extrinsic signaling reveals the importance of matrix remodeling on maintenance of embryonic stem cell self-renewal Proc Natl Acad Sci 109 3 2012 835 840 22215601
64 Mukherjee S. Carignano A. Seelig G. Lee S.-I. Identifying progressive gene network perturbation from single-cell rna-seq data 2018 40th annual international conference of the IEEE engineering in medicine and biology society (EMBC) 2018 IEEE 5034 5040
65 Amano H. Itakura K. Maruyama M. Ichisaka T. Nakagawa M. Yamanaka S. Identification and targeted disruption of the mouse gene encoding esg1 (ph34/ecat2/dppa5) BMC Dev Biol 6 1 2006 1 9 16412219
66 Niwa H. Miyazaki J.-i. Smith A.G. Quantitative expression of oct-3/4 defines differentiation, dedifferentiation or self-renewal of es cells Nat Genet 24 4 2000 372 376 10742100
67 Thompson J.R. Gudas L.J. Retinoic acid induces parietal endoderm but not primitive endoderm and visceral endoderm differentiation in f9 teratocarcinoma stem cells with a targeted deletion of the rex-1 (zfp-42) gene Mol Cell Endocrinol 195 1–2 2002 119 133 12354678
68 Akutsu H. Miura T. Machida M. Birumachi J.-i. Hamada A. Yamada M. Maintenance of pluripotency and self-renewal ability of mouse embryonic stem cells in the absence of tetraspanin cd9 Differentiation 78 2–3 2009 137 142 19716222
69 Li X. Chen Y. Schéele S. Arman E. Haffner-Krausz R. Ekblom P. Fibroblast growth factor signaling and basement membrane assembly are connected during epithelial morphogenesis of the embryoid body J Cell Biol 153 4 2001 811 822 11352941
70 Fujiwara H. Hayashi Y. Sanzen N. Kobayashi R. Weber C.N. Emoto T. Regulation of mesodermal differentiation of mouse embryonic stem cells by basement membranes J Biol Chem 282 40 2007 29701 29711 17690109
71 Hunt G.C. Singh P. Schwarzbauer J.E. Endogenous production of fibronectin is required for self-renewal of cultured mouse embryonic stem cells Exp Cell Res 318 15 2012 1820 1831 22710062
72 Network C.G.A. Comprehensive molecular portraits of human breast tumours Nature 490 7418 2012 61 23000897
73 Kanehisa M. Goto S. Kegg: Kyoto encyclopedia of genes and genomes Nucleic Acids Res 28 1 2000 27 30 10592173
74 VanKlompenberg M.K. Bedalov C.O. Soto K.F. Prosperi J.R. Apc selectively mediates response to chemotherapeutic agents in breast cancer BMC Cancer 15 1 2015 1 14 25971837
75 Zhu Y. Yang L. Chong Q.-Y. Yan H. Zhang W. Qian W. Long noncoding rna linc00460 promotes breast cancer progression by regulating the mir-489-5p/fgf7/akt axis Cancer Manag Res 11 2019 5983 31308741
76 Schnitt S.J. Classification and prognosis of invasive breast cancer: from morphology to molecular taxonomy Mod Pathol 23 2 2010 S60 S64 20436504
77 Ruan L. Qian X. Mir-16-5p inhibits breast cancer by reducing akt3 to restrain nf-κb pathway Biosci Rep 39 8 2019
78 De S. Dermawan J.K.T. Stark G.R. Egf receptor uses sos1 to drive constitutive activation of nfκb in cancer cells Proc Natl Acad Sci 111 32 2014 11721 11726 25071181
79 Cekin R. Arici S. Atci M.M. Secmeler S. Cihan S. The clinical importance of fibroblast growth factor 23 on breast cancer patients J Med Invest 4 2020 471 476
