Algorithm to Identify Frequent Coupled Modules from Two-Layered Network Series: Application to Study Transcription and Splicing Coupling

Lemma 1

The problem of maximizing \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$O ({\bf x}, {\bf y}, {\bf w}, Z) = H_{\cal G} ({\bf x}, {\bf w}) + \lambda H_{\cal E} ({\bf y}, {\bf w}) + \mu O_R ({\bf x}, {\bf y}, Z)$$ \end{document} with the constraints \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\sum\nolimits_{i = 1}^N x_i = K_1, \sum\nolimits_{j = 1}^M y_j = K_2, \sum\nolimits_{k = 1}^L w_k = K_3$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\sum\nolimits_{j = 1}^M z_{ij} \le K_2$$ \end{document} for all \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$i = 1, \ldots, N$$ \end{document} is NP-hard.

To address these two drawbacks, we instead solved a continuous optimization problem with the same objective by relaxing integer constraints to continuous constraints. That is, we looked for non-negative real vectors x, y, w and the non-negative real sparse matrix Z that jointly maximize O(x, y, w, Z). This optimization problem is formally expressed as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & \max \nolimits_{{\bf x}, {\bf y}, {\bf w}, Z \in {\mathbb R}^+} \quad O ({\bf x}, {\bf y}, {\bf w}, Z) = H_{\cal G} ({\bf x}, {\bf w}) + \lambda H_{\cal E} ({\bf y}, {\bf w}) + \mu O_R ({\bf x}, {\bf y}, Z) \\ & \hbox{subject to} \begin{cases} \hbox{Constraint I:} \ \ \quad \ \parallel {\bf x} \ \parallel_f = 1, \qquad 1 \le f < 2 \\ \hbox{Constraint II:} \quad \ \ \parallel {\bf y} \ \parallel_g = 1, \qquad 1 \le g < 2 \\ \hbox{Constraint III:} \quad \ \parallel {\bf w} \ \parallel_h = 1, \qquad h \ge 2 \\ \hbox{Constraint IV:} \quad \sum_{j = 1}^Mz_{ij} \le 1, \ \ \ i = 1, 2, \ldots, N\end{cases} & (1) \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\mathbb R}^{+}$$ \end{document} is a non-negative real space of vectors or matrix, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\ \parallel {\bf x} \ \parallel_p = \left(\sum \nolimits_{i} \mid x_i \mid^p \right)^{1 / p}$$ \end{document} is the L_p norm of the vector x. The L_p = 1 (1 ≤ p < 2) norm constraint of genes and exons (Constraints I and II above) is a well-known convex sparse coding scheme (Li et al., 2011; Ding et al., 2006) that can lead to the sparse solution where only a few elements of the gene (exon) membership vector x (y) are significantly different from zero and therefore their corresponding genes (exons) can be selected into the module. Therefore it is well suited to our purpose of selecting only a small number of genes (exons) into the module. On the contrary, the L_h = 1 (h ≥ 2) norm constraints of datasets (Constraint III) gives a “smooth” solution where the elements of the optimized vector w are approximately equal and therefore lead to the discovered coupled cluster occurring in as many datasets as possible. This problem definition and formulation can be viewed as the generalization of the recurrent heavy subgraph (RHS) problem defined in Li et al. (2011). However, requiring mapping relationships between gene subgraph and exon subgraph by introducing a relation criterion to the problem makes it much more complicated than the RHS problem. Such nature of simultaneously analyzing two different tensors makes it distinct from the problem defined in the single tensor.

Eq. (1) defines a tensor-based optimization formulation for the problem of identifying FCCs. By solving Eq. (1), users can easily identify the top-ranking datasets (after sorting elements of w in non-increasing order) and top-ranking genes/exons (after sorting elements x/y in non-increasing order) contributing to the objective function. After rearranging the tensors in this manner, the optimum FCC occupies the corner of the 3D tensor \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal G}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal E}$$ \end{document} . We then mask the edges contained in this cluster in both gene and exon networks with zeros and optimize Eq. (1) again to search for the next module. The constant weights λ and μ in the objective function provide us the opportunity to control the importance of criteria and data. To fully exploit such capability, we collect all FCCs discovered in all combinations ⁴ of different λ and μ, then remove duplicates. Details of the optimization procedure and pattern extraction process will be presented in the next section. Therefore, our problem formulation requires only three parameters f, g and h which can be fixed through simulation study.

We performed simulation studies to determine suitable values for the parameters f, g, and h by applying our tensor method to collections of random weighted networks. In subsets of these networks, we randomly placed FCC of varying size, occurrence, and heaviness. We then tested different combinations of f, g, and h, and adopted the combination (f = g = 1, h = 10) that led to the discovery of the most FCCs. More details on these simulations are provided in Appendix D.

2.2. Optimization algorithm and pattern extraction

We derived an iterative algorithm to maximize the objective function with the constraints defined in Eq. (1). This algorithm repeatedly applies the four updating steps (shown below) until the objective function converges to a fixed point, yielding the solution to Eq. (1). Each step of this procedure individually updates x (y or w or Z) while fixing the other three variables.

The following lemmas guarantee the correctness and convergence of this algorithm (see Appendices E and F for their proofs).

Lemma 2 (Correctness)

