TF-finder: A software package for identifying transcription factors involved in biological processes using microarray data and existing knowledge base

Identification of salt stress response and tolerance regulators

We applied ASCCA to 109 microarray data sets collected from seven salt stress microarray experiments. The input files contain the expression profiles of 159 positive target genes (non-TFs, Additional file 1) that are known to be involved in salt response and tolerance, 1638 Arabidopsis TFs present in Affymetrix ATH1 array, and 13 TFs (AT1G01520, AT2G40950-BZIP17, AT5G39610-ATNAC2, AT5G67450-AZF1, AT3G19580, AT1G52890-ANAC019, AT1G35515-HOS10, AT2G47190-MYB2, AT2G27300-NTL8, AT3G55980-SZF1, AT2G30250-WRKY25, AT2G38470-WRKY33, AT4G28110-MYB41) (Additional File 1) known to be involved in salt response and tolerance. The cluster analysis of the 159 target genes resulted in about 800 clusters that were used to hook TFs in a recursive manner. All TFs identified through this procedure were pooled for frequency calculation. The top 70 genes with highest occurrence frequencies are shown in Additional file 2. Among these genes, 17 TFs were clearly supported by existing evidence to be involved in salt response and tolerance (Additional file 2). For example, WRKY33, AZF2, and NATAC6 were among the list of 13 TFs used as guide genes (Additional file 1). Although the other 14 were all novel, indirect evidence suggests that they are likely involved in this stress response. For instance, CZF1, also known as SZF2, is the most homologous gene to SZF1, and it regulates salt stress responses in Arabidopsis[8]. ZAT6 is the most homologous gene to STZ (salt tolerance zinc finger) in Arabidopsis. RHL41 (also called ZAT12) is involved in hyperosmotic salinity response [9]. ANAC055 has been found to bind to the early responsive to dehydration (ERD1) stress gene promoter, and over-expression of this gene, together with ANAC019 and ANAC072, causes the expression of several stress-inducible genes that enhance drought tolerance [10]. Over-expression of SZF1 (Salt-inducible zinc finger 1) in transgenic plants caused reduced induction of salt responsive genes and increased tolerance to salt [8]. STZ (salt tolerance zinc finger) was found to increase salt tolerance of calcineurin mutants of wild-type yeast, which appears to be partially dependent on ENA1/PMR2, a P-type ATPase required for Li⁺ and Na⁺ efflux. ATAF1 is responsive to wounding and ABA. DREB2A and DREB2B (DRE/CRT-binding protein) are induced upon dehydration and high salinity [11]. ATMYC2 is a positive regulator of ABA signalling. MYBR1 is ABA-regulated and participates in mediating ABA effects [12]. CBF1 functions as a transcriptional activator that binds to the C-repeat/DRE DNA regulatory element in response to low temperature and water deficit [13]. Although CBF1 mainly responds to chilling, the expression of CBF1 also confers salt stress tolerance [14]. BZIP28, an ER-resident TF, serves as a sensor/transducer in Arabidopsis to mediate ER stress responses [15].

Identification of adaptive growth regulators under salt condition

We also used TF-finder to identify TFs controlling growth under the same stress condition. We used the expression profiles of 74 positive target genes that are involved in growth, 10 positive TFs (Additional file 1), and 1640 TFs. 10 positive TFs include AT5G02470-DPA, AT4G16110-ARR2, AT3G13960-GRF5, AT5G53660-GRF7, AT2G16720-MYB7, AT3G49690-MYB84, AT5G20730-NPH4, AT1G13260-RAV1, AT2G33880-HB3, and AT1G32640-MYC2. The resulting top 70 candidate TFs are shown in Additional file 2. Among these genes, 26 TFs have regulatory functions in growth. Three NAC domain-containing TFs: AT3G61910-ANAC066, AT1G60280-ANAC023, and AT5G04400-ANAC077, have been shown to be involved in the differentiation and expansion of petals, stamen, and roots [16–18]. Three closely related basic helix-loop-helix (bHLH) proteins, AT5G53210-SPCH, AT3G06120-MUTE and AT_FAMA, have been identified as positive regulators that direct three consecutive cell-fate decisions during stomatal development [19, 20]. AT3G13960-AtGRF5 is one of the nine members of GRF gene family that contain nuclear targeting domain, and is involved in root development [21]. AT2G13570-NF-YB7 encoding LEAFY COTYLEDON1-LIKE is a regulator essential for embryo development [22, 23]. KNAT6 is expressed in roots and is required for proper lateral root formation[24]. AT4G27330-SPL plays a central role in patterning of both the proximal-distal and the adaxial-abaxial axes in the ovule and is generally involved in cell differentiation [25]. AT2G35670-FIS2 and AT1G02580-MEA are involved in seed development [26]. CAL is floral homeotic gene encoding a MADS domain protein homologous to AP1 promoting the flower to shoot transformation in ap1 mutants [27]. AT3G15170-CUC1, together with CUC2 and CUC3, are responsible for shoot organ boundary and meristem formation throughout the different stages of Arabidopsis life cycle [28, 29]. NUB encodes a protein with a single C(2)H(2) zinc-finger domain and is involved in the growing of later organs [30]. DOT5 is involved in vein patterning, but dot5-1 mutants often have shorter roots, suggesting its functions in root development [31]. INO is involved in ovule development [32]. BLH8 encoded a BEL1 like protein, which was identified to play a role in shoot meristem [33] and ovule development[34]. B3 is differentially expressed in anther, and presumably involved in anther development and differentiation [35]. LBD10 encodes a protein that functions in defining the lateral organ boundaries [36]. AT5G58080-ARR18 encodes a type B response regulator that mediates cytokinins signaling transduction in Arabidopsis[37].

