Using graph models to find transcription factor modules: the hitting set problem and an exact algorithm

Data sets

Finding cooperative TFs for a set of co-regulated genes

Given a set of co-regulated genes, finding a set of cooperative TFs regulating them is a challenging problem in studying transcriptional regulation [2]. In this study, we cast the task as a graph problem, referred to as the weighted t-cover hitting set problem, and designed an efficient algorithm to solve it.

For a set of co-regulated genes, we induce a subgraph from the bipartite graph representing protein-DNA interaction. The new graph remains a bipartite graph, with one part being all TFs and the other being the given co-regulated genes. If a TF is connected to a gene, we say that the TF covers the gene. The task is to find a subset of TFs with minimum weight, called the t-TF cover, such that each gene is covered by at least t TFs in the cover. The weight of the t-TF cover is the sum of the weights of TFs in the cover. We defined the weight of a TF to be the reciprocal of the sum of edge weights, which is defined in the input bipartite graph by Huang et al.,[9, 10]. More specifically, if TF t₁ is connected to genes g₁, g₂, …, g _k , then $\mathit{\text{weight}}\left(t_1\right)=\frac{1}{\sum _{i=1}^k\mathit{\text{Edge}}\text{\_}\mathit{\text{Weight}}(t_1,g_i)}$. Recall that, in our case, an edge weight reflects the binding reliability between the TF and its target gene; as such, a TF which interacts with many genes with high reliability will have a small TF weight, and thus is more likely to be included in the solution of the t-TF cover problem.

When t = 1, the solution finds a set of TFs such that each gene must be covered by at least one TF, which will return a minimum set of TFs that covers all those co-regulated genes. Such a result is not interesting because we aim to find cooperative TFs. By constraining t to be greater than 1, we are able to try to find a set of TFs such that each gene must be covered by more than one TF in the set. By definition, a solution will preferentially search for and include TFs that co-cover as many genes in a target set as possible and with as high reliability as possible. From a biological viewpoint, the more genes a set of TFs co-covers, the more likely the TFs act cooperatively to regulate the genes. Thus, finding an optimal solution for the t-TF cover problem in this setting is a biologically sensible approach for identifying TF modules.

The problem of finding the t-TF cover is NP-hard, where the problem is equivalent to two well-known NP-hard problems, the set multicover problem and the t-cover hitting set problem [11, 12]. The t-TF cover problem can be reduced to the t-cover hitting set problem easily. The formal definition of the weighted t-cover hitting set is as follows:

weighted t-cover hitting set: Given a universal set X of m elements, a weight function w(x) : X → R⁺, a family $\mathcal{T}=\{S_1,\dots ,S_n\}$ of subsets of X, where the size of any subset in $\mathcal{T}$ is at least t, and an integer t, find a subset H ⊆ X of minimum weight such that every subset in $\mathcal{T}$ has at least t elements in H, where the weight of H is defined as $\sum _{x\in H}w\left(x\right)$. We denote an instance of the problem as $(X,\mathcal{T},w,t)$ and call H the minimum t-cover hitting set.^a

As mentioned above, the weighted t-cover hitting set problem is NP-hard and the previous best algorithm for the unweighted case of the problem has a time complexity of O((t + 1) ⁿ n m) [12]. In terms of our application, n is the number of co-regulated genes and m is the number of TFs; these two numbers are usually large enough to render the existing algorithms impractical for our problem. Thus, designing an efficient exact algorithm to solve the weighted t-cover hitting set problem becomes the main goal and therefore the main contribution of this study.

Identify repeatedly used TF modules

For each perturbation instance, a set of differentially expressed genes (considered to be co-regulated genes) can be decided. Using the method in the previous subsection, we can obtain a set of cooperative TFs for each perturbation instance. At this stage, a biologist may want to find modules of TFs that are repeatedly used under various experimental conditions, as they are more likely to be part of signal transduction pathways that are repeatedly involved in responses to different perturbations.

