Predicting gene expression in T cell differentiation from histone modifications and transcription factor binding affinities by linear mixture models

Regulatory signals predicts expression

As a first step, we want to determine which set of regulatory signals, X, can explain the observed gene expression data, Y , best. To this end, as a first step we supply our algorithm with a matrix X containing only predicted TF binding affinities, only histone modification data or both sets of regulatory signals and assess how well the resulting regression models can capture the data. As measure of quality for the different models we thereby compute the mean square error between the predicted gene expression values and the actual measurements (see Methods for details).

Predicting the gene expression from the four T-cell types based on only histone modification data and by means of only a single regression model yields MSEs of about 0.5 for HM and HM+TF on all data sets (see red bars in Fig. 2). A mixture of two regression models further reduces the MSEs to an average value 0.25 across all cell types. In all scenarios, the difference of MSE between one and two models were statistically relevant (t-test p-value < 0.01) indicating the advantage of using mixtures to predict expression. The model selection procedure (see Methods) indicates that the data is optimally explained by the combination of 2-4 regression models (see Fig. 2) and that gene expression data can be well predicted based on histone modification data alone.

In contrast, using a single regression model to predict gene expression data based on TF binding affinities alone yields considerably larger MSEs across all cell types (average MSE = 2, see blue bars in Fig. 2). Interestingly, supplying our algorithm with the combined data from both histone modifications and TF binding affinities yields MSEs similar to the ones obtained with only histone modification data alone (see Fig. 2). This indicates that the utilized histone modification and TF binding data cover rather redundant than complementary information about gene expression. As histone modifications and HM+TF affinities yield the solutions with the lowest MSEs, in the following we will continue the analysis with the models based on these data sets.

Control of Th1 gene expression

Having established that histone modification data together with a mixture of two regression models yields the most significant results we now investigate which type of modifications contribute most strongly to these models. We thereby restrict our analysis to the data from Th1 cells and the corresponding regulatory signals which obtain the largest absolute regression coefficients (results from other cell types closely resemble those from Th1 cells, see Additional File 1 for details).

For model 1, which explains the expression pattern of the most highly as well as moderately expressed genes, the histone modifications with largest influence are H3K4me3 and H3K27me3, with regression coefficients of +0.7975 and -0.4533, respectively. As shown in the top part of Fig. 3 these two modification form a gradient with H3K4me3 being most frequently found in the highly expressed genes while being absent in moderately to lowly expressed genes. In contrast, H3K27me3 is consistently detected in promoters of lowly expressed genes but appears weaker or even absent in promoters of highly expressed genes.

For model 2, which explains the transcriptional activity of a small subset of highly expressed as well as most of the lowly expressed genes, we again find H3K4me3 to have the strongest positive regression coefficient (b = 0.71). This is reflected by a strong association of this modification with the most highly expressed genes of this set (see Fig. 3). In contrast, H3K27me3 obtains a regression coefficient of close to zero (b = 0.03) in this model as this modification appears with the same intensity in nearly all genes assigned to model 2.

An alternative view of this results is presented at Fig. 4. There, we have the interpolated values of the histone modifications against the gene expression, the linear models for each component and the resulting mixture model. Clearly, dependence of H3K27me3 on gene expression is not linear, as low expressed genes all present a high presence of histone modifications. This non-linearity is captured by the mixture model (red line), and explains the lowest MSE errors obtained when more than one linear models is applied. To see whether the influence of TFs may contribute to this effect we next look in detail at the results obtained from combined histone and TF data together with a mixture of three correlation models. As shown in Fig. 5), we see similar results in respect to the histone markers: H3K4me3 as enhancer and H3K27me3 as inhibitor of expression for genes with high expression and H3K4me3 as enhancer for genes with low expression. Moreover, only for the genes with high expression, there are some TFs (Pax5, Stat5, Meis/Hox, Iscbp) promoting gene expression and a TF (MyB) inhibiting expression. For all TFs, regression coefficients were in the range of 0.1 to 0.15 (see Additional File 1 for additional results). For genes with low expression, we found no relation between the TFs binding affinities and expression. Th1 cells are known to be regulated by T-bet and Stat4 [1]. While our study lacked the PFM of T-bet, it listed the closely related Stat5 as a positive regulator of the genes with high expression. In relation to factor related to inhibition, there has been a recent implication of c-MyB to bind to H3 histone tails and to promote histone acetylations in Humans [22]. These results indicate a putative role of the MyB in down-regulating the expression of genes during CD4+ T differentiation by promotion of epigenetic changes. However, further acytilation modification data would be required for a better characterizaion of the role of this factor.

Comparison with previous studies