The fixed point of the updating steps in Algorithm 1 satisfies the Karush-Kuhn-Tucker (KKT) condition of Eq.(1).

Lemma 3 (Convergence)

The objective function O(x, y, w, Z) is nondecreasing under the updating rules in Algorithm 1. The objective function is invariant under these updates if and only if x, y, w, and Z are at a stationary point of the objective function.

As the relation matrix R is extremely sparse (with only M non-zeros), all element-wise matrix multiplication operations involving R and Z in Algorithm 1 could be implemented efficiently by only treating z_ij whose corresponding r_ij ≠ 0 and simply assigning the rest z_ij = 0 whose corresponding r_ij = 0. Thus, the complexity of Algorithm 1 is linear with respect to the total number of edges in two tensors \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal G}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal E}$$ \end{document} .

According to the definition of three aforementioned membership vectors x, y, and w, FCCs can be intuitively obtained by including those top-ranking genes, exons, and datasets with large membership values (after sorting x, y and w in descending order, respectively). In this article, we empirically define an FCC extracted from these top-ranking genes, exons and datasets as a triplet (G′, E′, D′), i.e., a set of ≥ 5 genes G′, a set of ≥ 5 exons E′ and a set of ≥ 2 active datasets D′ satisfying: (1) each gene subgraph formed by G′ and exon subgraph formed by E′ in each active dataset of D′ has an average edge weight (or so-called “heaviness”) ≥ a predefined threshold (i.e., ≥ 0.4) and (2) the fraction of host genes of E′ that are included in G′, so-called \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$coverage_{gene} (G^{\prime}, E^{\prime}) = \frac {\mid G^{\prime} \cap hostgene (E^{\prime}) \mid} {\mid G^{\prime} \mid}$$ \end{document} ≥ a predefined threshold (i.e., ≥ 0.7) and the fraction of exons of G′ which are included in E′, so-called \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$coverage_{exon} (G^{\prime}, E^{\prime}) = \frac {\mid E^{\prime} \cap exon (G^{\prime}) \mid} {\mid E^{\prime} \mid}$$ \end{document} ≥ a predefined threshold (i.e., ≥ 0.7). The first criterion gives a certain degree of how heavy and frequent the cluster should be, and the second criterion gives a measure of what degree of “coupling” should exist between a gene set and an exon set. Two examples of coupled gene set and exon set with coverage_gene > 0.6 and coverage_exon > 0.6 are shown in Figure 3. Based on this definition, we can generate one or more FCCs that are highly overlapping with each other from combinations of different sizes of top-ranking genes, exons, and datasets which satisfy above two criteria. Therefore, we call a group of FCCs derived from the same membership vectors x, y, and w as a family of FCCs.

OPEN IN VIEWER

FIG. 3.

Illustration of the coupling degree between a gene set and an exon set. Blue edges between genes and exons indicate that there are relationships between them. Red edges within genes or exons indicate that they are co-expressed or co-spliced. (A) Perfect coupling. (B) Coupling with coverage_gene > 0.6 and coverage_exon > 0.6.

3.1. Frequent coupled clusters are likely to represent functional modules and protein complexes

We evaluated the functional homogeneity of the merged genes and host genes of exons in a frequent coupling clusters using specific Gene Ontology (GO) terms (Ashburner et al., 2000). To ensure the specificity of GO terms, we filtered out those general terms associated with >300 genes. If the gene set of any FCCs within a family is significantly enriched in a GO term with a q-value <0.05 (the q-value is the hypergeometric p-value after a False Discovery Rate multiple testing correction), we declare this family as functionally homogeneous. According to this definition, 50.4% of the FCC families was functionally homogenous. In an ensemble of randomly generated FCC families with the same size distribution as our patterns, only 5.7% of the families was functionally homogenous. This enrichment fold ratio of 8.8 between real and random patterns demonstrates the strong biological relevance of the identified patterns. Furthermore, the enrichment fold ratio increases dramatically with increasing heaviness and recurrence of FCCs (Fig. 4). For example, when the FCCs are required to recur in at least 4 datasets, their enrichment fold ratio compared to random patterns increases to 48.0, highlighting the benefit of pursuing integrative analysis of multiple RNA-seq datasets. When the FCCs are required to be heavier than 0.6, their enrichment fold ratio raises to 65.0.

OPEN IN VIEWER

FIG. 4.

FCCs are highly likely to represent functional modules. The enrichment fold ratio, calculated by dividing the percentage of functionally homogenous FCC families by the percentage of functionally homogenous random patterns, increases significantly with the recurrence and heaviness of the FCCs. (A) The minimum number of active datasets in which FCCs recur is varied while fixing the minimum heaviness as 0.4. (B) The minimum heaviness of FCCs is varied while fixing the minimum number of active datasets as 2.

Functional modules represented by FCCs cover a large number of categories related to post-transcriptional processing, such as “regulation of mRNA stability,” “posttranscriptional regulation of gene expression,” “RNA polyadenylation,” “RNA splicing,” “mRNA processing,” “RNA transport,” “establishment of RNA localization,” “mRNA 3′-end processing,” “nucleic acid transport,” and “RNA transport.” This is not surprising, since recent experimental studies revealed the existence of factories in which individual machines for the transcription, post-transcription processing, and mRNA transport are functionally coupled (Reed, 2003). During expression of protein-coding genes, pre-mRNAs are transcribed in the nucleus and undergo several processing steps, including capping, splicing, 3′-end processing and polyadenylation, prior to be transported to the cytoplasm for translation. The coupling between those individual steps is believed to enable the proofreading and streamlining of the entire process of gene expression in vivo (Reed, 2003). Other functional categories over-represented by FCCs include “nucleosome organization,” linking to the recent evidences that the transcription elongation, an important step in transcription-splicing coupling, can be regulated by nucleosome (Lindstrom et al., 2003).