We refer to a set of TFs sharing a common set of target genes as a “hard” TF clique. The formal definition of a “hard” TF clique is as follows: a subset of TFs T forms a “hard” TF clique if there exists a subset of genes G such that T and G make a complete induced subgraph in the protein-DNA interaction graph. From a biological viewpoint, the more common target genes are shared by the TFs, the more likely the TFs truly cooperate to regulate the genes. Thus, we weigh a hard clique using the number of common target genes covered by the TFs to reflect the “goodness” of the clique. In such a setting, one way to find the TF modules that are repeatedly used in many perturbation instances is to find hard TF cliques that occur in multiple instances, where the overall weight of such hard cliques is the sum of weights of the clique from all instances. The higher the weight of a clique, the more likely the TFs in the clique will function as a module. However, as there are usually noise and errors in microarray data, it is difficult to consistently find hard cliques from multiple instances. To address this issue, we introduce a new formulation: search for a TF module whose member TFs tend to cooperatively regulate gene expression in multiple instances but not necessarily in all perturbation instances.

First, we use the cooperative TFs found in each instance to make one TF-TF-relation graph for each case of t = 1, 2, 3, 4. In such a graph, each TF is a node; a weighted edge between a pair of TFs is added if the TFs are found to cooperate in one or more instances, of which the weight is the sum of the number of common target genes from all instances. Hence, the higher the weight of an edge, the more cooperation instances the two TFs have. Then we search for all 3-cliques and 4-cliques from such a TF-TF graph and sort them by weights, where the weight of a clique is the sum of all edge weights for the clique. Since we have relaxed the requirement so that the TFs in these cliques do not need to form a hard clique in all individual instances, we refer to this approach as finding “soft” TF cliques. At this stage, we limit the search for cliques to 3-cliques and 4-cliques because our data [9, 10] indicated that only 15% of genes in yeast are connected to more than 5 TFs.

The clique finding problem is a well-known NP-hard problem. However, since we only need to find cliques of size 3 and size 4, we can search for cliques of these sizes in a graph with at most 214 (the total number of TFs in our case) vertices in a reasonable time. In the Results section, we will compare the results for “soft” 3-cliques for the cases of t = 1, 2, 3, 4 and hard 3-cliques for the case of t = 3, and “soft” 4-cliques for the cases of t = 1, 2, 3, 4 and hard 4-cliques for the case of t = 4.

Lemma 0.1

Let H₁, H₂, H^′ be three subsets of X such that H₁ ∩ H^′ = ∅ and H₂ ∩ H^′ = ∅. If h i t(H₁) = h i t(H₂), then h i t(H₁ ∪ H^′) = h i t(H₂ ∪ H^′).

Proof

As h i t(H₁) = h i t(H₂), for any $S_i\in \mathcal{T}$, h i t(H₁)[S _i ] = h i t(H₂)[S _i ], i.e., m i n(t, |S _i ∩H₁|) = m i n(t, |S _i  ∩ H₂|). Furthermore, because H₁ ∩ H^′ = ∅ and H₂ ∩ H^′ = ∅, we will have that, for any $S_i\in \mathcal{T}$, m i n(t, |S _i  ∩ (H₁ ∪ H^′)|) = m i n(t, |S _i ∩H₁| + |S _i  ∩ H^′|) = m i n(t, |S _i  ∩ H₂| + |S _i ∩H^′|) = m i n(t, |S _i ∩(H₂∪H^′)|). Therefore, h i t(H₁∪H^′) = h i t(H₂∪H^′) and the lemma is proved. □

The Lemma 0.1 guarantees that if any two sub-solutions cover in the same way, then keeping the sub-solution with the smaller weight is enough.

Lemma 0.2

