We constructed about 280,000 scatter diagrams of expression levels of TFs and their regulated genes using DNA microarray data from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) at the NCBI [19]. We selected 135 series of GEO DataSets (GDS) for Homo sapiens that were composed of more than 50 GEO samples. The list of the selected GDS records is in S1 Table. For details of the data preparation, see Methods. To match TFs to regulated genes, we used the TRANSFAC database provided by BioBase and selected 2,073 TF–regulated gene pairs from among 355 TFs and 647 regulated genes. Note that we analyzed TF–regulated gene pairs only with regulations are already confirmed.

We identified four typical types of scatter diagrams depending on the regulatory mechanisms between the TF–regulated gene pairs (Fig 1): constant expression levels for both a TF and the regulated gene, albeit with small fluctuations (no-change type); correlation between expression levels of a TF and the regulated gene (correlation type); no correlation between expression levels of a TF and the regulated gene because the gene has a constant expression level (horizontal-distribution type); and no correlation between expression levels of a TF and the regulated gene because the TF has a constant expression level (vertical-distribution type). We analyzed 2,073 (pairs) ×135 (GDS) to give approximately 280,000 scatter diagrams, and all the diagrams could be classified into these four types.

We studied how the four classification types are implemented. We also characterized each scatter diagram with six indicators to define the classification criteria. Four of them are standard variables: absolute value of slope (|s|) and coefficient of determination (R²) from a least squares approximation, and variance in TF distribution (V_TF) or regulated gene distribution (V_RG). The other two parameters, the uniformity count for a TF (U_TF) or its regulated gene (U_RG), were introduced to distinguish a uniform distribution from the no-change type with a few outliers. The uniformity count is defined as the number of filled units when the area between a maximum and a minimum value is divided into 10 units. These indicators are shown schematically in Fig 2.

After classifying the nearly 280,000 diagrams into the four types (Fig 1), we assigned one correlation type to each regulated gene. It should be noted that one TF can regulate multiple genes and one gene can be regulated by multiple TFs. In addition, one regulated gene can be classified into different types depending on the experimental conditions or cellular states. To avoid ambiguous classifications, we defined a logical rule (see Methods) and assigned one type for each regulated gene depending on the GDS.

Each regulated gene was classified into different relation types depending on the GDS, as shown in Fig 4. However, some genes were classified into one definite type in most GDS. To fix the type for each gene, we integrated the results from the 135 GDS by selecting a majority type from among the correlation, horizontal-distribution, and vertical-distribution types. We ignored the no-change type because our aim was to study how gene expression levels are controlled through the TF–regulated gene correlation. The no-change type occurs when there is no need to change expression levels under the experimental condition of a GDS. For example, PSMB9 was assigned into the correlation type as 53 GDS showed the correlation and only 3 GDS showed the vertical-distribution type (Table 2). In a similar way, CTNNB1 was assigned into the horizontal-distribution type and CCK was assigned into the vertical-distribution type according to the selecting a majority type rule (Table 2). Using the classification criteria, we successfully classified most of the 647 regulated genes: 111 into the correlation type, 178 into the horizontal-distribution type, and 318 into the vertical-distribution type; 40 genes could not be classified because they fell into two majority types (S2 Table).