Similar results were achieved by using the Comprehensive Resource of Mammalian protein complexes (CORUM) database (Ruepp et al., 2010) to assess the association between FCCs and known protein complexes. The merged sets of genes and exon-host-genes of 78.6% of FCC families are significantly enriched in genes from the same complexes with a q-value <0.05, compared to only 5.05% of randomly generated patterns with the same size distribution. The enrichment fold ratio dramatically increases with increasing heaviness and recurrence of the FCCs (Fig. 5).

OPEN IN VIEWER

FIG. 5.

FCCs are highly likely to represent protein complexes. The enrichment fold ratio increases significantly with the recurrence and heaviness of the FCCs. (A) The minimum number of active datasets in which FCCs recur is varied while fixing the minimum heaviness as 0.4. (B) The minimum heaviness of FCCs is varied while fixing the minimum number of active datasets as 2.

Protein complexes represented by FCCs generally fall into the following categories: (1) Splicing regulation, e.g. SNW1 complex, CDC5L complex, C complex spliceosome; (2) Transcription regulation, e.g., SNW1 complex, TLE1 complex, H2AX complex; (3) microRNA processing, e.g., Large Drosha complex, DGCR8 multiprotein complex; (4) DNA or chromatin processing, e.g., DNA-PK-Ku-eIF2-NF90-NF45 complex, Emerin complex; (5) Translation, e.g., eif complex, Nop56p-associated pre-rRNA complex; and (6) Protein folding, e.g., CCT micro-complex. Interestingly, the functions of all of those complexes are exclusively centered on the central-dogma-process, namely, the sequential information transfer from DNAs to proteins. Likely, all functions under coupled transcription and splicing regulation are those that need strict and precise regulation to ensure efficient and error-free information processing. Some of those complexes have been experimentally confirmed to couple splicing and transcription. For example, the SNW1 complex contains components of the spliceosome, and it is known to couple vitamin D receptor-mediated transcription and RNA splicing (Zhang et al., 2003). Since research on such coupling is still in its infancy, we expect that more protein complexes will be shown to link those individual regulatory steps into a coherent and smooth production line.

3.2. Frequent coupled clusters are likely to represent transcription and splicing modules

Because the gene set of an FCC are strongly co-expressed in multiple datasets generated under different conditions, it is likely to represent a transcription module. To assess this possibility, we used the 191 ChIP-seq profiles generated by the Encyclopedia of DNA Elements (ENCODE) consortium (Thomas et al., 2007). This dataset includes the genome-wide binding of 40 transcription factors (TF), 9 histone modification marks, and 3 other marks (DNase, FAIRE, and DNA methylation) on 25 different cell lines. For a detailed description of the signal extraction procedure, see Appendix H. These data provide the potential targets of regulatory factors that may or may not be active under a specific condition. However, if the gene set of an FCC is found to be highly enriched in the targets for a regulatory factor, then this factor is likely to actively regulate genes in this FCC, and we declare the family containing this FCC to be “transcriptionally homogenous.” According to this definition, 61.2% of FCC families are transcriptionally homogenous with an enrichment q-value <0.05, compared to 8.9% of the randomly generated patterns with the same size distribution. As shown in Figure 6, with the increasing heaviness and recurrence of FCCs, the enrichment fold ratio between FCCs and random patterns dramatically increases. The five most frequently enriched regulators are YY1, E2F4, c-Myc, MAX, and TAF, all of which have been implicated in cancer pathogenesis or progression (Gordon et al., 2005; Nevins, 2001; Little et al., 1983; Nair and Burley, 2003; Suzuki et al., 2010).

OPEN IN VIEWER

FIG. 6.

FCCs are highly likely to represent transcription modules. The enrichment fold ratio increases significantly with the recurrence and heaviness of the FCCs. (A) The minimum number of active datasets in which FCCs recur is varied while fixing the minimum heaviness as 0.4. (B) The minimum heaviness of FCCs is varied while fixing the minimum number of active datasets as 2.

Since the exon set of an FCC has highly correlated inclusion rate profiles across different experimental conditions, they are likely to be co-regulated by the same splicing factors. We collected the experimental RNA target motifs (2,220 RNA binding sites) of 62 splicing factors from the SpliceAid2 database (Piva et al., 2011). For each exon in an FCC, we retrieved the internal exon region and its 100bp flanking intron region which are enriched in the motifs of those 62 splicing proteins by performing BLAST search (E-score <0.05). If the exon set of any FCC within a family is significantly (q-value <0.05) enriched in the targets of a splicing factor, we consider this family to be “splicing homogenous.” Although the collection of known splicing motifs is very limited, we still observed that 8.9% of the FCC families are splicing homogenous, compared to 5.1% of randomly generated patterns with the same size distribution. Again, as shown in Figure 7, with the increasing heaviness and recurrence of the FCCs, the enrichment fold ratio increases as well. The five most frequently enriched splicing regulators are PSF, hnRNP-D, hnRNP-C1, SLM-2, and HuB.