Let $H^{\ell }=\{u_{i_1},u_{i_2},\dots ,u_{i_{\ell }}\}$, whose elements are in the same order as in the sorted X in Algorithm-1 (i.e., if j₁ < j₂ with respect to the index of H ^ℓ , then $i_{j_1}<i_{j_2}$ with respect to the index of X), be the minimum t-cover hitting set, and $H_j^{\ell }=\{u_{i_1},u_{i_2},\dots ,u_{i_j}\}$, 1 ≤ j ≤ ℓ. For any 1 ≤ j ≤ ℓ, if there is a $H\subset \{u_1,u_2,\dots ,u_{i_j}\}$ such that $\mathit{\text{hit}}\left(H\right)=\mathit{\text{hit}}\left(H_j^{\ell }\right)$, then $\mathit{\text{weight}}\left(H\right)\ge \mathit{\text{weight}}\left(H_j^{\ell }\right)$.

Proof

Let $H^{\prime }=\{u_{i_{j+1}},u_{i_{j+2}},\dots ,u_{i_{\ell }}\}$. Then H ∩ H^′ = ∅. If $\mathit{\text{weight}}\left(H\right)<\mathit{\text{weight}}\left(H_j^{\ell }\right)$, then by Lemma 0.1, H ∪ H^′ will be a t-cover hitting set with a smaller weight than the weight of H ^ℓ , which causes contradiction. Hence, the lemma is correct. □

The Lemma 0.2 shows that Algorithm-1 always keeps a sub-solution that will lead to the full-solution with minimum weight. Now, let us present and prove the main theorem.

Theorem 0.3

The weighted t-cover hitting set problem can be solved in O((t + 1) ⁿ m n t) time and in O((t + 1) ⁿ n t) space, where m is the size of the ground set and n is the number of subsets for the given instance. If, furthermore, the problem has at least $\frac{n}{1+d/\underset{2}{log}(t+1)}$ subsets whose sizes are upper bounded by d, then the problem can be solved in $O\left({\left({(t+1)}^{d/(d+\underset{2}{log}(t+1\left)\right)}\right)}^n\mathit{\text{mnt}}\right)$ time and in $O\left({\left({(t+1)}^{d/(d+\underset{2}{log}(t+1\left)\right)}\right)}^n\mathit{\text{nt}}\right)$ space.

Proof

We first prove the correctness of the algorithm.

Given an instance $(X,\mathcal{T},w,t)$ of the weighted t-cover hitting set problem, let X = {u₁, u₂, …, u _m }, where X is sorted as shown in Algorithm-1 such that the order of elements in X is as S₁, S₂ - X₁, …, S _n -X_n-1, where $X_j={\cup }_{i=1}^j{S}_i$. Let $H^{\ell }=\{u_{i_1},u_{i_2},\dots ,u_{i_{\ell }}\}$, whose elements are in the same order as in the sorted X (i.e., if j₁ < j₂ with respect to the index of H ^ℓ , then $i_{j_1}<i_{j_2}$ with respect to the index of X), be the minimum t-cover hitting set. Let $H_j^{\ell }=\{u_{i_1},u_{i_2},\dots ,u_{i_j}\}$ for all 0 ≤ j ≤ ℓ, where $H_0^{\ell }=\varnothing$.

To prove correctness, we claim that when the for loop in step 2 of Algorithm-1 is at loop i = i _j for all 1 ≤ j ≤ ℓ (Note: $u_{i_j}$ is the j th element in H ^ℓ and i _j th element in X), there exists a P = (h i t(H), H) in ${\mathcal{Q}}_{\mathit{\text{old}}}$ (loop in step 3) such that $\mathit{\text{hit}}\left(H\right)=\mathit{\text{hit}}\left(H_{j-1}^{\ell }\right)$ and $\mathit{\text{weight}}\left(H\right)=\mathit{\text{weight}}\left(H_{j-1}^{\ell }\right)$. We prove this claim by mathematical induction on j.

Therefore, when j = ℓ, we will save a (h i t(H), H) into ${\mathcal{Q}}_{\mathit{\text{new}}}$ such that $\mathit{\text{hit}}\left(H\right)=\mathit{\text{hit}}\left(H_{\ell }^{\ell }\right)=\mathit{\text{hit}}\left(H^{\ell }\right)$, and w e i g h t(H) = w e i g h t(H ^ℓ ), i.e., we will find the minimum t-cover hitting set. The correctness of Algorithm-1 is proved.

