 Research
 Open Access
 Published:
Networkbased cancer genomic data integration for pattern discovery
BMC Genomic Data volume 22, Article number: 54 (2021)
Abstract
Background
Since genes involved in the same biological modules usually present correlated expression profiles, lots of computational methods have been proposed to identify gene functional modules based on the expression profiles data. Recently, Sparse Singular Value Decomposition (SSVD) method has been proposed to bicluster gene expression data to identify gene modules. However, this model can only handle the gene expression data where no gene interaction information is integrated. Ignoring the prior gene interaction information may produce the identified gene modules hard to be biologically interpreted.
Results
In this paper, we develop a Sparse Networkregularized SVD (SNSVD) method that integrates a prior gene interaction network from a protein protein interaction network and gene expression data to identify underlying gene functional modules. The results on a set of simulated data show that SNSVD is more effective than the traditional SVDbased methods. The further experiment results on real cancer genomic data show that most coexpressed modules are not only significantly enriched on GO/KEGG pathways, but also correspond to dense subnetworks in the prior gene interaction network. Besides, we also use our method to identify ten differentially coexpressed miRNAgene modules by integrating matched miRNA and mRNA expression data of breast cancer from The Cancer Genome Atlas (TCGA). Several important breast cancer related miRNAgene modules are discovered.
Conclusions
All the results demonstrate that SNSVD can overcome the drawbacks of SSVD and capture more biologically relevant functional modules by incorporating a prior gene interaction network. These identified functional modules may provide a new perspective to understand the diagnostics, occurrence and progression of cancer.
Background
With the rapid development of (singlecell) RNASeq and microarray technologies, huge number of cancer genomic data have been generated [1–3]. The data provide some new opportunities to study on the gene cooperative mechanisms [4–8]. Based on the hypothesis that genes with similar functions may show similar expression patterns, clustering techniques have been used to identify coexpressed gene sets in which genes present similar expression patterns across all samples. However, these traditional clustering techniques face with the limitation that some genes can coregulate across some samples rather than all samples in the real biological systems [9]. Therefore, many biclustering methods [4, 10–13] are proposed to discover some coexpressed gene sets in which genes present similar expression patterns across some samples.
Recently, several Sparse Singular Value Decomposition (SSVD) based methods have been proposed for biclustering gene expression data to discover gene functional modules (biclusters) [14], such as ALSVD [4], L0SVD [15], and so on. However, most of them ignore the prior gene interaction network knowledge from a protein protein interaction (PPI) network, whereas such information is very useful to improve biological interpretability of discovered gene modules [16–18]. The PPI network has been used in many biological applications for accurate discovery or better biological interpretability [16, 19–22]. However, as far as we know, there is very little work to incorporate the gene interaction network knowledge from PPI network into the biclustering framework. To address it, we integrate the gene network in the SSVD model for biclustering gene expression [23].
In this paper, we develop a sparse networkregularized SVD (SNSVD) model to identify gene functional modules by integrating gene expression data and a prior gene interaction network from a PPI network (Fig. 1). To ensure the discovered gene modules in which genes are coexpressed and densely connected in the prior PPI network, we introduce a sparse networkregularized penalty [20] in the model. Compared with the traditional regularized penalties (e.g., LASSO [24]), the sparse networkregularized penalty can make the biclustering process tend to select correlated and interacted genes for enhancing biological interpretability of gene modules. We present an alternating iterative algorithm to solve the SNSVD model. We evaluate the performance of SNSVD using a set of artificial data sets and gene expression data from the Cancer Genome Project (CGP) [3], and compare its performance with other representative SSVD methods. We investigate the functionality of these identified modules from multiple perspectives. The results show that SNSVD can identify more biologically relevant gene sets and improve their biological interpretations.
Additionally, we present a framework based on SNSVD for analyzing matched miRNA and mRNA expression data to identify differentially coexpressed miRNAgene modules. Extensive results when applying onto TCGA breast cancer data demonstrate that the identified miRNAgene modules provide a new perspective for diagnosis and discrimination between two status of breast cancer.
Results
Simulation study
We evaluated the performance of SNSVD on the simulated data by comparing it with other sparse SVD based methods including L0SVD [15], ALSVD [4] and SCADSVD [4, 25]. Without loss of generality, we define a rankone true signal matrix as uv^{T} where u and v are vectors of size p×1 and n×1, respectively. The observed matrix is defined as X=uv^{T}+γε, where ε is a noise matrix each element in which is randomly sampled from a standard normal distribution and γ is a nonnegative parameter to control the signaltonoise ratio (SNR).
To generate the simulated data, we first generated two sparse singular vectors u and v with p=200,n=100 whose first 50 elements equal to 1 (nonzeros), and the remaining ones are zeros. Then we created a series of observation matrices X for each γ ranging from 0.02 to 0.06 in steps of 0.005. In addition, we created a prior “PPI” network for row variables of X, where any two nodes in first 50 vertices are connected with probability p_{11}=0.3, and remaining ones are connected with probability p_{12}=0.1. For each γ, we generated 50 different noise matrices ε to got 50 observed matrices X for testing. The average sensitivity, specificity and accuracy of u (or v) on the 50 matrices X were calculated. Moreover, we set σ=0.5 according to 5fold cross validation test, and forced u and v to contain 50 nonzero elements with same sparsity level for each method by tuning the parameters so that the results of different methods are comparable. The average sensitivities, specificities and accuracies of u (or v) with different γ were compared in Fig. 2. We found that the performance of our proposed method (SNSVD) is superior to that of other methods. The results illustrate that SNSVD model can enhance the power of variable selection by integrating the prior network knowledge.
Application to the CGP gene expression data sets
We further investigated the performance of SNSVD by using the gene expression data with 641 cell lines including diverse cancer types and tissues from the Cancer Genome Project (CGP) (http://www.cancerrxgene.org/downloads) [3], and a PPI network from the PathwayCommons database [26]. In total, there are 13,321 genes and 262,462 interactions in the PPI network. The 641 cell lines are from 16 tissues or 52 cancer types in the CGP data, where a tissue type contains about 40 cell lines and a cancer type contains about 12 cell lines.
Identifying functional modules
We set σ=100 according to 5fold cross validation test, and set k_{v}=50 (control the sample sparsity). In addition, we also selected a suitable λ to force the estimated u only containing 200 nonzero elements (control the gene sparsity). Using Algorithm 3, we identified the first 40 pairs of singular vectors {(u_{1},v_{1}),⋯,(u_{40},v_{40})}. Let U=[u_{1};⋯ ;u_{40}] and V=[v_{1};⋯ ;v_{40}], where the ith column of U and V correspond to the ith pair of sparse singular vectors. To reduce the false positive cases, we first calculated absolute zscore for each column of U (or V) according to Eq. (1). For each nonzero x_{ij}, we define the following formula:
where x_{ij} is ith element in u_{j} (or v_{j}), μ_{j} is the average value of all nonzero elements in u_{j} (or v_{j}), and σ_{j} is their standard deviation. If z_{ij} is greater than a given threshold, the corresponding gene (or sample) is then selected into the module j. Herein, we obtained 40 gene functional modules with 160 genes and 40 samples in average (Fig. 3).
Functional analysis of the genes in modules
Firstly, we investigated whether the genes within the same modules are significantly coexpressed by calculating the modularity score in Eq. (16), the result showed that all identified modules were statistically significant with pvalue <0.01 by using onesided Wilcoxon signed rank test (Fig. 4).
Secondly, we also investigated whether the genes within the same modules are connected with each other in the prior PPI network via the genegene interaction enrichment score. The result showed that 57% of the 40 modules were significantly interconnected with each other in the PPI network, illustrating that our method tend to cluster genes interacting with each other.
Finally, we also checked the biological relevance of all the identified gene modules using gene functional enrichment analysis via DAVID online web tool [27]. By selecting the GO BP (Gene Ontology Biological Process) and KEGG pathways with BenjaminiHochberg adjusted pvalues<0.05 as significant ones, we obtained 766 significant GO BP pathways and 70 significant KEGG pathways. By statistically, 62.5% modules are significantly related with at least one GO BP pathways and 42.5% modules are significantly related with at least one KEGG pathways.
Functional analysis of the samples in modules
To evaluate the subtypespecific of samples in the identified modules, we computed the overlapping significance level of between modulesamples and cancer/tissue specific samples. For each gene module, we first collected a sample set from the module. We then computed the overlapping significance levels between the sample set and any one tissuesample set using the right hypergeometric test (Fig. 5A), and the overlapping significance levels between the sample set and any one cancersample set (Fig. 5B). We found that most of the identified gene modules can be seen as subtypespecific gene functional modules, which provide insights into the mechanisms of the relationship between different tissues and cancers.
Additionally, we also found that the cancer/tissue types of some modules are consistent with their corresponding functional pathways. Some examples are listed in Table 1 (See also Fig. 5). For example, module 1 contains 47 cell lines significantly overlapping with blood tissue and some bloodrelated cancers (e.g., AML, B cell leukemia, B cell lymphoma, lymphoblastic leukemia, lymphoblastic T cell leukaemia, lymphoid_neoplasm other), while the top enriched GO/KEGG pathways of 174 genes in module 1 are related to the immune system. Some previous work have reported that the development of bloodrelated cancers are associated with immune pathway abnormalities [28, 29]. Similarly, these samples in module 2 are also significantly related with some bloodrelated cancers (B cell leukemia, B cell lymphoma, Burkitt lymphoma, lymphoblastic leukemia, and lymphoid_neoplasm other), while some genes in which are significantly enriched in some immunerelated pathways. These samples in module 4 are significantly related with central nervous system (CNS), while some genes in which are significantly enriched in nervous system related GO/KEGG pathways.
Finally, we also evaluated whether the identified 40 modules are greatly overlapped. Since each module contains a gene set and a sample set. To assess the overlapping relationship between two different modules. For any two gene modules, we computed the overlapping significance level p_{1} and p_{2} between their gene sets and sample sets respectively by using the righttailed hypergeometric test. If p_{1}<0.05 and p_{2}<0.05, then we considered that the two modules are significant overlapped. Among all 780 modulemodule pairs for the identified 40 modules, we found that only 17 out of the 780 modulemodule pairs are significantly overlapped (Fig. 6), showing that our method can find diverse functional modules.
Comparison with sparse SVD on the CGP gene expression data sets
Since L0SVD have shown good performances in simulation study compared to other sparse SVD methods, we compared it with our method to further illustrate the importance of integrating the PPI network. To this end, we also identified 40 gene modules on the CGP data by using L0SVD and Fig. 7 shows the comparing results. We found that the interaction enrichment scores of the identified modules by SNSVD were significantly higher than that by L0SVD (onesided Wilcoxon signed rank test pvalue <0.01) (Fig. 7A). These results demonstrate that SNSVD can find more tightly connected genes than L0SVD by integrating the PPI network. Furthermore, SNSVD obtains a greater number of significant GO BP terms at different levels than L0SVD (onesided Wilcoxon signed rank test pvalue <0.001) (Fig. 7B), showing that incorporating the PPI network does help SNSVD to discover more biological interpretable modules.
Application to the BRCA data sets
Data and preprocessing
We downloaded the processed RNAseq and miRNAseq data of Breast invasive carcinoma (BRCA) from TCGA database [30] (Broad GDAC Firehose: http://firebrowse.org/). We firstly filtered out the genes and the miRNAs which are not expressed in more than 70% samples and the raw gene/miRNA expression values were log2transformed. Secondly, we used the wilcoxon rank sum test to identify differentially expressed genes/miRNAs with bonferroni adjusted pvalue <0.05 between cancer and adjacent normal samples. It causes 9896 differentially expressed genes and 320 differentially expressed miRNAs to be preserved. Thirdly, we imputed the missing values of miRNA and gene expression data by using knearest neighbors [31]. Finally, we extracted the matched gene and miRNA expression matrices across cancer and adjacent normal samples, where A_{1} and B_{1} represent gene and miRNA expression data of cancer samples, respectively and (ii) A_{2} and B_{2} represent gene and miRNA expression data of adjacent normal samples, respectively (Fig. 8A). There are 9896 genes and 320 miRNAs, 760 cancer samples and 87 adjacent normal samples in the BRCA data sets.
Additionally, we also downloaded a PPI network from PathwayCommons database [26], and collected a set of cancer genes from the allOnco database (http://www.bushmanlab.org/links/genelists) which merges some different cancer genes from several databases, and a set of cancer miRNAs from the reference [32].
Identifying differentially coexpressed miRNAgene modules
Recent research revealed that some abnormal miRNAgene regulatory relationship plays key roles in tumor progression and development [33–35]. Some computational methods have been proposed for identifying miRNAgene coexpressed modules by using matched miRNA and mRNA expression data of cancer [13, 36–40]. Though power, these methods do not ensure that the miRNAs and genes in a module are differentially expressed between two biological conditions. Besides, some methods have already been developed for differential coexpression analysis [41–44]. However, these methods only focus on single gene expression data analysis. To this end, we proposed a new framework based on SNSVD for analyzing matched miRNA and mRNA expression data between two biological conditions to identify differentially coexpressed miRNAgene modules (Fig. 8).
Herein, we applied SNSVD to the BRCA data and empirically set λ,k_{v} in SNSVD to yield top ten differentially coexpressed modules for each σ. Each identified miRNAgene module contain about 10 miRNAs and 100 genes. Formally, a miRNAgene module contains a miRNA set and a gene set. We found that as σ becomes larger, the modules identified by SNSVD contain more edges (Table 2). The results showed SNSVD could overcome the drawbacks of sparse SVD (SNSVD with σ=0 in Table 2) to capture the modules with more edges by incorporating the PPI network.
Functional analysis of modules
Without loss of generality, the ten modules identified by SNSVD with σ=60 (See Table 2) were considered for further biological analysis. We found that (i) the average adPCC (absolute differential Pearson Correlation Coefficient) for the identified modules by SNSVD on the BRCA data is larger than the average of all absolute elements of X (Wilcoxon ranksum test, p<1e−16) (Fig. 9A); (ii) more than half of the miRNAs in the 70% (7 of 10) modules are cancer miRNAs, and 80% (8 of 10) modules are significantly enriched at least one KEGG/GO BP pathway (BenjaminiHochberg adjusted p<0.05); (iii) three modules (module 2, 3 and 8) contain significantly more cancer genes with hypergeometric test, p<0.05. More results are shown in Fig. 9B. Additionally, we obtained 39 miRNAs and 961 genes by combining the identified ten modules. We found that about 50% (19 of 39) miRNAs are cancer miRNAs, and about 21% (203 of 961) genes are cancer genes (hypergeometric test, p<4.3e−6).
Discussion
In our previous work, SSVD has been developed for module discovery and its effectiveness has been demonstrated [13]. However, it cannot integrate the gene network data from PPI network. To this end, we develop the SNSVD method that integrates gene expression data and a gene interaction network to identify underlying gene functional modules. In the SNSVD, we define a sparse network regularized function which is a combination of L_{1}regularized norm and networkregularized norm to make the biclustering process tend to select interacted genes in the prior gene interaction network. Experimental results on the CGP and BRCA data demonstrate that SNSVD can overcome the drawbacks of SSVD. Although SNSVD is an effective method, some further studies are deserved to investigate: (1) extend SNSVD to identify nonlinear relationships; (2) extend SNSVD to integrate other omics data, such as DNA methylation data; (3) apply SNSVD to other biological problems.
Conclusions
In this paper, we presented a Sparse Networkregularized SVD (SNSVD) model for networkbased cancer genomics data integration analysis and developed an alternating iterative algorithm to solve the model. By comparing with other representative methods on the simulated data and the real data, we found that SNSVD could find modules with high qualities by integrating the PPI interaction network. By investigating the modules identified by SNSVD on the CGP data, we found that all the genes within the same modules are coexpressed, and most genes in the same modules are connected with each other in the prior PPI network and enriched in at least one gene functional term. Besides, we also applied our method to the BRCA data from TCGA database for identifying ten differentially coexpressed miRNAgene modules. Some breast cancer related miRNAgene modules were discovered. To sum up, our work provides a promising way to integrate the network information into the sparse SVD framework, which can help to find biologically significant functional modules and makes the results easily interpreted. An R package of SNSVD is available at https://github.com/wenwenmin/SNSVD.
Methods
Sparse networkregularized SVD (SNSVD) model
Let \(\boldsymbol {X}\in \mathbb {R}^{p\times n}\) (p genes and n samples) be the gene expression data. Suppose \(\boldsymbol {A} \in \mathbb {R}^{p\times p}\) is an adjacency matrix of a PPI network, where A_{ij}=1 if vertex i and j is connected and A_{ij}=0 otherwise. Thus, the normalized Laplacian matrix L=(L_{ij})_{p×p} encoding the PPI network can be defined as:
where \(d_{i} = \sum _{j=1}^{p} \boldsymbol {A}_{{ij}}\). Correspondingly, we have \(\boldsymbol {u}^{T}\boldsymbol {L}\boldsymbol {u} = \frac {1}{2}\sum _{i}\sum _{j}\boldsymbol {A}_{{ij}}\left (\frac {u_{i}}{\sqrt {d_{i}}}\frac {u_{j}}{\sqrt {d_{j}}}\right)^{2}\), which encourages the estimated coefficients of u to be smooth over adjacent genes in the PPI network A [20]. To further force u to be sparse, we introduce a sparse networkregularized penalty:
where λ and σ are two parameters. In the penalty (3), the L_{1} norm (∥u∥_{1}) is to induce sparsity in u; and the quadratic Laplacian norm (u^{T}Lu) makes the selected genes tend to connect with each other in the PPI network.
To integrate the network information in SVD framework, we present a sparse networkregularized SVD (SNSVD) model as follows:
where c_{1} and k_{v} are two parameters to control the number of selected genes and samples separately. As for the samples, we simply use a L_{0}regularized penalty on v sample variables (corresponding to sample variables) to induce sparseness. Compared to L_{1}norm, L_{0}norm is known as the most essential sparsity measure and has nice theoretical properties [15, 45].
SNSVD algorithm
Since \(\boldsymbol {X}d\boldsymbol {u}\boldsymbol {v}^{T}_{F}^{2} = tr\left (\boldsymbol {X}\boldsymbol {X}^{T}\right)+d^{2}tr\left (\boldsymbol {u}\boldsymbol {v}^{T}\boldsymbol {v}\boldsymbol {u}^{T}\right)2d\boldsymbol {u}^{T}\boldsymbol {X}\boldsymbol {v}\), where tr(·) denotes the trace of a matrix; Both u and v are guaranteed to be two unit vectors, tr(uv^{T}vu^{T})=tr(u^{T}uv^{T}v)=1. Minimizing \(\boldsymbol {X}d\boldsymbol {u}\boldsymbol {v}^{T}_{F}^{2}\) in Eq. (4) is equivalent to minimizing −u^{T}Xv. Although there are three parameters u,v and d to be optimized in Eq. (4). It is notable that once u and v are fixed, then d can be determined d=u^{T}Xv in Eq. (4). Thus, to solve Eq. (4), we just need to optimize u and v. Inspired by Ref. [46], we present an alternating iterative strategy to solve u and v, i.e., fixing v to update u and fixing u to update v.
Fixing v in Eq. (4), it is equivalent to solve the following subproblem:
Let z=Xv, the optimization problem in (5) can be redefined as follows:
To solve it, we write its Lagrangian form as follows:
where λ≥0, η≥0,σ≥0. In order to facilitate the calculation without loss of generality, we use \(\frac {1}{2} \eta \) instead of \(\eta, \frac {1}{2} \sigma \) instead of σ, then Eq. (7) can be rewritten as:
It is a convex function with respect to u, therefore its optimal solution can be characterized by some subgradient equations (see e.g., [47]). Since L=I−D^{−1/2}AD^{−1/2} (based on Eq. 2). For convenience, let W=I−L=D^{−1/2}AD^{−1/2} (D is a diagonal matrix and \(\boldsymbol {D}_{{ii}} = \sum _{j} \boldsymbol {A}_{{ij}}\)), then we have the subgradient equations of Eq. (8) as:
where s_{j}=sign(u_{j}) if u_{j}≠0 and s_{j}∈{t,t≤1} otherwise; and W_{j} is the jth row of matrix W. Let the solution of (8) be \( \widehat {\boldsymbol {u}}=\left (\widehat {\boldsymbol {u}}_{1}, \widehat {\boldsymbol {u}}_{2}, \cdots, \widehat {\boldsymbol {u}}_{p}\right)\). By using a coordinate descent method [48, 49], we obtain the following coordinate update rule for \(\widehat {\boldsymbol {u}}_{j}\):
Define \(\mathcal {S}(a,\lambda) = \text {sign}(a)(a\lambda)_{+}\), we have \(\widehat {\boldsymbol {u}}_{j} = \mathcal {S}(\boldsymbol {z}_{j} + \sigma \boldsymbol {W}_{j} \widehat {\boldsymbol {u}}, \lambda)/(\eta +\sigma)\). Let \(\breve {\boldsymbol {u}}_{j} = \mathcal (\boldsymbol {z}_{j} + \sigma \boldsymbol {W}_{j}\widehat {\boldsymbol {u}}, \lambda)\) and \(\breve {\boldsymbol {u}}=\left (\breve {\boldsymbol {u}}_{1},...,\breve {\boldsymbol {u}}_{p}\right)^{T}\), we can obtain a normalized solution \(\boldsymbol {u}=\frac {\widehat {\boldsymbol {u}}}{\ \widehat {\boldsymbol {u}} \_{2}}=\frac {\breve {\boldsymbol {u}}}{\ \breve {\boldsymbol {u}} \_{2}}\). In a word, we use a coordinate descent method to minimize Eq. (8) and update one u_{j} at a time while keeping u_{k} fixed for all k≠j.
Fixing u in Eq. (4), it is equivalent to solve the following subproblem:
Let \(\boldsymbol {z}_{v} = \boldsymbol {X}^{T}\boldsymbol {u}, \boldsymbol {\widehat {v}} = d\boldsymbol {v}\), we thus have \(\\boldsymbol {X}  d\boldsymbol {uv}^{T}\_{F}^{2} = \\boldsymbol {z}_{v}\boldsymbol {\widehat {v}}\_{2}^{2} + c \), where c=tr(X^{T}X)−u^{T}XX^{T}u. Obviously c is a constant value with respect to v. Thus problem (11) is equivalent to:
Its optimal solution is \(\boldsymbol {\widehat {v}} = \boldsymbol {z}_{v} \bullet I(\boldsymbol {z}_{v} \geq \boldsymbol {z}_{v}_{(k_{v})})\) where I(·) is the indicator function and "∙" is point multiplication function, and z_{v}_{(i)} denotes the ith order statistic of z_{v}, i.e. z_{v}_{(1)}≥z_{v}_{(2)}≥,...,≥z_{v}_{(n)}. In other words, we only keep the k_{v} variables of z_{v} corresponding to its k_{v} largest absolute values. The normalized optimal solution of Eq. (11) is \(\boldsymbol {v} = \widehat {\boldsymbol {v}}/\\widehat {\boldsymbol {v}}\_{2}\), i.e.,
Finally, we develop an alternating iterative algorithm by alternately updating u and v to solve SNSVD model. The details of this algorithm is given in Algorithm 1, and its time complexity is \(\mathcal {O}\left (Tnp + Tp^{2} + Tn^{2}\right)\), where T is the number of iterations.
Convergence analysis of SNSVD algorithm
Next, we give the convergence analysis of Algorithm 1. In fact, Algorithm 1 is to solve the Lagrangian form of problem (4) as follows:
Let H(u,v)=−u^{T}Xv+σu^{T}Lu,f(u)=ρ(u)+λ∥u∥_{1} and g(v)=ρ(v)+τ(v,k_{v}) where
Therefore the Lagrangian form of problem (4) can be written as F(u,v)=H(u,v)+f(u)+g(v) which is a semialgebraic function [46]. Based on the Theorem 1 in [46], Algorithm 1 converges to a critical point of F(u,v).
Parameter selection of SNSVD algorithm
As to λ and k_{v}’s choice in Algorithm 1 when it is applied to the CGP gene expression data, we select a suitable λ to force the estimated u only containing 200 nonzero elements which is beneficial for further analysis of the biological function of the module and set k_{v}=50 (control the sample sparsity) which ensures that the number of samples within in the module is approximately the same as the number of samples of a subtype. As to σ’s choice in Algorithm 1, we present a 5fold crossvalidation framework (Algorithm 2). To this end, we define a binary matrix is.na(X) with the same size of X and the elements are 1 if they are missing in X, 0 otherwise.
Learning multiple pairs of singular vectors using SNSVD
It is notable that every run of Algorithm 1 can only obtain a pair of sparse singular vectors u and v (Fig. 1). In order to identify multiple modules, we can repeat running Algorithm 1. After each turn of the iteration, we use the obtained u,v and d to modify the gene expression data X, (X:=X−duv^{T}), the modified X is then used as the new input data for the next run to obtain the next pair of singular vectors. Moreover, we notice that Algorithm 1 may get different local optima with different initials, we run Algorithm 1 five times with different initials which are generated according to the multivariate standard normal distribution and choose the best one as the final solution of each turn. The detailed procedure is described in Algorithm 3.
Modularity score
To assess whether the genes within the same module are coexpressed/correlated, we use a modularity score to describe the overall coexpression of genes within the module. For a given module k containing p_{k} genes and n_{k} samples, we first calculate the correlation between gene i and j across the n_{k} samples, denoted as w_{ij}. For convenience, we force to set w_{ii}=0 for each i. Then the modularity score of the module can be defined as:
Intuitively, if a module has a high modularity score, then the genes within the module is highly coexpressed.
Genegene interaction enrichment score
In order to evaluate whether the genes within the same module are tightly connected in the prior PPI network, we use the right tailed hypergeometric test to compute a significance level of each module. Suppose that the PPI network contains n genes and m edges, and module i contains n_{i} genes and m_{i} edges, then the significance level of module k can be calculated via the following equation:
where \(N = \binom {n}{2}\) and \(N_{i} = \binom {n_{i}}{2}\). Accordingly, we can define the genegene interaction enrichment score s(i) of the module i by the following formula:
The higher the genegene interaction enrichment score is, the denser the genes connect with each other. If the score is higher than 1.3, then the genes are significantly interconnected with each other in the PPI network.
Availability of data and materials
The CGP datasets used during the current study are available in the Cancer Genome Project (CGP) repository (http://www.cancerrxgene.org/downloads). And the BRCA datasets used during the current study are available in the TCGA repository (Broad GDAC Firehose: http://firebrowse.org/).
Abbreviations
 SVD:

Singular value decomposition
 SSVD:

Sparse singular value decomposition
 SNSVD:

Sparse networkregularized singular value decomposition
 LASSO:

Least absolute shrinkage and selection operator
 SNR:

Signaltonoise ratio
 CGP:

Cancer genome project
 TCGA:

The cancer genome atlas
 PPI:

Protein protein interaction
 KEGG:

Kyoto encyclopedia of genes and genomes
 GO:

Gene ontology
References
Wan Q, Dingerdissen H, Fan Y, Gulzar N, Pan Y, Wu TJ, Yan C, Zhang H, Mazumder R. Bioxpress: An integrated RNAseqderived gene expression database for pancancer analysis. Database (Oxford). 2015; 2015:1–13.
Vallejos CA, Risso D, Scialdone A, Dudoit S, Marioni JC. Normalizing singlecell RNA sequencing data: challenges and opportunities. Nat Methods. 2017; 14(6):565–71.
Iorio F, Knijnenburg TA, Vis DJ, Bignell GR, Menden MP, Schubert M, Aben N, Gonçalves E, Barthorpe S, Lightfoot H, et al. A landscape of pharmacogenomic interactions in cancer. Cell. 2016; 166(3):740–54.
Lee M, Shen H, Huang JZ, Marron J. Biclustering via sparse singular value decomposition. Biometrics. 2010; 66(4):1087–95.
Liquet B, de Micheaux PL, Hejblum BP, Thiébaut R. Group and sparse group partial least square approaches applied in genomics context. Bioinformatics. 2015; 32(1):35–42.
Min W, Liu J, Zhang S. Edgegroup sparse pca for networkguided high dimensional data analysis. Bioinformatics. 2018; 34(20):3479–87.
Liu X, Chang X, Liu R, Yu X, Chen L, Aihara K. Quantifying critical states of complex diseases using singlesample dynamic network biomarkers. PLoS Comput Biol. 2017; 13(7):1005633.
Yu X, Zhang J, Sun S, Zhou X, Zeng T, Chen L. Individualspecific edgenetwork analysis for disease prediction. Nucleic Acids Res. 2017; 45(20):170.
Eren K, Deveci M, Küçüktunç O, Ümit V. Çatalyürek: A comparative analysis of biclustering algorithms for gene expression data. Brief Bioinforma. 2013; 14(3):279–92.
Sill M, Kaiser S, Benner A, KoppSchneider A. Robust biclustering by sparse singular value decomposition incorporating stability selection. Bioinformatics. 2011; 27(15):2089–97.
Oghabian A, Kilpinen S, Hautaniemi S, Czeizler E. Biclustering methods: Biological relevance and application in gene expression analysis. PLoS ONE. 2014; 9(3).
Chen S, Liu J, Zeng T. Measuring the quality of linear patterns in biclusters. Methods. 2015; 83:18–27.
Min W, Liu J, Luo F, Zhang S. A twostage method to identify joint modules from matched microRNA and mRNA expression data. IEEE Trans Nanobiosci. 2016; 15(4):362–370.
Yang D, Ma Z, Buja A. Rate optimal denoising of simultaneously sparse and low rank matrices. J Mach Learn Res. 2016; 17(1):3163–89.
Asteris M, Kyrillidis A, Koyejo O, Poldrack R. A simple and provable algorithm for sparse diagonal CCA. In: International Conference on Machine Learning: 2016. p. 1148–1157.
Sokolov A, Carlin DE, Paull EO, Baertsch R, Stuart JM. Pathwaybased genomics prediction using generalized elastic net. PLoS Comput Biol. 2016; 12(3):e1004790.
Hill SM, Heiser LM, et al. Inferring causal molecular networks: empirical assessment through a communitybased effort. Nat Methods. 2016; 13(4):310–8.
Enrico G. Using prior knowledge from cellular pathways and molecular networks for diagnostic specimen classification. Brief Bioinforma. 2016; 17(3):440–52.
Lee E, Chuang HY, Kim JW, Ideker T, Lee D. Inferring pathway activity toward precise disease classification. PLoS Comput Biol. 2008; 4(11):e1000217.
Li C, Li H. Networkconstrained regularization and variable selection for analysis of genomic data. Bioinformatics. 2008; 24(9):1175–82.
Sun H, Feng R, Lin W, Li H. Networkregularized highdimensional cox regression for analysis of genomic data. Stat Sin. 2013; 24(3):1433–59.
Chen J, Zhang S. Integrative analysis for identifying joint modular patterns of geneexpression and drugresponse data. Bioinformatics. 2016; 32(11):1724–32.
Zhu F, Liu J, Min W. Gene functional module discovery via integrating gene expression and ppi network data. In: International Conference on Intelligent Computing: 2019. p. 116–126. https://doi.org/10.1007/9783030269692_11.
Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Ser B (Methodological). 1996; 58(1):267–88.
Fan J, Li R. Variable selection via nonconcave penalized likelihood and its oracle properties. J Am Stat Assoc. 2001; 96(456):1348–60.
Cerami EG, Gross BE, et al. Pathway commons, a web resource for biological pathway data. Nucleic Acids Res. 2011; 39(Database Issue):685–90.
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protocol. 2009; 4(1):44–57.
Leeksma OC, de Miranda NF, Veelken H. Germline mutations predisposing to diffuse large Bcell lymphoma. Blood Cancer J. 2017; 7(2):532.
Disis ML. Immune regulation of cancer. J Clin Oncol. 2010; 28(29):4531–8.
Lander ES, Park PJ. The cancer genome atlas pancancer analysis project. Nat Genet. 2013; 45(10):1113–20.
Troyanskaya O, Cantor M, et al. Missing value estimation methods for DNA microarrays. Bioinformatics. 2001; 17(6):520–5.
Xie B, Ding Q, Han H, Wu D. miRCancer: a microRNA cancer association database constructed by text mining on literature. Bioinformatics. 2013; 29(5):638–44.
Garzon R, Calin GA, Croce CM. MicroRNAs in cancer. Ann Rev Med. 2009; 60:167–79.
Adams BD, Kasinski AL, Slack FJ. Aberrant regulation and function of microRNAs in cancer. Curr Biol. 2014; 24(16):762–76.
Iorio MV, Croce CM. MicroRNAs in cancer: small molecules with a huge impact. J Clin Oncol. 2009; 27(34):5848.
Zhang S, Li Q, Liu J, Zhou XJ. A novel computational framework for simultaneous integration of multiple types of genomic data to identify microRNAgene regulatory modules. Bioinformatics. 2011; 27(13):401–9.
Zhang S, Liu CC, Li W, Shen H, Laird PW, Zhou XJ. Discovery of multidimensional modules by integrative analysis of cancer genomic data. Nucleic Acids Res. 2012; 40(19):9379–91.
Bryan K, et al. Discovery and visualization of miRNAmRNA functional modules within integrated data using bicluster analysis. Nucleic Acids Res. 2013; 42(3):17.
Li Y, Liang C, Wong KC, Luo J, Zhang Z. Mirsynergy: detecting synergistic miRNA regulatory modules by overlapping neighbourhood expansion. Bioinformatics. 2014; 30(18):2627–35.
Jin D, Lee H. A computational approach to identifying genemicroRNA modules in cancer. PLoS Comput Biol. 2015; 11(1):1004042.
Tesson BM, Breitling R, Jansen RC. DiffCoEx: a simple and sensitive method to find differentially coexpressed gene modules. BMC Bioinformatics. 2010; 11:497.
Ideker T, Krogan NJ. Differential network biology. Mol Syst Biol. 2012; 8:565.
Ha MJ, Baladandayuthapani V, Do KA. Dingo: differential network analysis in genomics. Bioinformatics. 2015; 31(21):3413–20.
Zhu L, et al. MetaDCN: metaanalysis framework for differential coexpression network detection with an application in breast cancer. Bioinformatics. 2016; 33(8):1121–9.
Yang F, Shen Y, Liu ZS. The proximal alternating iterative hard thresholding method for L_{0} minimization, with complexity \(\mathcal {O}(1/\sqrt {k})\). J Comput Appl Math. 2017; 311:115–29.
Bolte J, Sabach S, Teboulle M. Proximal alternating linearized minimization for nonconvex and nonsmooth problems. Math Program. 2014; 146(12):459–94.
Nesterov Y. Primaldual subgradient methods for convex problems. Math Program. 2009; 120(1):221–259.
Friedman J, Hastie T, Höfling H, Tibshirani R, et al. Pathwise coordinate optimization. Ann Appl Stat. 2007; 1(2):302–332.
Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010; 33(1):1–22.
Acknowledgements
Wenwen Min would like to thank the support of National Center for Mathematics and Interdisciplinary Sciences, Academy of Mathematics and Systems Science, CAS during his visit.
About this supplement
This article has been published as part of BMC Genomic Data Volume 22 Supplement 1, 2021: Proceedings of the 2019 International Conference on Intelligent Computing (ICIC 2019): genomic data. The full contents of the supplement are available online at https://bmcgenomdata.biomedcentral.com/articles/supplements/volume22supplement1.
Funding
This work was supported by the National Science Foundation of China [No. 61802157] and Natural Science Foundation of Jiangxi Province of China [No. 20192BAB217004]. The funding bodies had no role in the design of the study or collection, analysis and interpretation of the data.
Author information
Authors and Affiliations
Contributions
The authors wish it to be known that, in their opinion, F.Z. and W.M. should be regarded as Joint First Authors. F.Z., W.M., J.Liu and J.Li developed the methodology and executed experiments. F.Z. and W.M. drafted the manuscript, J.Liu and J.Li revised the manuscript. W.M., J.Liu and J.Li supervised this study. The authors read and approved the final manuscript.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Zhu, F., Li, J., Liu, J. et al. Networkbased cancer genomic data integration for pattern discovery. BMC Genom Data 22 (Suppl 1), 54 (2021). https://doi.org/10.1186/s1286302101004y
Published:
DOI: https://doi.org/10.1186/s1286302101004y
Keywords
 Gene coexpression analysis
 Differentially coexpression analysis
 Gene interaction network
 Sparse SVD
 Structured sparse learning