OPEN IN VIEWER

FIG. 7.

FCCs are likely to represent splicing modules. The enrichment fold ratio increases significantly with the recurrence and heaviness of the FCCs. (A) The minimum number of active datasets in which FCCs recur is varied while fixing the minimum heaviness as 0.4. (B) The minimum heaviness of FCCs is varied while fixing the minimum number of active datasets as 2.

3.3. Exploring the mechanisms of the transcription and splicing coupling

Recently, individual case studies have suggested that the strong connection between transcription and splicing might be due to (a) the coupled recruitment of transcription and splicing factors and (b) the kinetic elongation control by RNA Pol II (Kornblihtt et al., 2004). However, global analyses of these hypotheses have not been conducted. Our catalogues of coupled transcription-splicing modules provide a unique opportunity for us to systematically assess the mechanism (a). In particular, we evaluate the hypothesis that the functionally coupled recruitment of both types of factors could be mediated by the direct or indirect protein-protein interactions between them. This hypothesis is supported by the following biological evidences: (1) the presence of an intron or simply a 5′ splice site immediately downstream from a promoter, i.e., the short distance between the splicing and transcription factors, greatly enhances transcription, both in mammalian and yeast genes (Kornblihtt et al., 2004); and (2) human spliceosome contains at least 30 proteins with known or putative roles involved in transcription, presenting direct interactions between the transcription and splicing (Zhou et al., 2002; Rappsilber et al., 2002).

We used protein-protein interaction (PPI) repository BioGRID (Stark et al., 2006) to examine the PPIs between enriched transcription and splicing factors in coupled module families. We consider not only the direct PPIs between two factors, but also the indirect interactions through a mediator protein (one-hop interaction), because experiments showed that transcription and splicing factors could be recruited by the same other proteins (Nayler et al., 1998; Blencowe et al., 1999). In order to have a broad coverage of transcription factors, we used 109 human transcription factors in ENCODE (Thomas et al., 2007) and JASPAR database (Sandelin et al., 2004), where 10,278 DNA binding motifs are downloaded from Tabach et al. (2007), with which we performed BLAST search (E-score <0.05) in 1000bp regions of the transcription start sites of all genes, and then used hyper-geometry test (q-value <0.05) to identify enriched transcription factors for genes in an FCC. The putative splicing factors for each FCC were identified based on the description in Section 3.2. The total number of PPIs between transcription and splicing factors within the same FCC families (including direct and one-hop interactions) is 105, compared to only 14.8 in random families with the same size of gene and exon sets, with the fold ratio 7.1. This evidence supports that PPI mediated association can be an important mechanism of transcription-splicing coupling.

An interesting example is the FCC family associated with the splicing regulator PSF and the transcription regulator TAT-SF1. PSF is a part of the human spliceosome and essential in RNA splicing (Peng et al., 2006; Garcia-Jurado et al., 2011). TAT-SF1 functions as a general transcription factor playing a role in the process of transcriptional elongation (Zhou and Sharp, 1996). Based on the kinetic model that transcription rate affects splicing outcomes during transcription elongation (Kornblihtt et al., 2004), previous studies suggested that TAT-SF1 couple with PSF for co-transcriptional splicing. In fact, the spliceosome has been shown to be associated with the transcription elongation complexes (Kameoka et al., 2004), and TAT-SF1 was also characterized to be part of the spliceosome (Rappsilber et al., 2002). We found that TAT-SF1 and PSF both interact with the protein WBP4, a general spliceosomal protein that may play a role in cross-intron bridging of U1 and U2 snRNPs (Bedford et al., 1998). Thus, we hypothesize that the elongation factor TAT-SF1 may functionally couple with PSF through the mediator WBP4 in spliceosome.

As another example, the co-spliced exons of an FCC family is enriched in the RNA binding motif of the splicing factor hnRNP D, and the co-expressed genes of the same FCC family is over-represented by the DNA binding motif of the transcription factor TDP43. hnRNP D belongs to heterogeneous nuclear ribonucleoprotein family, which generally suppress splicing through binding to exonic splicing silencer (Matlin et al., 2005). TDP43 was originally identified as a transcriptional repressor that binds to chromosomally integrated TAR DNA and represses HIV-1 transcription (Ou et al., 1995). Interestingly, TDP43 was later shown to also have splicing regulatory functions, wherein the recruitment of TDP43 inhibits exon recognition (Buratti and Baralle, 2001). We suspect that TDP43 may couple with hnRNP D through hnRNP A1, because there is a direct protein-protein interaction between hnRNP A1 and hnRNP D/TDP43 (Stark et al., 2006), whereas TDP43-hnRNP A1 interaction is suggested to be mediated by C-terminal region of TDP43 (Ayala et al., 2010). Furthermore, experiments have shown that TDP43, hnRNP D and hnRNP A1 are all shuttled between nucleus and cytoplasm (Loflin et al., 1999; Ayala et al., 2008). Thus, it is very possible that TDP43 and hnRNP couple the transcription-splicing process via hnRNP A1.

A. RNA-seq datasets selection and processing, and network construction