Next, we consider the time complexity and space complexity of the algorithm. Step 2 loops |X| = m times. Step 3 loops $\left|{\mathcal{Q}}_{\mathit{\text{old}}}\right|$ times. As ${\mathcal{Q}}_{\mathit{\text{old}}}$ only remember different combinations of [c₁, c₃, …, c _n ] and each c _i is between 0 and t, it is obvious that $\left|{\mathcal{Q}}_{\mathit{\text{old}}}\right|\le {(t+1)}^n$. Steps 4.1 to 4.4 take O(n t) time. Steps 4.6 to 4.11 can be finished in O(log2(t + 1) ⁿ ) = O(n log2(t + 1)) time if we use AVL tree to implement ${\mathcal{Q}}_{\mathit{\text{new}}}$ and ${\mathcal{Q}}_{\mathit{\text{old}}}$. Hence the total time complexity is O((t + 1) ⁿ m n t).

In the case that $\mathcal{T}$ has at least $\frac{n}{1+d/\underset{2}{log}(t+1)}$ subsets whose sizes are bounded from above by d, then when $i=\frac{\mathit{\text{dn}}}{1+d/\underset{2}{log}(t+1)}$, both ${\mathcal{Q}}_{\mathit{\text{old}}}$ and ${\mathcal{Q}}_{\mathit{\text{new}}}$ have at most $2^{\frac{\mathit{\text{dn}}}{1+d/\underset{2}{log}(t+1)}}={\left({(t+1)}^{\frac{1}{1+\underset{2}{log}(t+1)/d}}\right)}^n$ elements. Furthermore, when $i=\frac{\mathit{\text{dn}}}{1+d/\underset{2}{log}(t+1)}$, for any P = (h i t(H), H) in ${\mathcal{Q}}_{\mathit{\text{old}}}$ or in ${\mathcal{Q}}_{\mathit{\text{new}}}$, if let h i t(H) = [c₁, c₂, …, c _n ], then c _j  = t for $1\le j\le \frac{n}{1+d/\underset{2}{log}(t+1)}$. Hence, when $i>\frac{\mathit{\text{dn}}}{1+d/\underset{2}{log}(t+1)}$, all elements in ${\mathcal{Q}}_{\mathit{\text{old}}}$ or in ${\mathcal{Q}}_{\mathit{\text{new}}}$ have at most ${(t+1)}^{n-\frac{n}{1+d/\underset{2}{log}(t+1)}}$ combinations of h i t(H), i.e., the size of ${\mathcal{Q}}_{\mathit{\text{old}}}$ or ${\mathcal{Q}}_{\mathit{\text{new}}}$ is always bounded from above by ${(t+1)}^{n-\frac{n}{1+d/\underset{2}{log}(t+1)}}={\left({(t+1)}^{\frac{1}{1+\underset{2}{log}(t+1)/d}}\right)}^n$. Therefore, the total time complexity is $O\left({\left({(t+1)}^{\frac{1}{1+\underset{2}{log}(t+1)/d}}\right)}^n\mathit{\text{mnt}}\right)$.

It is obvious that the space complexity is $O\left(\right|{\mathcal{Q}}_{\mathit{\text{old}}}|·max$ (lengthes of elements in ${\mathcal{Q}}_{\mathit{\text{old}}}\left)\right)=O\left(\right|{\mathcal{Q}}_{\mathit{\text{new}}}\left|\right)·max$ (lengthes of elements in ${\mathcal{Q}}_{\mathit{\text{new}}}\left)\right)$. The lengthes of elements in both ${\mathcal{Q}}_{\mathit{\text{old}}}$ and ${\mathcal{Q}}_{\mathit{\text{new}}}$ are bounded from above by O(n t). Therefore, in the general case, the space complexity is O((t + 1) ⁿ n t) and in the case sizes of many subsets in $\mathcal{T}$ are bounded from above by d, the space complexity is $O\left({\left({(t+1)}^{d/(d+\underset{2}{log}(t+1\left)\right)}\right)}^n\mathit{\text{nt}}\right)$. □