Several computational biology methods have been previously proposed for the use of linear models for predicting gene expression in the context transcription factor binding [10, 11] or histone modifications [12]. In all cases, distinct datasets were used and results are not directly comparable. In relation to [12], the analysis were based on Human naive CD4+ T cells and included 38 histone modifications. Their predicted model obtained a correlation coefficient of 0.72 on genes with low CpG content with HM H3K4me3 and H3K79me1, while our method had a coefficient at the range 0.64 – 0.68 for one model and 0.85 – 0.87 for two models for H3K4me3 and H3K27me3 data. The increase of the correlation coefficient from single linear models to two linear models is an indication that all these approaches would profit from the use of the mixture of linear regressions framework.

Mixture of linear regressions

In the following we want to model the observed expression level of all N genes, using different linear combinations of the M different regulatory signals associated with the promoters (i.e. binding affinities for various TFs and different histone modifications). To this end let y _i be the gene expression level of gene i (the dependent variable) and x _i be a corresponding vector of M regulatory signals (the regressor variables). The single linear regression model is then defined as

y _i =b₀+ x _i B ^T + ∊_i, (1)

where B is a vector (b₁, …,b _M ) representing regression coefficients and ∊_i is an error term. For mathematical convenience, we redefine the vector with the regressor variables to be x _i = (1,x _{i 1}, …,x _iM ) and include the bias parameter b₀ in the beginning of B, that is B = (b₀, b₁,…, b _M ). Assuming the error ∊ follows a Normal distribution with standard deviation σ², the linear regression model has the following distribution

ℙ(y _i |x _i , B,σ² ) = N(y _i |x _i B ^T ,σ² ). (2)

A mixture of linear regression models is defined as a convex summation of K distributions

(3)

where Π = (π₁, …,π _K ) are the mixture coefficients, which respect π _k ≥ 0 and , and Θ are the model parameters (Π, B₁,…, B _K , , …, ).

For a given data X and Y, where X is a set on N observations x _i and Y a vector with N observations y _i , the mixture of linear regression models can be estimated with the Expectation-Maximization algorithm [14, 24]. We resort to Maximum-a-posteriori (MAP) estimates of the parameters, as described in the next section, to avoid over-fitting [25]. The EM works by finding estimates Θ maximizing the posterior distribution over the data X and Y

ℙ(Θ|X, Y, Z) ≈ ℙ(Y, Z|X,Θ)ℙ(Θ) (4)

where Z is the vector of hidden variables with z _i ∈ {1,…, K} indicating which linear model an observation i belongs to. ℙ(Θ) is the prior distribution over the model parameters (see next section for the definition of the prior). ℙ(Y, Z|X, Θ) is the complete data likelihood and is given by:

(5)

where r _ik is the posterior probability (or responsibility) [25] that observation i belongs to the linear model k and is given by:

(6)

For further details on mixture models we refer the reader to [25].

The EM algorithm works by iteratively estimating the model assignments (r _ik ) and the model parameters Θ until some convergence criteria is reached. In the context of the mixture of linear regression models, we need estimates of the linear regression parameters (B _k , ) for a particular model k, and all other parameters (r _ik ,Π) follow the usual EM algorithm [25].

Once the mixture model is estimated, the predicted value ŷ _i for a particular regressor observation x _i is given by

(7)

That is, the linear regression prediction is a mixture of the predictions of each individual component times the posterior probability of the observation i to belong to the model k. In our particular application problem, we are interested in estimating the models which corresponds to an unsupervised learning problem, that is, the coefficients indicating whether a regulatory signal plays an important repressive or activating role. The predictions ŷ can thereby be used for evaluating the fit of our model. In cases where one wants estimate the expression level of genes, that is, estimation of ŷ (supervised learning problem), the above equation should not be used, as the posterior probabilities are based on the response variable y, which is usually unknown in a predictive scenario. In such a context, methods for combinations of predictors, such as [26], are required.

Bayesian linear regression estimates

We resort to Bayesian approach for obtaining MAP estimates of the linear regression models as proposed in [27]. Therefore, we avoid problems related to over-fitting which usually occur with the EM algorithm and mixture models [25]. More formally, the prior distribution in Eq. 4 can be decomposed as

(8)

We use the following conjugate prior for the regression coefficient B _k

ℙ(B _k ) = N(B _k |0, β _k I), (9)

where 0 is a vector with M zeros, I is a M x M identity matrix and β _k is the hyper-parameter.

Let r _k be an N dimensional vector (r_{1 k} , …,r _Nk ) containing the posterior probabilities of the observations belonging to model k and let W _k = diag(r _k ), then the estimates from model k maximizing Eq. 4 are defined as

(10)

with

(11)

From Eq. 10, we can see that β _k works by shrinking the regression coefficients. Small β _k imposes a higher shrinkage on the regression coefficients. Furthermore, for β _k → ∞ we have a non-informative prior and the regression coefficients are the maximum likelihood estimates.