From NCBI's Sequence Read Archive (SRA), we selected all human RNA-seq datasets, each of which contains at least six samples (the minimum for robust correlation estimation). This results in a total of 38 datasets. For each dataset, we used the Tophat (Trapnell et al., 2009) tool to map short reads to the hg18 reference genome, setting the program options to report only the optimal alignment and discard those reads that map equally well to multiple positions. Next, we applied the transcript assembly tool Cufflinks (Trapnell et al., 2010) to estimate expressions for all transcripts with known UCSC transcription annotations (Fujita et al., 2011). We calculated the inclusion rate of each exon in every sample, as the ratio between its expression (i.e., sum of FPKM over all transcripts that cover the exon) and the expression of the host gene (i.e., sum of FPKM over all transcripts of the gene). It is worth noting that in RNA-seq experiments, a gene expression with low FPKM is usually not precisely estimated because the number of reads mapped to the gene is quite small. In order to work with reasonably accurate estimates of exon inclusion rates, as pointed out by Jiang and Wong (2009), we calculated inclusion rates only for those exons whose host genes' expressions are above 80^th percentile across at least six samples. Throughout this study, we only considered the genes containing cassette exons whose inclusion rate profiles met the above criterion. This resulted in 12,287 exons covering 8,196 genes. The 38 datasets that met these criteria on January 30, 2011 were used for the analysis described herein.

For each RNA-seq dataset containing a set of samples, a weighted gene co-expression network can be constructed by representing genes as nodes and edges weights as correlations between the expression profiles of two genes. The same procedure is applied to build a weighted exon co-splicing network for the same RNA-seq dataset, in which nodes represent exons and edge weights represent correlations between the inclusion rates of two exons. To determine the weights, we first compute the correlation between two genes (exons) as the leave-one-out Pearson correlation coefficient estimate (Zhou et al., 2005). The resulting correlation estimate is conservative and sensitive to similarities in the patterns, yet robust to single experimental outliers. To make the correlation estimates comparable across datasets, we then applied Fisher's z transform (Anderson, 2003). Given a correlation estimate r, Fisher's transformation score is calculated as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$z = \frac {\sqrt {n - 3}} {2} \ln \left(\frac {1 + r} {1 - r} \right)$$ \end{document} . Note that the sample size n may be different for different datasets and even for different gene (exon) pairs due to missing values. Practically, we observed the distributions of z-scores may still vary from dataset to dataset, we standardized the z-scores to enforce zero mean and unit variance in each dataset by following the normalization procedure introduced in Xu et al. (2008) and Li et al. (2011). Then, the “normalized” correlations r′ are obtained by inverting the z-score with the same virtual sample size n′ = 10: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$r^{\prime} = \frac {\exp (\frac {2} {\sqrt {n^{\prime} - 3}} \,z) - 1} {\exp (\frac {2} {\sqrt {n^{\prime} - 3}} \, z) + 1}$$ \end{document} . Finally, the absolute value of r′ is used as the edge weight of networks.

B. Non-uniform sampling for fast computation

Even though our optimization method is efficient, its computation time can still be long for large sets of networks with many edges. In such cases, edge sampling can provide an efficient approximation to many graph problems (Tsay et al., 1999; Achlioptas and McSherry, 2007). From the perspective of matrix or tensor computation, such sampling methods can also be viewed as matrix/tensor sparsification (Arora et al., 2006). As FCCs predominately contain edges with large weights, we designed a non-uniform sampling method that preferentially selects edges with large weights. Specifically, given a tensor \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal A}$$ \end{document} which could be \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal G}$$ \end{document} (or \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal E}$$ \end{document} ), each edge a_ijk is sampled with probability p_ijk: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} p_{ijk} = \begin{cases} 1, \qquad \qquad \ {\rm if} \ a_{ijk} \tilde {a} \\ p \left(\frac {a_{ijk}} {\tilde {a}} \right)^b, \qquad \; {\rm if} \ a_{ijk} < \tilde {a} \end{cases} \tag {2} \end{align*} \end{document}

where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{a} \in (0, 1), b \in [1, \infty)$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$p \in (0, \tilde{a}^b]$$ \end{document} are constants that control the number of sampled edges. Note that Eq. (2) always samples edges with weights ≥ ã. It selects an edge of weight a_ijk < ã with probability p_ijk proportional to the b^th power of the weight. We choose \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde{a} = 0.4$$ \end{document} , b = 3, and p = 0.1 as a reasonable tradeoff between computational efficiency and the quality of the sampled tensor, meanwhile satisfying the conditions of Theorem 4, i.e., \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$p \, \tilde{a}^{b - 1}$$ \end{document} and b ≥ 1.

To correct the bias caused by this sampling method, the weight of each edge is corrected by its relative probability: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{a}_{ijk} = a_{ijk} / p_{ijk}$$ \end{document} . The expected weight of the sampled network, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$E (\widehat{a}_{ijk})$$ \end{document} , is therefore equal to the weight of the original network. However, in practice, when the adjusted edge weight \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{a}_{ijk} > \tilde{a}$$ \end{document} (but the original edge weight \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$a_{ijk} < \tilde{a}$$ \end{document} ), we enforced it to be \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{a}_{ijk} = \tilde{a}$$ \end{document} to avoid too large edge weights. The overall edge sampling procedure adopts the simple random-sampling based single-pass sparsification procedure introduced in Arora et al. (2006). Details of the sampling procedure is given below. This single-pass sampling procedure's time complexity is O(N²L). It is obvious that the sparsification procedure in Arora et al. (2006) is a special case of our sampling procedure when all entries of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal A}$$ \end{document} are non-negative and p = 1, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\tilde {a} = \frac {\epsilon} {N + N + L}$$ \end{document} , b = 1. After edge sampling, the procedure described above will use the corrected tensor \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{\cal A} = (\widehat{a}_{ijk})_{N \times N \times L}$$ \end{document} instead of the original tensor \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal A}$$ \end{document} .