The Algorithn-1 only reports one solution with the minimum weight, even problems in application have multiple solutions with the minimum weight. The setting of weights of TFs increases the probability that any solution with the minimum weight includes most correct TFs that regulate differently expressed genes. However, in some cases of the application, we may also want to study other top weight solutions (as the data error, the actual solution may not have the minimum weight). By modifying the algorithm such that for each distinct cover way, save k top weight sub-solutions, then the new algorithm can output k top weight solutions. It is also easy to prove that the time complexity and space complexity of the new algorithm will only increase by a ratio k.

Before we finish this section, we briefly summarize the time complexity of our algorithm and compare it with the previously reported one [12]:

TFs for sets of co-regulated genes

For each yeast genetic/pharmacological perturbation experiment, we identified genes that were differently expressed. Then, we applied our weighted t-cover hitting set algorithm to each bipartite component induced by connecting the differentially expressed genes and TFs to identify a set of cooperative TFs from the component. In addition, we applied our algorithm by setting t to 2, 3, and 4 respectively to investigate the impact of setting this parameter. Empirically, we found that setting t = 2 was sufficient to force our algorithm to find a set of cooperative TFs. We recommend setting t=2 as a beginning parameter and exploring other settings based on the degree of the genes in the organism of interest. After inspecting the resulting TF modules, we found many of them to be already well-known. Due to page limitations, we will discuss just one of the well-known modules identified by our algorithm below.

The pheromone response pathway of S. cerevisiae, which consists of more than 20 proteins [13], is a well-studied signal transduction pathway in yeast. When this pathway is activated by pheromone, a well-studied transcription program is initiated which is known to be cooperatively regulated by TFs: S t e 12p, D i g 1p/D i g 2p[13] or S t e 12p, D i g 1p/D i g 2p, and M c m 1p[9]. In the experiments by Hughes et al., 12 gene-encoding proteins in the pathway were perturbed, and many canonical pheromone response genes were differentially expressed. Hence, one may expect that S t e 12p, D i g 1p or S t e 12p, D i g 1p, and M c m 1p are involved in mediating these responses. Indeed, our results show that, when we set t = 2, S t e 12p, D i g 1p were returned as members of the set of cooperative TFs identified in all those 12 perturbations and S t e 12p, D i g 1p, M c m 1p were found in the set of cooperative TFs in 9 out of those 12 perturbation experiments. In comparison, when we set t = 1, S t e 12p, D i g 1p were found in 5 instances, and S t e 12p, D i g 1p, M c m 1p were found in only 1 instance. These results indicate that our algorithm is capable of identifying cooperative TFs from individual experiments in a context-specific manner.

We further studied the coverage of 7 robust pheromone response genes (i.e., their expression levels change significantly in almost all 12 perturbations of genes on the pheromone response pathway), namely F A R 1, F U S 1, G P A 1, S S T 2, S T E 2, S T E 6, T E C 1, to investigate how they were covered by TFs in the results obtained from the t-cover hitting set algorithm. Figure 3 shows the coverage of these genes by the TFs in the results returned by t-TF cover algorithm. As expected, when t = 1, the algorithm failed to find cooperative TFs but returned S t e 13p as the regulating TF. On the other hand, when we set t = 2, the algorithm returned all three TFs, S t e 12p, D i g 1p, and M c m 1p, as the members of solution TFs covering these genes, which form a dense graph.

Above examples show that if we know a set of genes that are co-regulated, our new program can find correct TFs that regulate them when we set t to 2 or 3. There exist more examples to indicate that the algorithm is working as expected.

Finding TF modules involved in multiple instances

Using the t-cover hitting set results (for t = 1, 2, 3, 4) from 300 microarray data, we constructed a TF-TF relation graph and searched for 3-cliques and 4-cliques in the graph. The cliques are ranked according to their weights; the larger the weight, the higher the rank is. Our assumption is that, because high-ranking cliques consist of the TFs deemed to function cooperatively in multiple instances by our algorithm, they may truly function as partners in real biological settings. We evaluated this assumption looking at two factors: 1) whether members of the cliques are known to function as a module based on literature, and 2) whether the members of the cliques have shown physiological or genetic interactions in previous experiments.

