Model-based comparative prediction of transcription-factor binding motifs in anabolic responses in bone.

Understanding the regulatory mechanism that controls the alteration of global gene expression patterns continues to be a challenging task in computational biology. We previously developed an ant algorithm, a biologically-inspired computational technique for microarray data, and predicted putative transcription-factor binding motifs (TFBMs) through mimicking interactive behaviors of natural ants. Here we extended the algorithm into a set of web-based software, Ant Modeler, and applied it to investigate the transcriptional mechanism underlying bone formation. Mechanical loading and administration of bone morphogenic proteins (BMPs) are two known treatments to strengthen bone. We addressed a question: Is there any TFBM that stimulates both "anabolic responses of mechanical loading" and "BMP-mediated osteogenic signaling"? Although there is no significant overlap among genes in the two responses, a comparative model-based analysis suggests that the two independent osteogenic processes employ common TFBMs, such as a stress responsive element and a motif for peroxisome proliferator-activated receptor (PPAR). The post-modeling in vitro analysis using mouse osteoblast cells supported involvements of the predicted TFBMs such as PPAR, Ikaros 3, and LMO2 in response to mechanical loading. Taken together, the results would be useful to derive a set of testable hypotheses and examine the role of specific regulators in complex transcriptional control of bone formation.

Introduction

One of the goals of systems biology is to raise testable hypotheses based on the large amount of experimental data generated from high-throughput technologies such as microarray. With the availability of genomic DNA sequences, computational tools have been developed to predict the cis-acting elements—transcription-factor binding motifs (TFBMs)—in the promoter and enhancer regions that cause the global gene expression patterns to emerge 1., 2..