Based on the well-known Chernoff-Hoeffding bounds (Motwani and Raghavan, 1995), we gave the bound of the non-zero entries in the corrected tensor \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{\cal A}$$ \end{document} after sampling in Lemma 4. Therefore, the computational complexity of the proposed algorithm on the tensors \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{\cal G}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{\cal E}$$ \end{document} after sampling is linear to the number of the non-zero entries of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{\cal G}$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\widehat{\cal E}$$ \end{document} .

C. Proof of Lemma 1

Proof

We can reduce the well-known NP-complete K-clique problem ⁶ (i.e., is there a clique with K nodes in a graph?) to this problem and therefore prove its NP-hardness.

Let G(V, E) be an undirected unweighted graph. We can copy this graph L times to generate a graph set consisting of the L graphs \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal G} = \{G_1 (V, E), \ldots, G_m (V, E) \}$$ \end{document} in which all graphs have the same vertex set V and edge set E. Then we consider our problem with λ = μ = 0 in the set of L graphs \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal G}$$ \end{document} . In other words, our problem with λ = μ = 0 can be rewritten as follows, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & \max \quad H_{\cal G} ({\bf x}, {\bf w}) = \sum \nolimits_{i = 1}^N \sum \nolimits_{j = 1}^N \sum \nolimits_{k = 1}^L{g_{ijk}x_ix_jw_k} \\ & \hbox{subject to}\; \begin{cases} \sum_{i = 1}^N{x_i} = K_1 \\ \sum_{i = 1}^L{w_j} = K_3 \\ x_i \in \{0, 1 \} \quad \qquad \hbox{for any}1 i \le N \\ w_j \in \{0, 1 \} \qquad \quad \hbox{for any}1 j \le L\end{cases} & (8) \end{align*} \end{document}

The solution vectors x* and w* of this problem defines a module having K ₁ nodes whose corresponding x_i = 1 and K₃ networks whose corresponding w_j = 1. Therefore, this module can be called as “recurrent dense subgraph” because its K₁ nodes are most densely interconnected in the K₃ networks. Then the question of “is there a clique with K vertices in G?” can be answered by solving above defined problem, because we can easily claim:

• If the “recurrent dense subgraph” with K nodes and K₃ networks found in the graph set \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal G}$$ \end{document} is a clique with K vertices recurring in the K₃ graphs, then there exist a clique with K vertices in G. This claim is obvious and straightforward.

Since the K-clique problem is NP-complete, the above defined problem is NP-hard and therefore our problem is NP-hard.   ▪

Proof

The Lagrangian function incorporating constraints is formed as follows, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} {\cal L} = H_{\cal G} ({\bf x}, {\bf w}) + \lambda H_{\cal E} ({\bf y}, {\bf w}) + \mu O_R ({\bf x}, {\bf y}, Z) - \tau_1 (\ \parallel {\bf x} \ \parallel_f^f - 1) - \tau_2 (\parallel {\bf y} \ \parallel_g^g - 1) - \tau_3 (\parallel {\bf w} \parallel_h^h - 1) - \sum_{i = 1}^N{\psi_i \left(\sum_{j = 1}^Mz_{ij} - 1 \right)} \tag{9} \end{align*} \end{document}

where τ₁, τ₂, τ₃ and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\psi_i (i = 1, \ldots, N)$$ \end{document} are Lagrangian multipliers for enforcing constraints. According to the Karush-Kuhn-Tucker (KKT) conditions, the partial derivatives of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $${\cal L}$$ \end{document} with respect to x_i, y_j, w_k, and z_ij should be equal to zero: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\frac {\partial {\cal L}} {\partial x_i} = 0, \frac {\partial {\cal L}} {\partial y_j} = 0, \frac {\partial {\cal L}} {\partial w_k} = 0$$ \end{document} and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\frac {\partial {\cal L}} {\partial z_{ij}} = 0$$ \end{document} . The non-negativeness of x, y, w, and Z is enforced by the complementarity condition, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\frac {\partial {\cal L}} {\partial x_i} x_i = 0, \frac {\partial {\cal L}} {\partial y_j} y_j = 0, \frac {\partial {\cal L}} {\partial w_k} w_k = 0$$ \end{document} , and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\frac {\partial {\cal L}} {\partial z_{ij}} z_{ij} = 0$$ \end{document} . This leads to, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & \left(\sum_{j, k}{g_{ijk}x_jw_k} + \mu [(Z \cdot R) {\bf y}]_i - \tau_1fx_i^{f - 1} \right) x_i \qquad \qquad = 0 & (10) \\ & \left(\lambda \sum_{i, k}{e_{ijk}y_iw_k} + \mu [(Z \cdot R)^T{\bf x}]_j - \tau_2gy_j^{g - 1} \right) y_j \qquad \quad = 0 & (11) \\ & \left(\sum_{i, j}{g_{ijk}x_ix_j} + \lambda \sum_{i, j}{e_{ijk}y_iy_j} - \tau_3hw_k^{h - 1} \right) w_k \qquad \ \quad = 0 & (12) \\ & \left(\mu r_{ij}x_iy_j - \psi_i \right) z_{ij} \qquad \qquad \qquad \qquad \qquad \qquad \qquad \ \ \quad = 0 & (13) \end{align*} \end{document}