We estimate the hyper-parameter β _k in an Empirical Bayes approach with

(12)

where

(13)

and λ _j is the jth eigenvalue of the PCA decomposition of matrix (see [27] for details).

Note that β _k requires the definition of B _k , which in our context is taken from the previous iteration of the EM algorithm.

For the mixture mixing coefficients, we use a symmetric Dirichlet distribution as prior

ℙ(Π) = Dirichlet(Π|α), (14)

where α is the hyper-parameter. Hence, the mixing coefficients estimates used by the EM algorithm are

(15)

We use a prior of α = 2, which avoids models with a low number of observations assigned to it.

Transcription factor affinity

TF binding motifs are traditionally described in the form of position frequency matrices (PFMs). PFMs show how often a certain base occurs at a given position in the alignments of known binding sites of the TF. To predict the binding strength of a given TF to a promoter sequences we utilize the TRAP method [17]. In contrast to motif matching algorithms which make a binary distinction between binding sites and non-binding sites, TRAP avoids this artificial separation and instead computes the probability of a TF to bind site i in the sequence using the following equation

(16)

where δE _i (λ) is the energy difference between the state in which the factor is bound to site i and the state in which the factor is bound to its consensus site. This so called mismatch energy is scaled by a parameter λ which was previously determined to have an optimal value of 0.7 [17]. The second transcription factor dependent parameter R₀ determines both, the binding energy between the factor and its consensus site as well as the TF concentration. R₀ is derived for each PFM individually as

R₀ = exp(0.6 • W – 6), (17)

where W is the number of columns in the PFM with information content exceeding 0.1 bits. Matrix positions which fall below this entropy cutoff also do not contribute to the mismatch energy in Eq. 16. The nucleotide dependent mismatch energies for each site in the promoter sequence are computed by

(18)

where v _{i , max} is the frequency of the consensus base at position i in the PFM and v _{i , α} , is the frequency of the observed base α at position i in the PFM. Eventually, TRAP obtains the expected number N of TFs bound to the promoter by summing over the individual probabilities from all L sites in the sequence:

(19)

As input, TRAP requires for each TF a PFM suitable for computing the mismatch energies and a DNA sequence of interest (see [17] for details).

For our study we use a selection of 102 PFMs from the Transfac database version 11.1 [28], which correspond to TFs involved in lymphoid development (see Additional File 1 for TF list). As we are mainly interested in binding sites near the promoter, the analysis was based on the 200 base pairs upstream of the transcription start site (TSS) of the genes. We restrict the analysis to genes with normalized CpG content < 0.5 in their promoter sequence [21], as such genes tends to be expressed in a tissue and stage specific way. In the end, we calculate the affinity (Eq. 19) for all the selected genes and PFMs. This yields the matrix X containing the TF binding data, where x _{i , j} corresponds to the affinity of TF j to the promoter of gene i.

T-cell gene expression and histone modification data

We use the gene expression and histone modification data from Th1, Th2, Th17 and iTreg cells published by [16]. The histone modification data was measure with the Chip-Seq Illumina platform. We used the Cisgenome tool [29] to align sequence data and to detect peaks. As we are only interested in the modifications near the promoter, we consider the region of 8000 bps upstream and 2000 bps downstream of the TSS and kept the tag counts of the highest peak. Finally, we added a pseudo count to avoid zero values and applied a log transform. This yields the matrix X containing the histone modification data, where x _{i , j} corresponds to the number of ChIP-seq tags derived from a particular histone modification j that are being mapped to the promoter of gene i.

The expression data was measured with Affymetrix 430 chips. The raw data has been normalized using the variance stabilization method of [30] and normalized the tissues to have mean expression equal to zero.Microarray probes were mapped to ENSEMBL gene identifiers with the help of the biomart tool [31]. We thereby kept all genes that had their expression measured by multiple probe sets. In the following, we restrict our analysis to those 6154 genes with low CpG content for which both, gene expression as well as histone modification data is available. The final data sets used in this analysis can be found at http://www.cin.ufpe.br/~igcf/MixLin.

Experimental design

We model gene expression from four different T helper cell types (Th1, Th2, Th17 and iTreg) with the use of either transcription factor affinities (TF), histone modifications (HM) or both regulatory signals combined (HM+TF). As parameter of our method, we vary the number of linear models, K, from 1 to 6. In order to select the optimal model for each cell type, we first perform 10 fold cross-validation on each parameter setting and then estimate the Mean Square Errors (MSE) from the validation sets. As the MSE tends to decrease with higher K[32], we use a model selection procedure, the Bayesian Information Criteriun (BIC) [25], to indicate the optimal number of models. The method has been implemented with Pymix [33] and is freely available at http://www.pymix.org.