We previously reported an ant algorithm to predict functional TFBMs from the microarray data for chondrogenesis 3., 4.. The ant algorithm 5., 6. is a branch of swarm intelligence techniques (, and it was inspired by the foraging behavior of ant colonies of Linepithema humile that communicate each other through deposition of a chemical named pheromone. The concentration of pheromone on trails is considered “distributed information” and is constantly evaporated and re-deposited by ants to reflect their experience while searching for food. In ant algorithms, artificial pheromone, which is given by ant-like agents to a potential solution, is used as a clue to obtain sub-optimal solutions in combinatorial problems. The identification of TFBMs is formulated as a combinatorial optimization problem with the assumption that the mRNA level is determined by multiple transcription factors and each transcription factor controls multiple genes. In the current study, each ant represents a particular selection of putative TFBMs and deposits a varying amount of artificial pheromone based on how well the occurrence of this choice of motifs correlates with the observed mRNA levels. The algorithm is suited to identify a heterogeneous set of motifs whose length varies from 4 bp to more than 20 bp.

Here we developed a set of web-based software, Ant Modeler, in order to facilitate model-based comparative prediction of TFBMs. In many biological systems, more than one treatment can induce a common outcome such as cellular proliferation and differentiation as well as apoptosis. In order to identify TFBMs in more than one treatment, a model-based comparative analysis is useful. Among many TFBMs involved in those treatments, the analysis allows us to distinguish a group of TFBMs that are commonly utilized from other TFBMs that are specific to each treatment. Ant Modeler is built to provide a user-friendly platform through which putative TFBMs are identified from microarray data in a public or local domain. It receives a list of genes and their fold changes in response to any treatment, and generates a report for a group of putative TFBMs, where a concentration of pheromone deposited to TFBMs by artificial ants represents a predicted contribution of each TFBM to the observed expression profile.

In this study, we applied Ant Modeler to predict TFBMs involved in two independent treatments for bone formation: mechanical loading ( and administration of bone morphogenic proteins (BMPs) (. Bone is a dynamic system that is constantly destroyed by osteoclasts and rebuilt by osteoblasts 10., 11.. It is therefore important, particularly in the aging population, to enhance bone formation and prevent bone loss. Although the two treatments are known to stimulate bone formation, their patterns of gene expression are significantly different. An intriguing question is whether the two treatments share any common TFBMs essential for bone formation. With Ant Modeler, a comparative analysis was conducted using two sets of microarray data in the public domain for “in vivo mechanical loading” and “in vitro administration of BMPs”.

Results and Discussion

Two datasets were employed from the microarray data in response to mechanical loading published by Xing et al. ( and BMP administration by Peng et al. (. We included all genes that were biologically identified with clear annotation (69 genes in mechanical loading and 53 genes in BMP administration). To our surprise, only two genes (lysyl oxidase and neuropilin) are common to the two datasets. We searched TFBM candidates using 1,000-bp DNA sequences flanking to the 5’-end of the genes. The prediction of TFBMs was conducted using 100 artificial ants with 1,000 iterations (see Materials and Methods for details).

Prediction of 5-bp TFBM candidates

We first applied Ant Modeler to evaluate 5-bp TFBM candidates in the two datasets. There are 512 candidates without considering the polarity of DNA sequences. The plot of pheromone distributions for 512 TFBM candidates exhibits multiple peaks (Figure 1). Among them, two local peaks (“AAGGG” and “GGGCA”) appeared to be common to the two examples. They are part of the consensus sequences of stress responsive element (STRE) ( and peroxisome proliferator-activated receptor (PPAR) (.

Fig. 1

Spectrum of pheromone concentrations for 5-bp TFBM candidates (512 in total starting from “AAAAA” along the x-axis). A. Spectrum of pheromone concentrations for the dataset linked to mechanical loading. B. Spectrum of pheromone concentrations for the dataset linked to BMP administration

Linkage to biologically known consensus sequences using 6-bp TFBM candidates

The results based on 5-bp TFBM candidates indicate a common usage of some TFBMs in the two responses, but most biologically known motifs are longer than 5 bp. In order to investigate the role of regulatory sequences longer than 5 bp, we selected 6-bp TFBM candidates and mapped them to TRANSFAC database version 9.2 (. Firstly, using Ant Modeler, the 6-bp DNA sequences whose pheromone concentration is significantly higher [> 3 standard deviation (s.d.)] than the average were chosen as potential core TFBMs (Figure 2). There were 29 and 25 outliers in the prediction for mechanical loading and BMP administration, respectively. Secondly, these outliers were mapped to biologically known TFBMs in TRANSFAC database and a cumulative score for each of the known TFBMs was defined: where K = the number of 6-bp outliers that match to the known TFBM of interest; S = similarity score of the predicted motif to the position-specific scoring matrix of the biologically known binding sites in the TRANSFAC database (; L = total length (in bp) of the known TFBM of interest; and L = length (in bp) covered by 6-bp outliers. Besides STRE/Nrf2 and PPAR, which were predicted in the analysis with 5-bp candidates, the linkage analysis through mapping to TRANSFAC also identified other known motifs and transcription regulators, such as HEN1 (basic helix-loop-helix protein), Ikaros 3 (zinc finger protein), Helios A (Ikaros-related protein), NRSF (neuron-restrictive silencer factor), major T-antigen binding site, LMO2 (LIM domain transcription regulator 2), GCNF (germ cell nuclear factor), STAT6 (signal transducer and activator of transcription 6), ER (estrogen receptor), and ERRα (estrogen-related receptor α) (Table 1). Since the TRANSFAC database often includes a few stretches of nonspecific nucleotides as part of consensus sequences, the predicted percent coverage in Table 1 is underestimated.

Fig. 2

Distribution of the total number for each of the 6-bp TFBM candidates (2,080 in total) to be selected by Ant Modeler. A. Distribution for the dataset linked to mechanical loading. B. Distribution for the dataset linked to BMP administration.

Table 1

TFBM candidates predicted for mechanical loading and BMP administration

Motif*	Total score	Treatment#	Score	Coverage	Sequence
HEN1	26.5	ML	19.8	73% (16/22)	nngGGNCGCAGCTGCGNCCcnn
		BMP	6.6	36% (8/22)	nngggncGCAGCTGCgncccnn
PPAR	23.5	ML	14.9	62% (13/21)	nnwgRGGTCAAAGGTCAnnnn
		BMP	8.65	33% (7/21)	nnWGRGGTCaaaggtcannnn
COUP	23.4	ML	15.2	100% (13/13)	TGACCTTTGACCC
		BMP	8.2	54% (7/13)	tGACCTTTgaccc
HNF4	20.4	ML	15.9	87% (13/15)	nRGGNCAAAGGTCAn
		BMP	4.5	40% (6/15)	NRGGNCaaaggtcan
LXR	8.7	ML	9.0	71% (12/17)	YGAMCTnnasTRACCYn
		BMP	8.7	71% (12/17)	yGAMCTNnastRACCYN
Ikaros 3	22.4	ML	4.5	40% (6/15)	NRGGNCaaaggtcan
		BMP	12.1	85% (11/13)	tNYTGGGAATACc
Helios A	17.3	ML	12.0	82% (9/11)	nNTWGGGANNn
		BMP	5.3	55% (6/11)	nNTWGGGannn
NRSF	22.0	ML	8.1	57% (12/21)	ttcagCACCACGGACAGmgcc
		BMP	13.9	71% (15/21)	ttcaGCACCACGGACAGMGcc
major T-antigen	21.1	ML	16.8	84% (16/19)	GGGAGGCAGAGGCAGGygg
		BMP	4.3	32% (6/19)	gggagGCAGAGgcaggygg
LMO2	18.2	ML	11.4	75% (9/12)	cNNCAGGTGBnn
		BMP	6.8	50% (6/12)	cnncAGGTGBnn
GCNF	17.8	ML	11.1	67% (12/18)	ntcaAGKTCAAGKTCAnn
		BMP	6.7	44% (8/18)	ntcAAGKTCAAgktcann
STAT6	16.8	ML	10.4	88% (7/8)	NNYTTCCy
		BMP	6.4	75% (6/8)	NNYTTCcy
ER	14.1	ML	7.6	55% (6/11)	nAGGTCAnnny
		BMP	6.5	55% (6/11)	NAGGTCannny
ERRα	12.9	ML	6.9	43% (6/14)	nnntnaAGGTCAnn
		BMP	6.0	43% (6/14)	nnntnAAGGTCann
STRE	13.6	ML	6.3	75% (6/8)	TMAGGGgn
		BMP	7.2	75% (6/8)	TMAGGGgn

HEN1: helix-loop-helix protein 1; PPAR: peroxisome proliferator-activated receptor; COUP: chicken ovalbumin upstream promoter; HNF4: hepatocyte nuclear factor 4; LXR: liver X receptor; NRSF: neuron-restrictive silencer factor; LMO2: LIM domain transcription regulator 2; GCNF: germ cell nuclear factor; STAT6: signal transducer and activator of transcription 6; ER: estrogen receptor; ERRα: estrogen related receptor α; STRE: stress responsive element.

ML: mechanical loading; BMP: BMP administration.

Biological considerations

Out of the top ten candidates in Table 1, four motifs, namely PPAR, STAT6, ER, and STRE/Nrf2, are reported to be related to bone formation or resorption. The suppression of PPAR is shown to stimulate differentiation of osteoblast cells through Wnt signaling pathway (. STAT6 is known to inhibit osteoclast differentiation and thereby bone resorption (. ERs are demonstrated to mediate skeletal growth and differentiation and their effects are modulated by genders 17., 18.. STRE/Nrf2 is reported to interact with activating transcription factor 4, which is an essential transcription factor in bone formation (. Six other motifs, including HEN1 (, Ikaros 3 and Helios A (, NRSF (, Major T-antigen (, LMO2 (, and GCNF (, have not been shown with any linkage to bone formation or resorption. All these six motifs, however, play critical roles in the processes of development and regeneration in other tissues 26., 27., 28., 29., 30., 31..

In vitro biological evaluation in response to mechanical shear stress supports some of the prediction by Ant Modeler (Figure 3). For instance, potential involvements of PPAR, Ikaros 3, and LMO2 are suggested (Table 1). Real-time PCR revealed that their mRNA levels in mouse osteoblasts were significantly altered by 1 h exposure to mechanical loading both at 10 dyn/cm2 and 20 dyn/cm2. Five hours after loading, the level was reduced to 0.22±0.04 and 0.11±0.04 (mean±s.d.) (PPAR), 0.33±0.04 and 0.24±0.04 (Ikaros 3), and 0.35±0.06 and 0.14±0.02 (LMO2) for 10 dyn/cm2 and 20 dyn/cm2, respectively. In summary, the model-based analysis herein would be useful to generate a series of new hypotheses that can be experimentally tested to unravel a novel pathway to stimulate bone formation.

Fig. 3

mRNA expression levels of PPAR, Ikaros 3, and LMO2 in mouse MC3T3 osteoblast cells (C4 clone) in response to fluid shear for 1 h. A. Responses to 10 dyn/cm2 fluid shear. B. Responses to 20 dyn/cm2 fluid shear.

Materials and Methods

Two datasets for bone formation

The two datasets analyzed in the present study were employed from the microarray data in response to mechanical loading published by Xing et al. ( and BMP administration by Peng et al. ( (Figure 4). In the first dataset, mechanical loads were applied to the right tibia of the mice with the left tibia used as unloaded control. Agilent mouse development oligonucleotide microarray slides (containing approximately 22,000 markers) were used, and the results were presented in the paper as well as in the supplementary data of Xing et al. (. In the second dataset, mouse pluripotent mesenchymal precursor line C2C12 was infected with three viruses consisting of BMP-2, BMP-6, and BMP-9 genes and the effects of these BMPs were assayed using the Affymetrix’s mouse genechip U74Av2 (containing approximately 12,000 markers). The genes with significantly altered mRNA levels were presented (.

Fig. 4

Flow chart for the application of Ant Modeler. Two mRNA expression datasets (mechanical loading and BMP administration) were used to predict TFBMs involved in bone formation.

Ant Modeler

The web-based software Ant Modeler was constructed in combination of the script language Perl (version 5.8.5), rational database mySQL (version 4.1.20), and statistical computing package R (version 2.3.1) (Table 2). It is applicable for three mammalian species (Homo sapiens, Mus musculus, and Rattus norvegicus) commonly employed in expression analyses (Figure 5). The start page requests input in a form of a tab-delimited or comma-delimited file. The file documents a list of genes whose expression levels are altered with their fold changes, the species of the sample of interest, and the length of promoter sequences (upstream of the transcription starting site) for prediction of TFBMs. In addition, a set of parameters are needed to be defined for the prediction, including the number of artificial ants, number of iterations, pheromone preference factor (ɛ), pheromone evaporation factor (δ), and power factor for error evaluation (ɑ). For the detailed definition of each parameter, please refer to Liu et al. (. Default parameter values, which seem appropriate for most applications, were provided.

Table 2

Configuration of three subsystems in the web interface

Subsystem	Function	Software
User interface	Input microarray data and set up parameters	Perl
Database annotation*	Retrieve DNA sequences of the genes of interest	mySQL
Computation	Predict a set of TFBMs	R

The following genes, which were included in the two original datasets, were not included in the present analysis. Dataset 1: BC015839, BC030010, BG071710, BG071952, BG072471, NM021584, NM008788, and 10 genes starting with “XM”; Dataset 2: AI851750, AF004874, AV109962, U68267, AV093331, AV359510, and L10076.

Fig. 5

Web-based interface of Ant Modeler. The required inputs include a data file (gene accession number and fold change), an organism of interest (Mus musculus, Rattus norvegicus, or Homo sapiens), and 9 parameters to run the ant algorithm.

Formulation with multiple TFBMs

The fold change values of gene expression were expressed in a piecewise logarithmic model (:where Z = logarithmic ratios of differential gene expression levels, h in H× = the j-th TFBM in the regulatory region of the i-th gene, X = functional levels of TFBMs, m = number of TFBMs, and n = number of genes. For each gene, we searched through the contig file and identified the position of the transcription start site. A 1,000-bp region upstream of the transcription-starting site was used as the regulatory region. To estimate the functional level of each motif, the least-square method was employed:

Ant Modeler was used to minimize the cost function, which was defined as: (differences between the experimental and equation-predicted expression levels), ɑ > 1.

Strategy for pheromone deposition, pheromone-guided selection, and pheromone evaporation

Pheromone deposition

Each ant was assigned a set of m random TFBM candidates. The cost function for each ant was evaluated, and the amount of pheromone deposited on each candidate motif was defined as the inverse of the cost function.

Pheromone-guided selection

The candidate selection was based on the pheromone concentration of each potential TFBM. The probability of selecting potential motif j in the (i+1) iteration step was defined as:where ε = pheromone preference factor (ε > 0), F, = cumulative pheromone concentration of motif candidate j at iteration I, and M = total number of TFBM candidates. The variable ɛ controls weight to the current pheromone concentration. Greater values of ɛ will produce results with heavy preference to pheromones, while a value of ɛ = 0 would conduct a random selection without any preference to pheromones.

Pheromone evaporation

At each iteration, the pheromone concentration was updated, including a process of evaporation with δ as the pheromone evaporation factor (0 ≤ δ ≤ 1). Note that with δ = 0 the pheromone concentration is conserved without evaporation, whereas with δ = 1 the previous pheromone concentration is totally lost.

Biological evaluation in response to mechanical loading using mouse osteoblast cells

MC3T3 osteoblast cells (C4 clone) were grown on a glass slide coated with 40 μg/mL type I collagen (BD Biosciences) in ɑMEM containing 10% FBS and antibiotics. Prior to mechanical loading, the cells were incubated in the medium containing 1% FBS for 24 h. They were then exposed to 1 h flow shear at intensity of 10 dyn/cm2 or 20 dyn/cm2 in a Streamer Gold flow device (Flexcell International) (. Total RNA was extracted 1 h, 3 h, and 5 h after the onset of flow shear using an RNeasy plus mini kit (Qiagen). Reverse transcription was conducted, and real-time PCR was performed using ABI 7500 with SYBR green PCR kits (Applied Biosystems). The PCR primers were PPARγ (5’-GGAAAGACAACGGACAAATCA-3’ and 5’-TACGGATCGAAACTGGCAC-3’), Ikaros 3 (5’-ATGGATGTCGATGAGGGTCAAG-3’ and 5’-TTAGCTCAGGTGGTAACGATGC-3’), LMO2 (5’-TCAGCTGTCACCTCTGTGG-3’ and 5’-CACCCGCATCGTCATCTC-3’), and GAPDH (5’-TGCACCACCAACTGCTTAG-3’ and 5’-GGATGCAGGGATGATGTTC-3’), where GAPDH was used for internal control.

Authors’ contributions

WX and SM conducted microarray experiments. HY and YL developed the ant algorithm and drafted the manuscript. AC and GW coded and run the programs. KH conducted biological interpretations and shared writing the manuscript. All authors read and approved the final manuscript.

Competing interests

The authors have declared that no competing interests exist.