Summing over index i in Eq. (10), over index j in Eq. (11), over index k in Eq. (12), over index j in Eq. (13) for each index i, we obtain the values of Lagrangian multipliers: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \tau_1 & = \frac {\sum_{i, j, k} {g_{ijk} x_ix_jw_k} + \mu {\bf x}^T (Z \cdot R) {\bf y}} {f} & \qquad\qquad(14) \\ \tau_2 & = \frac {\lambda \sum_{i, j, k} {e_{ijk} y_iy_jw_k} + \mu {\bf x}^T (Z \cdot R) {\bf y}} {g} & \qquad\qquad (15) \\ \tau_3 & = \frac {\sum_{i, j, k} {g_{ijk} x_ix_jw_k} + \lambda \sum_{i, j, k} {e_{ijk} y_iy_jw_k}} {h} & \qquad\qquad(16) \\ \psi_i & = \mu [R \cdot ({\bf x} {\bf y}^T)]^T_{i \ast} Z_{i \ast} \qquad \qquad \qquad \qquad \qquad \qquad i = 1, \ldots, N & \qquad\qquad (17) \end{align*} \end{document}

Note that Eq. (17) is obtained by enforcing the constraint \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\sum\nolimits_{j = 1}^M z_{ij} = 1$$ \end{document} for all i. This doesn't impact the correctness of the algorithm because it can be proved that the constraint \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\sum\nolimits_{j = 1}^M z_{ij} = 1$$ \end{document} does not hurt the optimization when using the constraint \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes {10} {9} {7} {6} \begin{document} $$\sum\nolimits_{j = 1}^M z_{ij} \le 1$$ \end{document} when x and y are enforced by sparsity constraints.

Therefore, it is obvious to verify that the fixed point of the updating steps in Algorithm 1 satisfies Eq. (10,11,12,13). So it is proved.

Proof

We use the auxiliary function approach to prove the theorem. G(x, x′) is an auxiliary function of the function O(x) if \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} G ({\bf x}, {\bf x}^{\prime}) O ({\bf x}) \quad {\rm and} \quad G ({\bf x}, {\bf x}) = O ({\bf x}) \end{align*} \end{document}

We have the following lemma for the auxiliary function.

Lemma 5

If G(x) is an auxiliary function of O(x), then O(x) is nondecreasing under the updates x^(t+1) = arg max _x G(x, x^(t)), because of the fact O(x^(t)) = G(x^(t),x^(t)) ≤ G(x^(t+1),x^(t)) ≤ O(x^(t+1)).

To design a suitable auxiliary function of our Lagrangian function, we used the inequality of x ≥ 1 + log(x) for any positive x which was introduced by Ding et al. (2006) for proving convergence. So the auxiliary function of our Lagrangian function in Eq. (9) is introduced as follows, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} G ({\bf x}, {\bf y}, {\bf w}, Z, {\bf x}^{\prime}, {\bf y}^{\prime}, {\bf w}^{\prime}, Z^{\prime}) & = \sum_{i, j, k} {g_{ijk} x^{\prime}_ix^{\prime}_jw^{\prime}_k \left(1 + \log \frac {x_ix_jw_k} {x^{\prime}_ix^{\prime}_jw^{\prime}_k} \right)} + \lambda \sum_{i, j, k} {e_{ijk} y^{\prime}_iy^{\prime}_jw^{\prime}_k \left(1 + \log \frac {y_iy_jw_k} {y^{\prime}_iy^{\prime}_jw^{\prime}_k} \right)} \\ & \quad + \mu \sum_{r_{ij} \neq 0} {z^{\prime}_{ij} r_{ij} x^{\prime}_iy^{\prime}_j \log \frac {z_{ij} x_iy_j} {z^{\prime}_{ij} x^{\prime}_iy^{\prime}_j}} - \tau_1 \left(\sum_i {x_i^f} - 1 \right) - \tau_2 \left(\sum_j {y_j^g} - 1 \right) \\ & \quad - \tau_3 \left(\sum_k {w_k^h} - 1 \right) - \sum_{i = 1}^N {\psi_i \left(\sum_{j = 1}^Mz_{ij} - 1 \right)} & (18) \end{align*} \end{document}