1. High scores for known cooperative TFs: We first inspected the top 20 high-ranking 3 cliques identified in the TF-TF graph derived from the results of the t-cover hitting set with t=2. We found that the majority of them are well known cooperative TF modules. For example, TFs participating in cell cycle checkpoint [9], S w i 4p, S w i 6p, and M b p 1p, ranked 2nd and, interestingly, the 3-clique including D i g 1p, M c m 1p, and S t e 12p (the pheromone responding TF module) ranked 4th. Similarly, many of the high ranking 4-cliques constitute TF modules that are also supported by existing knowledge. However, it is interesting to note that the 3-clique that ranked the highest consists of D i g 1p, S t e 12p and T e c 1p, for which we could not find any published results reporting that they work together; thus it would be an interesting case to investigate if our approach has found previously unknown interactions.

We further investigated if the input data, the TF modules returned using the t-cover hitting sets with different settings of t, has an impact on the results of TF cliques. We found that the ranking and the members of cliques were similar when t was set to 2, 3, and 4, respectively. However, the results for when t = 1 and for hard cliques were totally different. For example, when the TF-TF graph is built using the results with t = 1, the first 3-clique containing S w i 4p, S w i 6p, and M b p 1p ranks 254th. Furthermore, if we rank this commonly used TF module by finding common “hard” cliques, the highest ranking 3-clique containing D i g 1p, M c m 1p, and S t e 12p ranks 273rd. Thus, the results indicate that, by setting t > 1 and finding “soft” cliques, our methods enhance the capability of finding TF modules that are repeatedly utilized in multiple conditions, and are thus likely to play key roles in cellular signal transduction systems. Due to page limitation, we made the results available by listing the top 100 soft 3-cliques and 4-cliques on the supplement.

2. Interactions between proteins inside the cliques: From a biological point of view, if a set of TFs work together to regulate the transcription processes, there should be physical or genetic interactions among those TFs. If the results of the t-cover hitting set algorithm truly capture the cooperations, one would expect that TFs in the repeatedly utilized cliques have a high probability of interacting. The following section evaluates the TF cliques through analyzing their protein-protein interactions, where we define that there exists an interaction between a pair of TFs if they have a physical interaction, a genetic interaction determined by synthetic lethality [14], or both.

We constructed multiple TF-TF relation graphs using the results from the t-cover hitting set algorithm with t = 1, 2, 3, 4 to assess the impact of setting t. From each graph, we identified the top 100 4-cliques and assessed the percent (aka, the probability) that the cliques have at least a given number of interactions. As a control, we also randomly sampled 100 groups of 4 TFs for comparison. The results are plotted as a cumulated distribution in Figure 4. From the figure, we can see that the majority of randomly picked TF groups have no interaction, with less than 20% of groups containing 1 or more interaction(s). Similarly, the curve for the cliques that are derived from the results with t = 1 —a setting that is not optimized for finding cooperations among TFs—is located close to that curve of the random groups. On the other hand, with t = 2, 3, 4, over 50% of 4-cliques contain at least 3 interactions. These results indicate that when one sets t > 1, the t-cover hitting set algorithm indeed strives to find the TFs that cooperatively regulate transcriptions; the results also indicate that setting t = 2 is sufficient to find cooperative TFs.

To illustrate the difference in performance between the “soft” and “hard” cliques, we also identified the top 100 “hard” 4-cliques that were obtained directly from the t-cover hitting set for individual perturbation instances by setting t = 4 and studied interactions among the TFs in these cliques. The results show that the performance of the “hard” cliques is superior to those of the random TF groups and the “soft” cliques derived from t = 1 results, indicating that the t-cover hitting set algorithm is capable of revealing cooperation among TFs. However, “hard” cliques consistently underperform in comparison with the “soft” cliques that are derived from integrating results in multiple instances (i.e., except the case of t = 1), which indicates that the information integration approach further enhances the quality of the TF modules.