Identification of adaptive growth regulators under drought condition

Identification of regulatory genes from water stress data using ASCCA

To test if TF-finder can identify growth regulators from water stress data, we used three files that contained the profiles of 120 target genes, 9 positive TFs (AT3G57600, AT1G75490, AT5G05410-DREB2A, AT2G47190-MYB2, AT1G54160-NF-YA5, AT2G38880-NF-YB1, AT4G27410-RD26, AT1G69600-ZFHD1, and AT4G28110-MYB41) (Additional file 1) and all 1640 TFs detected to be expressed in the water stress data set. The resulting top 70 genes were found to contain 9 novel TFs (Additional file 2) that are supported by existing evidence to be involved in root growth under water stress condition. Again we are not going to elaborate these genes' functions at length. All pieces of evidence that support these genes are positive were shown in Table 1.

ASCCA (Adaptive Sparse Canonical Correlation Analysis)

where each x _i (i = 1, ..., p) is a variable in the set X, and each y _i (j = 1, ..., q) is a variable in the set Y. Then by applying ASCCA on these two sets of data X and Y, we could get a pair of p and q entry weight vectors (canonical vectors), a and b, such that the correlation between the two linear combinations (canonical variates), X a and Y b, is maximized. The canonical vectors a and b are sparse due to many of their entries being zero, which is achieved by introducing L₁ penalties into the criterion that constrains a and b (Witten and Tibshirani, 2010). Specifically, in our study, we only focus on the first component, i.e. first pair of canonical vectors, of the ASCCA solution. To facilitate the ASCCA implementation, Parkhomenko et al. (2009) has developed an iterative algorithm as described below (we assume X and Y are standardized to have columns with zero means and unit variances).

Consider the singular vectors u and v, which are related with the canonical vectors a and b by $a={\sum }_{XX}^{-1/2}u$ and $b={\sum }_{YY}^{-1/2}v$, where Σ _XX and Σ _YY are the variance matrices of X and Y respectively. Given the penalization parameters, λ _u , λ _v and γ, as well as the initial values u⁰ and v⁰, the singular vectors u and v could be approximated iteratively by the following two steps until convergence:

Where X _j and Y _j (the j th subset of the k - fold CV) are the testing sets, and ${\stackrel{\wedge }{a}}^{-j}$ and ${\stackrel{\wedge }{b}}^{-j}$ are the canonical vectors estimated for the training set, in which subset j was removed; Since increasing the penalization parameters decreases the number of non-zero terms in u and v, and for our data u and v would become zero vectors if λ _u and λ _v are greater than 0.4. Consequently, we screen the λ _u and λ _v values from 0 to 0.4 with a step of 0.01, and trace γ from 0 to 2 with a step of 0.1. Finally, 35301 (41*41*21) combinations of the three parameters λ _u , λ _v and γ are examined.

To further evaluate the set of TFs identified by ASCCA, we take the advantage of known positive TFs by examining if the set of identified TFs contains "enough" number of known TFs (Figure 2), which has the similar reasoning as the enrichment test (Rivals et al., 2007) but uses more straightforward and computationally efficient criterion. Denote N as the total number of TFs in X, N _pos as the total number of known positive TFs involved in the same biological process (original input), N _ASCCA as the number of TFs fished out by ASCCA, and N_pos∩ASCCAas the number of known positive TFs that are fished out by ASCCA. Then based on the ratio of positive TFs to total TFs (N _pos / N), the expected number of positive TFs identified by ASCCA is (N _pos / N)*N _ASCCA which is an ideal criterion to be compared with N_pos∩ASCCA, the actual number of positive TFs identified by ASCCA. To make the above criterion more stringent such that only the true significant TF sets being retained, we multiply the expected number of positive TFs by an enrichment factor (EF) which varies from 1~5. That is, if N_pos∩ASCCA> EF * (N _pos / N)*N _ASCCA , the hooked TF set is saved and discarded otherwise. For all the sample results shown in this study, we set EF = 3. In this way, we integrate prior biological knowledge into our mathematical model to deciding if the hooked set of TFs should be retained for further investigation or not.