We performed pathway analysis of the gene functions in the correlation, horizontal-distribution and vertical-distribution types using the curated gene sets in the Canonical pathways (C2:CP) from the Molecular Signatures Database (MSigDB; http://www.broadinstitute.org/gsea/msigdb) [20] with the hypergeometric test at the 1% level of significance (see Methods). We obtained 25 significant pathways for the correlation type, 19 significant pathways for the horizontal-distribution type, and 14 significant pathways for the vertical-distribution type (S3 Table). To compare the different relation types, we categorized the pathways according to the hierarchical framework by denoting a pathway by the top-class entity of its hierarchical framework (Fig 5).

To our surprise, we found that some of the genes in each relation type were associated with type-specific functions: cellular regulation (e.g., Cell Cycle and DNA Replication) for the correlation type, Human Diseases for the horizontal-distribution type, and Metabolism or Signal Transduction for the vertical-distribution type. It is interesting that serious diseases, such as cancers, Parkinson’s disease, Alzheimer’s disease, and Huntington’s disease, were observed in the horizontal-distribution type, i.e., the degradation dominant type. The implications of this observation are considered in the Discussion. To summarize, the scatter diagrams for the TF–regulated gene pairs were characterized systematically according to the functions of the regulated genes. The results indicate that the mechanisms that regulate gene expression levels correspond to gene functions at the whole-cell level.

We used the DNA microarray data from GDS as the expression data in this study. First, we normalized the expression data and removed measurement specificity, generally involving different DNA microarray instruments, to compare the different GDS. Several normalization procedures are available and each has its own advantages [32–37]. In this study, we needed a general-purpose method applicable to various measurement platforms and used a popular method, Z scores [33, 37, 38], as follows.

For each sample in a GDS, we first transformed the original expression data given as {x₁,x₂,…,x_s} to the log-scale, {ln(x1),ln(x2),⋯,ln(xs)}.(2)

Then, we normalized the values using the Z-score method by defining E=1s∑i=1sln(xi),V=ss−1{1s∑i=1s(ln(xi))2−E2}.(3)

The normalized value is given as ln(xi)-EV(4) for every x_i.

We constructed scatter diagrams for expression levels of each TF (TF_i) and its regulated gene (RG_i) from a GDS. Suppose a GDS contains N_i(⩾50) samples, then the diagram has N_i points as explained below. When a sample had only one data point for TF_i (RG_i), we used this value for plotting, and when a sample contained more than one data point for TF_i (RG_i), we used the average value. Thus, one sample yielded one point, and the diagram had N_i points in total.

One GDS normally contains subclasses such as an experimental group and a control group. It could be that each subclass produces a different domain structure in the diagram and falsely represents an imaginary correlation. Namely, if the samples in one subclass show smaller TF_i and RG_i and the samples in another subclass show higher TF_i and RG_i because of the experimental conditions, a correlation may be observed between TF_i and RG_i even if there is no real correlation. We confirmed that such imaginary correlations appeared rarely and did not influence the results.

We considered a general and simple mathematical model of a transcription process. We analyzed two situations: TF promoting gene expression (up-regulation), and TF suppressing gene expression (down-regulation). First, we explain the up-regulation case in detail and next the down-regulation case in brief.

For the up-regulation case, suppose a TF molecule stochastically binds to or dissociates from a promoter sequence, and transcription takes place only when the TF binds to the sequence. Then, P_b is the probability of the TF’s binding to the promoter sequence, and is defined as a fraction of bound TF molecules among all TF molecules. By assuming that the binding process and dissociation process are in equilibrium, we get the following equation: kbTF(1-Pb)=kuPb.(5)

Here, [TF] represents the TF concentration, and k_b and k_u are the reaction coefficients for the binding and dissociation processes, respectively. The left side of Eq (5) represents the reaction rate of the TF binding process proportional to the product of the TF concentration ([TF]) and the unbound promoter sequence (1−P_b). On the other hand, the right side represents the TF dissociation reaction rate regulated by the bound promoter sequence(P_b). From Eq (5), we obtain Pb=TF/K1+TF/K(6) where the dissociation constant K = k_u/k_b. Because transcription occurs only when the TF binds to the promoter region, the mRNA synthesis rate of the regulated gene should be proportional to P_b. Therefore, we can write the time dependence of the expression of the regulated gene mRNA ([RG]) as dRGdt=aTF/K1+TF/K-bRG.(7)

Here, a and b are reaction coefficients for the synthesis and degradation. By considering a steady state of Eq (7), d[RG]*dt=0, we finally obtain the steady-state mRNA level of the regulated gene ([RG]*) as a function of the steady-state TF concentration ([TF]*) as shown in Eq (1), RG*=1γTF*/K1+TF*/K.(1)

For the down-regulation case, we obtain the following equation from a similar analysis except that the production rate is proportional to the dissociation probability (1−P_b): dRGdt=a′(1-TF/K′1+TF/K′)-b′RG(8) and therefore RG*=a′b′11+TF*/K′.(9)

Eq (9) describes the same two types of characteristic behaviors as Eq (1) depending on the dissociation constant K′, although [RG]* shows a strong negative correlation with [TF]* when K′ ≪ 1 and [RG]* remains at a constant level regardless of [TF]* when K′ ≫ 1. In the up-regulation or down-regulation cases, the correlation between a TF and the regulated gene indicates fine-tuned rate regulation, whereas the horizontal distribution indicates the absence of regulation.

We defined the criteria for classifying the scatter diagrams into the four types as follows. First, we excluded the data with V_TF = 0 or V_RG = 0 because they probably originate from a measurement flaw. It is virtually impossible for all the samples in a GDS to show exactly the same expression level of a gene. We also assumed that the expression level of a TF (or a regulated gene) changed significantly when V_TF > 0.25 (V_RG > 0.25), whereas it is constant, albeit with small fluctuations, when V_TF ⩽ 0.25 (V_RG ⩽ 0.25). After that, we classified those with 0.1 < |s|<10 and R² > 0.3 into the correlation type when V_TF > 0.25 and V_RG > 0.25. We then classified diagrams with |s|<0.1, V_TF > 0.25 (V_RG > 0), V_TF/V_RG > 3, and U_TF ⩾ 7 into the horizontal-distribution type, and diagrams with |s|>10, V_RG > 0.25 (V_TF > 0), V_TF/V_RG < 1/3, and U_RG ⩾ 7 into the vertical-distribution type. The remaining diagrams were assigned to the no-change type. The baseline values used here were set arbitrarily, but the discussion will not change if the values are changed to some extent.

In many cases, TFs and regulated genes do not have a one-to-one correspondence. When a regulated gene has several TFs, some of the TFs finely regulate the transcription process whereas others simply switch the process on or off. The former TF–regulated gene pairs may match the correlation type, whereas the latter often correspond to the vertical-distribution type. In such a mixed case, the regulated gene should be classified into the correlation type not into the vertical-distribution type. Similarly, when TF–regulated gene pairs match the horizontal-distribution type and others correspond to the no-change type, the regulated gene should be classified into the horizontal-distribution type. Using these rules, we classified every regulated gene as follows.

Suppose a gene has M TFs (M ⩾ 1) and among them, M_n TFs are of the no-change type, M_c of the correlation type, M_h of the horizontal-distribution type, and M_v TFs are of the vertical-distribution type. When M_h > M_c and M_v, the regulated gene is classified into the horizontal-distribution type; when M_v ⩾ M_h ⩾ M_c = 0, the regulated gene is classified into the vertical-distribution type; and when M_c > M_h ⩾ 0, the regulated gene is classified into the correlation type, regardless of M_v. This is because a TF, even if it serves as a single TF, can determine the correlation type as explained above. We classified a regulated gene into the no-change type only when M_n = M. For the remaining cases, we aborted the classification because there was not sufficient evidence. We thus assigned one classification type to each regulated gene depending on the GDS.

We determined whether a list of genes (genes of each relation type) over-represents a biological process (gene sets for a pathway from MsigDB) using the hypergeometric test. Suppose we listed n genes from a total of N genes; i.e., we selected n genes without replacement from the N genes. M genes among the total of N genes are involved in the biological process, and m genes among the listed n genes are involved in the same process. Then, the probability distribution of m (p(m)) is described by the hypergeometric distribution as p(m)=(Mm)(N-Mn-m)(Nn).(10)

Our goal is to determine whether the case of m genes being involved in the biological process (out of the n listed genes) is statistically significant or happened by chance. Because we are testing whether our gene set corresponds to over-representation, the hypergeometric p value (p) is calculated as the probability of random involvement of m or more genes in the biological process (out of n genes) and is expressed as p=∑x=mnp(x)=∑x=mn(Mx)(N-Mn-x)(Nn).(11)

When the p value is less than the value we set as the level of significance (1%), we conclude that our set of genes is over-represented, i.e., the m genes occurred non-randomly. However, when the p value is greater than the threshold value, we conclude that the m genes are selected by chance.