where x′ denotes x^(t). The maximum for Eq. (18) is obtained by setting its gradient to zero, \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & \frac {\partial G} {\partial x_i} = \frac {x^{\prime}_i \sum_{j, k} {g_{ijk} y^{\prime}_jw^{\prime}_k}} {x_i} + \mu \frac {x^{\prime}_i \sum_{j} r_{ij} z^{\prime}_{ij} y^{\prime}_j} {x_i} - \tau_1 fx_i^{f - 1} = 0 & (19) \\ & \frac {\partial G} {\partial y_j} = \lambda \frac {y^{\prime}_j \sum_{i, k} {e_{ijk} x^{\prime}_iw^{\prime}_k}} {y_j} + \mu \frac {y^{\prime}_j \sum_{i} r_{ij} z^{\prime}_{ij} x^{\prime}_i} {y_j} - \tau_2 gy_j^{g - 1} = 0 & (20) \\ & \frac {\partial G} {\partial w_k} = \frac {w^{\prime}_k \sum_{i, j} {g_{ijk} x^{\prime}_iy^{\prime}_j}} {w_k} + \lambda \frac {w^{\prime}_k \sum_{i, j} {e_{ijk} x^{\prime}_iy^{\prime}_j}} {w_k} - \tau_3 hw_k^{h - 1} = 0 & (21) \\ & \frac {\partial G} {\partial z_{ij}} = \mu \frac {r_{ij} z^{\prime}_{ij} x^{\prime}_iy^{\prime}_j} {z_{ij}} - \psi_i = 0 & (22) \end{align*} \end{document}

The above equations can be re-written as \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} & x_i = \left(x^{\prime}_i \frac {\sum_{j, k} {g_{ijk} x^{\prime}_jw^{\prime}_k} + \mu [(Z^{\prime} \cdot R) {\bf y}^{\prime}]_i} {\tau_1f} \right)^{1 / f} & (23) \\ & y_j = \left(y^{\prime}_j \frac {\lambda \sum_{i, k} {e_{ijk} y^{\prime}_iw^{\prime}_k} + \mu [(Z^{\prime} \cdot R)^T {\bf x}^{\prime}]_j} {\tau_2g} \right)^{1 / g} & (24) \\ & w_k = \left(w^{\prime}_k \frac {\sum_{i, j} {g_{ijk} x^{\prime}_ix^{\prime}_j} + \lambda \sum_{i, j} {e_{ijk} y^{\prime}_iy^{\prime}_j}} {\tau_3h} \right)^{1 / h} & (25) \\ & z_{ij} = z^{\prime}_{ij} \frac {[R \cdot ({\bf x}^{\prime} {\bf y}^{{\prime}^T})]_{ij}} {\psi_i \frac {1} {\mu}} & (26) \end{align*} \end{document}

Since above variables must satisfy the constraints in Eq. (1), we can tune τ₁, τ₂, τ₃, ψ_i like Eqs. (14–17) so that they meet the constraints. Therefore, we can recover the updating rules in Algorithm 1, by letting x^(t+1) = x, y^(t+1) = y, w^(t+1) = w, Z^(t+1) = Z, and x^(t) = x′, y^(t) = y′, w^(t) = w′, Z^(t) = Z′.

To ensure the solution gives a maxima, we show that the second order derivatives (Hessian matrix) of the auxiliary function are negative definite. We have \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland,xspace}\usepackage{amsmath,amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6} \begin{document} \begin{align*} \frac {\partial^2 G} {\partial x_i \partial x_i} & = - \left[\frac {x^{\prime}_i \sum_{j, k} {g_{ijk} y^{\prime}_jw^{\prime}_k}} {x_i^2} + \mu \frac {x^{\prime}_i \sum_{j} r_{ij} z^{\prime}_{ij} y^{\prime}_j} {x_i^2} + \tau_1 f (f - 1) x_i^{f - 2} \right] \\ \frac {\partial^2 G} {\partial y_j \partial y_j} & = - \left[\lambda \frac {y^{\prime}_j \sum_{i, k} {e_{ijk} x^{\prime}_iw^{\prime}_k}} {y_j^2} + \mu \frac {y^{\prime}_j \sum_{i} r_{ij} z^{\prime}_{ij} x^{\prime}_i} {y_j^2} + \tau_2 g (g - 1) y_j^{g - 2} \right] \\ \frac {\partial^2 G} {\partial w_k \partial w_k} & = - \left[\frac {w^{\prime}_k \sum_{i, j} {g_{ijk} x^{\prime}_iy^{\prime}_j}} {w_k^2} + \lambda \frac {w^{\prime}_k \sum_{i, j} {e_{ijk} x^{\prime}_iy^{\prime}_j}} {w_k^2} + \tau_3 h (h - 1) w_k^{h - 2} \right] \\ \frac {\partial^2 G} {\partial z_{ij} \partial z_{ij}} & = - \left[\mu \frac {r_{ij} z^{\prime}_{ij} x^{\prime}_iy^{\prime}_j} {z_{ij}^2} \right] \\ \frac {\partial^2 G} {\partial x_i \partial y_j} & = \frac {\partial^2 G} {\partial x_i \partial w_k} = \frac {\partial^2 G} {\partial x_i \partial z_{ij}} = \frac {\partial^2 G} {\partial y_j \partial w_k} = \frac {\partial^2 G} {\partial y_j \partial z_{ij}} = \frac {\partial^2 G} {\partial w_k \partial z_{ij}} = 0 \end{align*} \end{document}

Because f, g, h ≥ 1, this Hessian matrix is a diagonal matrix with negative quantities on the diagonals. Therefore G(x, y, w, Z, x′, y′, w′, Z′) is a concave function in x, y, w, and Z and has a unique global maximum. This completes the proof of the theorem.   ▪