Cluster analysis

Before applying ASCCA to extract candidate TFs, we applied k-means clustering method (Eisen et al., 1998) to partition the positive target genes into several clusters (Figure 2) and then use each cluster as an input (Y) for ASCCA to bait TFs. The k-means algorithm was selected because: first, target genes in the same cluster are assumed to be co-regulated under the same regulatory machinery and thus each cluster can serve as an ideal bait for ASCCA; second, the result of ASCCA is subject to considerable instability from one input to another, i.e. including or excluding one target gene in Y would possibly result in two quite different sets of TFs. This is not surprising because on one hand the sparse canonical vectors (a and b) are derived from both the greatest correlation between two sets (X and Y) and correlations among variables within each set; while on the other hand, when there is so much information in the datasets (TFs across whole-genome), there exist several alternative solutions that are almost equally good (Waaijenborg et al., 2008). Consequently, because we aimed to identify TFs by the virtue of their true regulatory causality rather than by chances or due to extraneous factors, we performed ASCCA using many target gene clusters and finally averaged the outcomes to minimize the effect of instability.

Because the optimal number of clusters may not exist since the genes involved in different functional domains are co-regulated in varying sizes, we ran cluster analysis several times by varying the number of clusters from a lower to an upper boundary. At the lower boundary, the average number of target genes in each cluster is 20, and at the upper boundary, the average number of genes in each cluster is 4. For instance, given 100 target genes, k-mean clustering analysis is run 17 times with the average size of clusters varying from 4 to 20, and totally ${\displaystyle \sum _{i=4}^{25}{n}_i}$ (n _i represents the number of clusters when the average size of each cluster is i) clusters are processed by ASCCA.

The application of ASCCA on each target gene cluster results in a set of candidate TFs who cooperatively regulate the target genes in this cluster. To extract the truly important TFs from all of the resulting TFs sets, we calculate how many times a TF has been identified by ASCCA. Then the TFs are ranked by the frequency of their occurrence. The more frequent a TF has been identified, the more important is its role in the corresponding biological process. Therefore the list of ranked TFs can provide new hypotheses for further experimental testing. Below is a step-by-step summary of our algorithm:

Step-by-Step Summary of TF-finder

Comparison of TF-finder with ICE

The principle of ICE is based on "guilty by association". It was implemented in such a way that if a candidate gene is associated with a group of positive genes more often than the others [7], this candidate should be selected. In this study, we used a group of positive TFs to judge if a candidate TF is associated multiple times with the members in this group. Due to the presence of multigenic regulation in the gene network, it is usual that transcription regulators controlling the same set of target genes are coordinated or co-expressed. Therefore, we employed Spearman correlation to 'associate' a candidate TF to a number of positive TF. Detail for ICE implementation is described as following. Let Y = {y₁,y₂, ..., y _m } is a set of known positive TFs controlling a biological process, and X = {x₁,x₂, ..., X _n } is a set of TFs across genome with X∩Y = ϕ. A Spearman rank correlation ρ _ij is calculated between any pair of x _i and y _j (i = 1, ..., n, j = 1, ..., m), and X _i and y _j are considered linked when ρ _ij is larger than a pre-specified threshold ρ₀. In our study, we set ρ₀ = 0.6. Then all TFs in X are sorted by the number of links to Y. The genes at the top of the list have more links to Y, and thus are the candidate regulating genes involved in the biological process. Since each selected x _i is located at the "intersection" of multiple elements from Y in a network, we termed this approach as "the intersection of coexpression (ICE)".

Preparation of microarray data sets

Microarray data sets were downloaded from multiple resources. Salt stress experimental data set contains108 chips from 6 experiments (GSE7636, 7639, 7641, 7642, 8787, 5623) and was downloaded from NCBI GEO http://www.ncbi.nlm.nih.gov/geo/. Water stress data sets were downloaded from European Arabidopsis Stock Centre's website http://arabidopsis.info/ and include 62 chips from 3 experiments of AtGenExpress: Stress Treatments (Drought stress) contributed by AtGenExpress Consortium. All data mentioned above are derived from hybridization of Affymetrix 25 k ATH1 microarrays [42]. The original CEL files were processed by the robust multiarray analysis (RMA) [43] algorithm using the Bioconductor package. For quality control we used methods that were previously described [44]

Availability of software package

The ASCCA package was written in R. A wrapper for calling ASCCA, and a number of parsers were written in Perl. To facilitate use of this package, we release for public use the original codes rather than executables. The users need to use an Unix/Linux environment where R and Eisen's k-means clustering package [45] are installed. Installation of Perl is not necessary because it is usually carried by the Linux/Unix operating system. Interested users can receive the package by sending email to: hairong@mtu.edu.