Hierarchical Modularity in ERα Transcriptional Network Is Associated with Distinct Functions and Implicates Clinical Outcomes

Recent genome-wide profiling reveals highly complex regulation networks among ERα and its targets. We integrated estrogen (E2)-stimulated time-series ERα ChIP-seq and gene expression data to identify the ERα-centered transcription factor (TF) hubs and their target genes, and inferred the time-variant hierarchical network structures using a Bayesian multivariate modeling approach. With its recurrent motif patterns, we determined three embedded regulatory modules from the ERα core transcriptional network. The GO analyses revealed the distinct biological function associated with each of three embedded modules. The survival analysis showed the genes in each module were able to render a significant survival correlation in breast cancer patient cohorts. In summary, our Bayesian statistical modeling and modularity analysis not only reveals the dynamic properties of the ERα-centered regulatory network and associated distinct biological functions, but also provides a reliable and effective genomic analytical approach for the analysis of dynamic regulatory network for any given TF.

MCF7 cells were treated with DMSO (as 0 hour time point) or E2 (108 mol/L) for 1, 4 and 24 hours. 5 x 10 7 cells were cross-linked with 1% formaldehyde for 10 min, at which point 0.125 M glycine was used to stop the crosslinking. In brief, after crosslinking, cells were treated by lysis buffers and sonicated to fragment the chromatin to a size range of 500bp-1kb. Chromatin fragments were then immunoprecipitated with 10ug of antibody/magnetic beads. The antibodies against ERα were purchased from Santa Cruz Biotechnology (Santa Cruz, sc-8005 X). After immunoprecipitation, washing, and elution, ChIP DNA was purified by phenol:chloroform:isoamyl alcohol and solubilized in 70 μl of water. The ChIP DNA sample was run in 12% PAGE and the 100-300bp DNA fraction was excised and eluted from the gel slice overnight at 4 °C in 300 μl of elution buffer (5:1, LoTE buffer (3 mM Tris-HCl, pH 7.5, 0.2 mM EDTA)-7.5 M ammonium acetate) and was purified using a QIAquick purification kit (Qiagen, Cat#28104). The library was constructed using Illumina genomic DNA prep kit by following its protocol (Illumina, cat# FC-102-1002), clusters were generated on the Illumina cluster station (Illumina, cat# FC-103-1002), DNA samples (20 nM per sample) quantified by an Agilent Bioanalyzer, were loaded onto Illumina Genome Analyzer IIx (GAIIx) for sequencing according to the manufacturer's protocol. Reads generated from the Illumina GAIIx pipeline were aligned to the Human Genome Assembly (NCBI build 36.1/hg18) using ELAND algorithm.

ChIP-seq data peak-calling analysis:
After mapped to the corresponding reference genome (HG18), the processed data sets need peak-calling and further identification of ER binding sites, currently there are quite a few programs and scripts available for such tasks, such as FindPeaks 1 , Quest 2 , and MACS 3 , etc. Our ChIP-seq data sets were performed peak-calling by the in-house program, wBELT, developed in our laboratory 4 .
And more details on performance comparison and implementation can be retrieved from the corresponding references.

Bayesian multivariate statistical modeling of genetic transcription rates:
In most cases, Bayesian statistical models are concerned with learning the parameter set θ=(θ 1 , …, θ l ) of the dimension l, containing uncertain quantities (fixed and random effects), hierarchical parameters (hyperparameters), unobserved and/or latent variables, and even missing data 5 .
In the Bayesian statistical modeling schema, prior information about the model parameter is denoted by the probability density function p(θ), the likelihood function is represented as p(X|θ), and the inference purpose is to derive the posterior density function p(θ|X). According to the Bayes theorem, the general inference can be formulated as, For the multivariate case, p(X) can be specified by the decomposition, which can be regarded as a normalized constant by integrating over all values of θ in the product p(X|θ)p(θ). Thus, the equation above can be formulated as, The strength of the observation data and corresponding prior knowledge influence those diverse weights on the beliefs inferred from the multiple sources.
Due to a relatively small sample size of E2-stimulated time-series gene expression data which contain no knowledge of transcription factors, binding information and direct target genes, the conventional Bayesian modeling has its limitation. In order to infer ERα-centered regulatory network, we integrate the time-series E2-stimulated ERα ChIP-seq data, where it can detect transcription factors and hubs, and facilitate the further reverse-engineering of the regulatory network by means of inferring parameters in the Bayesian statistical framework.
Herein we propose a Bayesian multivariate statistical approach for modeling the time-variant ERα transcriptional regulatory network. The basic model framework is illustrated as follows, , ( ) ( ) , 1,..., , 1,..., where ( ) i y t  denotes the ith gene's transcription rate, x j (t) for the jth gene's expression level at the investigated time, α ij for the corresponding regulatory argument or strength of the jth gene which has any possible transcription regulatory activity on the ith gene, and ε represents the potential stochastic effects during the transcription regulatory process, which normally follows a normal distribution, i.e.
Thus for a genetic regulatory network containing M transcription factors at T time points, the above equation can be organized as, and the error term  is specified as an independent and normally distributed random vector with Ndimensional zero mean and a N N  covariance  . And the multivariate N-dimensional normal distribution for errors is represented as, where ε i is the N-dimensional error vector. And from the multivariate normal error specification, the observation vector also follows a multivariate normally distribution, denoted as, then with extension to matrix model representation, the above equation can be formulated as, where the trace operator tr defines the sum of the diagonal entries of its matrix argument. Here, we denote the conditional prior distribution for A as p(A|Σ), which follows normal distribution, denoted as, And the prior distribution for the error covariance matrix Σ as p(Σ), and normally the inverse Wishart distribution (IWD) 6,7 is considered as the conjugate prior for the covariance matrix Σ of the multivariate normal distribution, depicted as, where M, d and S are hyperparameters. According to the Bayes theorem, the joint posterior distribution of the coefficient matrix can be denoted as, For the multivariate modeling, inference of the coefficient matrix contains two sub-problems, i.e. to determine the parameters' posterior distributions for the model, and then to apply the marginal posterior mean estimation for the parameter inference. Thus to determine the marginal posterior distributions for coefficient matrix A, the joint posterior distribution can be marginalized with integration on Σ, illustrated as, , X 0 is the prior input matrix, δ=(n-p+q)/2>0 and Q is a high-order matrix independent of A. The integration for the posterior distribution is taken as the same form as an inverse Wishart distribution in Bayesian multivariate statistical estimation 6,7 . And denote where A 0 is the prior mean matrix. Thus the maximum of posterior coefficient mean estimation A should satisfy the following partial derivative, Meanwhile, due to the unobservability of the transcriptional rate, Y  , hence the argument can be approximated as follows, where Δt denotes the sampling interval. Thus the first-order differential of transcription rate can be replaced with the observable vector, Y(t).
Furthermore, due to proportionality characteristics of coefficient matrix inferred by Bayesian statistical analysis, we normalize those coefficients by scaling them with the range within -1 to 1.

Identification of recurrent regulatory motif patterns:
From the those regulatory modularity analysis, we mainly find the four recurrent regulatory motif patterns, i.e. self-loop, interactive loop, feed-forward loop, feedback loop and bi-span loop ( Figure S7). Actually, such motif patterns implement specific biological functions in genetic regulatory processes, i.e. (i) Self-regulation: the regulatory activities include simple self-activation or selfinhibition on those genes under investigation, see Figure S7 (a) self-loop; (ii) Cooperative activation and inhibition regulation: more than two regulators cooperatively regulate (activate or inhibit) one or several target genes, see Figure S7 (c) feed-forward loop and (d) bi-span; (iii) Feedback mechanism: one or more feedback loops exist in source regulators and target genes.
Such feedback loops directly reflect the regulatory effects on target genes to source regulators, thus simultaneously the source regulators self-tune their regulatory activities. Normally there exists a dynamic equilibrium between source and target genes' activities, see Figure S7 (b) interactive loop.

The Kaplan-Meier survival analysis based on diver patient cohorts:
For the clinical outcome analysis, the gene signatures are the specific genes that characterized the diverse clustered patient subgroups.
In determining the gene signature and clinical outcome analysis, we firstly adopted the k-means clustering approach to analyze the three patient cohorts containing 337, 251 and 137 patients, respectively. Since initially we have no concrete group information about the signatures, actually the k-means clustering here is an unsupervised approach.    :   Table S1. The network elements' SNR statistics (dB) in Module I.