Gamma-Normal-Gamma Mixture Model for Detecting Differentially Methylated Loci in Three Breast Cancer Cell Lines

Abstract

With state-of-the-art microarray technologies now available for whole genome CpG island (CGI) methylation profiling, there is a need to develop statistical models that are specifically geared toward the analysis of such data. In this article, we propose a Gamma-Normal-Gamma (GNG) mixture model for describing three groups of CGI loci: hypomethylated, undifferentiated, and hypermethylated, from a single methylation microarray. This model was applied to study the methylation signatures of three breast cancer cell lines: MCF7, T47D, and MDAMB361. Biologically interesting and interpretable results are obtained, which highlights the heterogeneity nature of the three cell lines. This underlies the premise for the need of analyzing each of the microarray slides individually as opposed to pooling them together for a single analysis. Our comparisons with the fitted densities from the Normal-Uniform (NU) mixture model in the literature proposed for gene expression analysis show an improved goodness of fit of the GNG model over the NU model. Although the GNG model was proposed in the context of single-slide methylation analysis, it can be readily adapted to analyze multi-slide methylation data as well as other types of microarray data.

Keywords

CpG islands mixture modeling methylation/epigenetic signature microarrays breast cancer cell lines

1. Introduction

Mammalian DNA methylation occurs when a methyl group is added to the cytosine residue of a CpG dinucelotide (Jaenisch and Bird, 2003) and indeed a large portion of genomic DNA is methylated in CpG dinucleotides. Essential exceptions are the CpG islands (CGIs) within gene promoters that are generally unmethylated in normal cells. For more than two decades it has been known that DNA methylation is associated with gene activity and chromatin structure. Recent studies have demonstrated that aberrant CpG methylation is connected to human diseases such as cancer (Laird, 2005; Yang et al. 2006; Wang et al. 2006), Niemann-Pick Disease (Simonaro et al. 2006), lupus (Ballestar et al. 2006; Patole et al. 2005), and Rett syndrome (Kriaucionis and Bird, 2003).

Implementation of state-of-the-art microarray technologies has made possible the mea surement of the methylation signature of multiple genes simultaneously (Piotrowski et al. 2006; Shen et al. 2006; Yan et al. 2002). Methylation signatures, and more generally epi genetic signatures, play an important role in furthering the system-wide understanding of the biological mechanism of disease. Recently, groups have conceptualized, and in a small part realized, the idea of a human epigenome project (Jones and Martienssen, 2005; Rakyan et al. 2004). As the human genome project has advanced discoveries in both basic and translational science, it is expected that the human epigenome project combined with high-throughput methylation assays will facilitate our fundamental understanding of aberrant epigenetic mechanisms.

Two aspects of microarray data make the analysis exceptionally challenging. First, even with recent advances in the microarray technology, there exists experimental variability, which makes detection of genuine biological differences dificult to accomplish as the two may be confounded (Verducci et al. 2006). Second, if p genes are measured in the microar ray experiments, then the microarray data represents a p-dimensional random vector with mutually dependent components. Cost, both monetary and temporal, make it inhibitive to generate more than a small number (with respect to p) of replicates. This “large-p-small-n” scenario prevents the use of classical statistical approaches to data analysis (Verducci et al. 2006).

In some circumstances microarray analysis is further hampered by the lack of any ex perimental or biological replicates. Note that a single-slide analysis (n = 1), when there is no replication at all, is the “large-p-small-n” scenario taken to the extreme. This may be due to lack of resources (such as time, biological samples, arrays, etc.) for a larger experimental study. For example, a pilot study may be interested in the epigenetic differences between tumors in different stages of progression; however, due to limited financial resources only one sample from each stage is assayed. Once the human epigenome project is completed, another conceivable situation is the clinical application of high-throughput arrays for assaying the methylation signature of DNA markers associated with different diseases. Gene expression data from microarrays has recently been used to model interactions of multiple risk factors associated with clinical outcome for use in individual patient care (Potti et al. 2006; Li, 2006; Pittman et al. 2004). The integration of genome wide methylation information from a single slide patient sample will potentially enhance and strengthen the clinical predictive power of such risk models (West et al. 2006).

A number of methods have been suggested for finding genes that are differentially ex pressed between two variants for single-slide data, all of which are proposed for gene ex pression microarray data. The most common approach is to flag those with an intensity ratio above a given cut-off, typically between 1.5 and 3. Such a selection criterion is ad-hoc in nature. More formal statistical treatments are available in the literature. Newton et al. (2001) suggested modeling the Cy5 and Cy3 signals on an array as independent random variables each following a gamma distribution. Such an approach benefits from considering the thousands of measurements on an array simultaneously rather than individually (Efron and Morris, 1973; Carlin and Louis, 1996). A different single-slide applicable approach was recently proposed by Dean and Raftery (2005), which models the normalized log intensity ratios as coming from a two-component, Normal and Uniform (NU), mixture of distributions.

In this paper, we present a methodology for the detection of differentially methylated loci between the two samples co-hybridized onto a methylation microarray. The log ratio of the Cy5 to Cy3 intensities of the probes on the array are modeled as random variables sampled from one of three distributions. These three groups are biologically interpreted as hypomethylated: the assayed DNA from the Cy3 sample is methylated while the Cy5 sample is not; hypermethylated: the Cy5 sample is methylated and the Cy3 sample is not; and undifferentiated: there is no observable difference in methylation status between the two samples. This model was applied to three breast cancer cell lines, MCF7, T47D, and MDAMB361. These analyses yielded biologically meaningful results, which are consistent with the current understanding of the role of aberrant methylation in breast cancer development. A boot strap resampling technique was used to evaluate the reliability of the probe classifications from the fitted model. We further assessed the goodness of fit of this model by comparing it to that of the NU model (Dean and Raftery, 2005) based on the Kullback-Leibler distance measure.

2. Methods

In this section we first introduce a mixture model for identifying differentially methylated loci (probes) in the sample. An algorithm is described for estimating the parameters of the model. Then, we outline the usage of the model for classifying the probes into the three groups: hypomethylated, undifferentiated, and hypermethylated.

2.1. Gamma-Normal-Gamma mixture model

Consider the Red (R; Cy5) and Green (G; Cy3) intensity measurements. Our interest is to model Y, the normalized log ratio, log₂(R/G), by the Gamma-Normal-Gamma (GNG) model. Our assumption here is that certain steps have already been taken to normalize the data first. Briefly, loess normalization transform as described in (Yang et al. 2002) is applied to the data, however, the loess fit is based on a rank-invariant subset as in (Speed, 2003).

We begin by modeling probes as three different groups: hypomethylated, undifferentiated, and hypermethylated. Since the undifferentiated probes have a true log ratio of zero, we model the observed log ratio for these probes by a Normal distribution. For the hyper-methylated probes, the log ratio will be positive and hence a Gamma distribution will be a reasonable choice for modeling this group of probes. The choice of the Gamma model is mainly due to the perception that the characteristics of the distribution match well with the biological data it models. In fact, the Gamma model has a long history in statistical ecology (e.g. Fisher et al. 1943) due to both its analytical convenience and its possession of a potential deeper biological interpretation, and it has recently been used in microrray analysis (Newton et al. 2001). The log ratios corresponding to the hypomethylated probes are negative, and we choose to use the mirror image of the Gamma distribution used for the hypermethylated probes for this group of probes. However, we choose a different set of parameters for the mirror image Gamma distribution. Thus, overall, a mixture of three distributions, that is, one Normal and two different Gamma distributions will be used to model the normalized log ratio Y. Mathematically, the model is given as follows.

Consider the indicator function I_{y≫0} = 1, for y > 0, and zero otherwise. Our proposed mixture model is then given by

\begin{matrix} f (y; Ψ) = π_{1} G (- y; α_{1}, β_{1}) I_{{- y > 0}} \\ + π_{2} N (y; μ, σ^{2}) \\ + π_{3} G (y; α_{2}, β_{2}) I_{{y > 0}} \end{matrix}

(1)

where π_k > 0, k = 1, 2, 3, are called mixing proportions, and they represent the proportion of each of the three above mentioned groups, and further $\sum_{k = 1}^{3} π_{k} = 1.$ The G(·; α, β) and N(·; μ, σ²) stand for the Gamma and Normal density functions, and

Ψ = (α_{1}, β_{1}, α_{2}, β_{2}, μ, σ^{2}, π_{1}, π_{2}, π_{3})

is the vector of all unknown parameters in the model, but note that the π's are interdependent. In (1), the π₁ and π₃ show the probabilities of each probe being hypomethylated and hypermethylated, respectively. The π₂ is the probability that a probe is not differentially methylated. Although the theoretical mean of the normal component is zero, we opted to allow this value to be estimated jointly with the other parameters to capture any residual bias that has not been removed by normalization.

2.2. Parameter Estimation

Since our modeling approach is fully parametric, we use the most popular parameter estimation method, that is maximum likelihood estimation, for estimating the parameter vector Ψ.

Let y_i, i = 1, 2, …, n, be the normalized log ratios. The log-likelihood function of the parameter vector Ψ is given by

l_{n} (Ψ) = \sum_{i = 1}^{n} \log f (y_{i}; Ψ)

where f(y_i; Ψ) is the mixture density function in (1). The Maximum Likelihood Estimate (MLE) of Ψ is then given by

{\hat{Ψ}}_{n} = \arg \max_{Ψ} l_{n} (Ψ) .

Due to the complexity of the log-likelihood function l_n(Ψ), there is no apparent analytical form for the estimator ${\hat{Ψ}}_{n}$ . Thus a numerical method is needed to find the solution to the above maximization problem. In the context of mixture models, the Expectation-Maximization (EM) algorithm of Dempster, Larid and Rubin (1997) provides a convenient approach to obtain the MLE ${\hat{Ψ}}_{n}$ . The EM algorithm in our problem is outlined as follows.

Let z_ik, k = 1, 2, 3, be indicator variables showing the component membership of each observation y_i in the mixture model (1). Note that z_ik's are unobservable (missing) variables. The EM algorithm works with the complete data log-likelihood function

\begin{matrix} l_{n}^{c} (Ψ) = \sum_{i = 1}^{n} \sum_{k = 1}^{3} z_{i k} \log π_{k} \\ + \sum_{i = 1}^{n} {z_{i 1} \log [G (- y_{i}; α_{1}, β_{1})] I_{{- y_{i} > 0}} \\ + z_{i 2} \log [N (y_{i}; μ, σ^{2})] \\ + z_{i 3} \log [G (y_{i}; α_{2}, β_{2})] I_{{y_{i} > 0}}} . \end{matrix}

The EM algorithm finds ${\hat{Ψ}}_{n}$ iteratively in two steps as follows.

E-Step

Let $Ψ^{(m)}$ be the estimate of Ψ after the mth iteration. The E-step computes the conditional expectation of the function $l_{n}^{c} (Ψ) l_{n}^{c} (Ψ)$ with respect to $P (z | y, Ψ^{(m)})$ , where z = {z_ik, i = 1, …, n, k = 1, 2, 3.}, and y = {y₁, …, y_n}. The conditional expectation is found to be

\begin{matrix} Q (Ψ; Ψ^{(m)}) = \sum_{i = 1}^{n} \sum_{k = 1}^{3} w_{i k}^{(m)} \log π_{k} \\ + \sum_{i = 1}^{n} {w_{i 1}^{(m)} \log [G (- y_{i}; α_{1}, β_{1})] I_{{- y_{i} > 0}} \\ + w_{i 2}^{(m)} \log [N (y_{i}; μ, σ^{2})] \\ + w_{i 3}^{(m)} \log [G (y_{i}; α_{2}, β_{2})] I_{{y_{i} > 0}}} . \end{matrix}

Where

\begin{array}{l} w_{i 1}^{(m)} = \frac{π_{1}^{(m)} G (- y_{i}; α_{1}^{(m)}, β_{1}^{(m)}) I_{{- y_{i} > 0}}}{f (y_{i}; Ψ^{(m)})}, \\ w_{i 2}^{(m)} = \frac{π_{2}^{(m)} N (y_{i}; μ^{(m)}, σ^{2 (m)})}{f (y_{i}; Ψ^{(m)})}, \\ w_{i 3}^{(m)} = \frac{π_{3}^{(m)} G (y_{i}; α_{2}^{(m)}, β_{2}^{(m)}) I_{{y_{i} > 0}}}{f (y_{i}; Ψ^{(m)})}, i = 1, …, n, \end{array}

are the conditional expectations of the z_ik's.

M-Step

The M-step on the (m+1)th iteration maximizes the function $Q (Ψ; Ψ^{(m)})$ with respect to Ψ. By maximizing this function with respect to π_k's, μ and σ², we have

\begin{array}{l} π_{k}^{(m + 1)} = \frac{1}{n} \sum_{i = 1}^{n} w_{i k}^{(m)}, k = 1, 2, 3, \\ μ^{(m + 1)} = \frac{\sum_{i = 1}^{n} w_{i 2}^{(m)} y_{i}}{\sum_{i = 1}^{n} w_{i 2}^{(m)}}, \\ σ^{2 (m + 1)} = \frac{\sum_{i = 1}^{n} w_{i 2}^{(m)} {(y_{i} - \bar{y})}^{2}}{\sum_{i = 1}^{n} w_{i 2}^{(m)}} . \end{array}

Maximization with respect to (α_k, β_k), k = 1, 2, needs to be done by using a numerical method. In the current study we used the function optim() in R for this task. This function implements the method of Nelder and Mead (1965) to maximize the objective function.

Starting from an initial value Ψ⁽⁰⁾, the EM algorithm iterates between the E and M-steps until some convergence criterion is satisfied. For example, for a pre-specified value ∊ > 0, the algorithm will stop if

∥ Ψ^{(m + 1)} - Ψ^{(m)} ∥ < ϵ .

This was the criterion used in the analyses carried out in this paper, with ∊ = 10^–12 and the use of the L² norm.

Due to the sensitivity of the results to the choice of starting point, it is recommended that multiple sets of random and non-random initial parameter values be used to increase the chance of finding the global maximum, especially if the likelihood surface is multimodal. For our applications, we used 20 sets of starting values; the estimates that led to the maximum likelihood were taken as our MLE.

2.3. Classification

Finite mixture models provide a model-based approach for classification (McLachlan and Peel, 2000). Consider a finite mixture model with K components which in fact corresponds to K classes. In the model-based classification one first calculates the posterior probability that an observation belongs to any of the classes k = 1, 2, …, K. Then the observation is classified into the class with highest posterior probability. Accordingly, the classification based on the GNG mixture model can be performed as follows.

Consider the log ratios Y_i, i = 1, 2, …, n, that are assumed to follow the GNG mixture model in (1). The posterior probabilities that Y_i belongs to any of the classes 1, 2 or 3 are given by

\begin{array}{l} p_{i 1} = P (Y_{i} \in 1 | Y_{i} = y_{i}) = \frac{π_{1} G (- y_{i}; α_{1}, β_{1}) I_{{- y_{i} > 0}}}{f (y_{i}; Ψ)}, \\ p_{i 2} = P (Y_{i} \in 2 | Y_{i} = y_{i}) = \frac{π_{2} N (y_{i}; μ, σ^{2})}{f (y_{i}; Ψ)}, \\ p_{i 3} = P (Y_{i} \in 3 | Y_{i} = y_{i}) = \frac{π_{3} G (y_{i}; α_{2}, β_{2}) I_{{y_{i} > 0}}}{f (y_{i}; Ψ)} . \end{array}

Using the MLE ${\hat{Ψ}}_{n}$ , let ${\hat{p}}_{i k}$ k = 1, 2, 3, be the MLE's of the posterior probabilities. A probe with log ratio Y_i is classified to class k if

{\hat{p}}_{i k, \max} = \max {{\hat{p}}_{i 1,} {\hat{p}}_{i 2,} {\hat{p}}_{i 3,}} .

To increase one's confidence on the probes classified, especially those labeled as differentially methylated, and to guard against false positives, we decided to classify a probe with log ratio Y_i to class k only if

\frac{{\hat{p}}_{i k, \max}}{{\hat{p}}_{i j}} > C

for some constant value C ≥ 1, for all j ≠ k; otherwise that particular probe is left as unclassified. We suggest to set the value of C to be inversely related to the confidence of the technology generating the data. That is, a small value for C should correspond to a high confidence in the sensitivity and reproducibility of the microarray measurement. In the applications to the three cell lines, a value of 2 was chosen for C to reflect our degree of confidence in the microarray experiments, which was also partly motivated by ease of interpretation: a two-fold difference in probabilities was required to determine classification. In the situation where a priori information about the quality of the microarray measurement is unavailable, a bootstrap technique as described below may be used as a means to select an appropriate C that will lead to high consistency rates.

3. Results

3.1. Breast cancer data

Three cell lines

Breast cancer cell lines have been used as research models for a long time. A large number of them were established in the late 1970s. Two of the cell lines we chose to analyze, MCF7 and T47D, are among the most often studied breast cancer cell lines. They share distinctive features such as isolated originally from pleural effusions, little to no invasiveness, hormone-receptor (estrogen-and progesterone receptor) positive, Her-2/neu negative, expression of WNT7B oncogene, and high expression of genes associated with the luminal epithelial-like phenotype. Although there are many similarities between these two cell lines, there also exist differences. First, the age of the patients from which these two cell lines derived from were 69 and 54 for MCF7 and T47D, respectively. As DNA methylation increases with age, we would expect more methylation events in samples from older patients. Second, though both of these cell lines express estrogen receptor, it was noted by Karey and Sirbasku (1988) that they responded differently to the challenge of 17 beta-estradiol. Third, due to known (caspase 3 is expression in T47D but not in MCF7) and unknown gene expression profiles associated with the cell lines, they were shown to respond differently to gene perturbation (Yamashita et al. 2003) and to drug- and chemical treatments (Green, 2003; Bhat, 2001). Although the MDA-MB-361 cell line shares only one common feature with the other two cell lines (positive hormone receptor status) and has many distinct features (isolated originally from a brain metastasis of a 40 years old patient, moderate in vitro invasiveness, Her-2/neu overexpression, and expression of WNT7H + oncogene), it still tends to be closer to the luminal epithelial-like MCF7 and T47D than to the mesenchymal-like breast cancer cell lines. In all, we expect these three cell lines to share some commonly methylated loci and harbor aberrant methylation in loci that were uniquely theirs.

Generation of methylation profiles

The array employed to generate the methylation profiles for the three cell lines was custom designed to assay the promoter regions located within a CpG island (regions of the DNA with higher than expected CpG dinucleotides). The oligonucleotide probes printed on the array were selected from a larger library of probes comprising Agilent's location analysis array (Boyer et al. 2005). Oligonucleotide probes were considered if their target was located within the predicted promoter region of a gene. This subset was further reduced so that the resulting array contained 40659 unique probes with targets located within the predicted promoter region of 10834 unique genes. For each cell line, DNA was isolated from the cells, digested by the restriction enzyme BfaI in order to reduce genomic complexity, and ligand linked for subsequent PCR. Methylation sensitive enzymes, HpaII and HinP Il, were used to interrogate the methylation status of the DNA fragments. The remaining intact DNA was then PCR amplified, dye coupled with either Cy5 (cancer samples) or Cy3 (normal samples), and hybridized onto the CpG promoter methylation array as described above. Thus, we have a total of three slides, each captures the methylation signature of a particular cell line. In all three slides, the normal samples come from the immortalized normal breast epithelial cell line MCF10A.

Data preprocessing

By utilizing the ink jet technology, Agilent has eliminated spatial effect often found correlated with print-tip regions on other oligo arrays. Upon inspection of the heat map for the normalized log-ratio of the probes with respect to their position on the array, see Figure 1, no obvious spatial effects were apparent. Therefore, no spatial correction was made. However, Figure 2, A–C show a significant, yet typical, dye effect, and thus the data were normalized as described above. Since the spread of the resulting mean-corrected log-ratios was not significant (0.33 ± 0.03), no methods for adjusting the spread of the data was used (see Figure 2, D–F).

Figure 1.

Heat map for the normalized log-ratio of the Cy5 to Cy3 intensity for the three analyzed data sets. Color scheme represents relative ranking of the log-ratio with green denoting a low ranking, or relatively small ratio, and red denoting a high ranking, or relatively large ratio. The x and y axis of the plots denote the position of the probe on the microarray.

Figure 2.

Scatter plots for the average log intensity versus the log ratio for each of the three data sets. A–C: The un-normalized data is plotted. D–F: The loess normalized data is plotted with the hyper- and hypomethylated probes highlighted in red and green, respectively. The unclassified probes are highlighted in blue. The black bullets in D are some of the genes that we had validated as hypermethylated.

3.2. Model Fitting

We fitted the GNG model to the normalized data form the three microarray experiments. The parameter estimates corresponding to each experiment are given in Table 1. For all three cell lines, the mixing proportions follow approximately a 15-70-15% split, although deviations among one another is obvious. Similar observations apply to the estimates of the other parameters. In particular, we note that the mean estimates for the normal component are all close to 0, as one would expect with the normalized data, but the three estimates deviate from one another. To check the performance of the fitted GNG model empirically, we plotted the density function of the fitted models over the histograms of the normalized log ratios. Furthermore, QQ-plots are also provided. As can be seen from Figure 3, apart from a few points at the end of the QQ plots not falling right on the straight lines, the plots show that the GNG model provided reasonably good fits to each of the microarray data. Using the classification criterion discussed in 2.3, probes were classified as hypermethylated (red dots in Figure 2, D–F), undifferentiated (black dots), hypomethylated (green dots), or left unclassified (blue dots). This resulted in a total of 965, 1051, and 375 genes, in MCF7, T47D, and MDAMB361, respectively, labeled as hypermethylated, the group that we will focus on in discussing the biological relevance in the next sub-section.

Table 1.
Parameter Estimates in the Three Experiments.

Estimates MCF7 T47D MDAMB361

$({\hat{α}}_{1}, {\hat{β}}_{1})$ (0.94, 0.42) (1.07, 0.31) (0.95, 0.37)

$({\hat{α}}_{1}, {\hat{β}}_{1}) ({\hat{α}}_{1}, {\hat{β}}_{1})$ (0.002, 0.24) (–0.01, 0.21) (0.01, 0.20)

$({\hat{μ}}_{1}, {\hat{σ}}_{1})$ (1.11, 0.31) (0.93, 0.36) (0.91, 0.24)

$({\hat{π}}_{1}, {\hat{π}}_{2}, {\hat{π}}_{3})$ (0.16, 0.68, 0.16) (0.13, 0.73, 0.15) (0.16, 0.71, 0.14)

Estimates	MCF7	T47D	MDAMB361
$({\hat{α}}_{1}, {\hat{β}}_{1})$	(0.94, 0.42)	(1.07, 0.31)	(0.95, 0.37)
$({\hat{α}}_{1}, {\hat{β}}_{1}) ({\hat{α}}_{1}, {\hat{β}}_{1})$	(0.002, 0.24)	(–0.01, 0.21)	(0.01, 0.20)
$({\hat{μ}}_{1}, {\hat{σ}}_{1})$	(1.11, 0.31)	(0.93, 0.36)	(0.91, 0.24)
$({\hat{π}}_{1}, {\hat{π}}_{2}, {\hat{π}}_{3})$	(0.16, 0.68, 0.16)	(0.13, 0.73, 0.15)	(0.16, 0.71, 0.14)

Figure 3.

Density plots and QQ plots of fitted model.

Despite the fact that we know there are three groups of probes: hypomethylated, hyper-methylated and undifferentiated, one may also consider fitting the NU model (Dean and Raftery, 2005), which was proposed to identify differentially expressed genes from non-differentially expressed ones. To compare the relative goodness of fits of the GNG and the NU models to each of the dataset, we use the Kullback-Leibler (K-L) distance. In probability and information theory, the Kullback-Leibler distance (or relative entropy) is a natural distance measure between two probability distributions with densities functions, say, f(x) and g(x), and is denned as

K L (f, g) = \int_{- \infty}^{+ \infty} f (y) \log \frac{f (y)}{g (y)} d y .

(2)

The K-L distance is always non-negative, and KL(f, g) = 0 if and only if f = g. In applications, typically f represents the empirical distribution of the observed data, and g represents a proposed approximation of f. In the current application, the proposed GNG or the NU mixture plays the role of g, after the fitted density is discretized, and thus (2) becomes a summation. Table 2 shows the K-L distance based on each of the GNG and NU mixture models as well as a measure of reduction of distance, for each the three experiments under consideration. As we can see from the table, the KL distances are much larger under the NU model than under the GNG model.

Table 2.

Kullback-Leibler Distance Between the Fitted Models (GNG or NU) and the Observed Data.

K-L Distance	MCF7	T47D	MDAMB361
KL(GNG, Obs.)	0.00055	0.00047	0.00086
KL(UN, Obs.)	0.022	0.017	0.026
(KL(UN, Obs.) – KL(GNG, Obs.))/KL(UN, Obs.)	0.97	0.97	0.97

We further investigated the reliability of the GNG model for classification using a bootstrap resampling technique. For each of the three experiments under study, we generated a bootstrap sample in each experiment and re-classified the probes flagged in the original data based on the model fitted with the bootstrap sample. We then calculated the proportion of the hypermethylated loci that retained the same classification as the originals, and referred to it as the rate of consistency. We generated B = 200 bootstrap samples for each of the three experiments, which resulted in a mean consistency rate of 94.6, 98.4, and 99.0%, respectively.

3.3. Biological relevance of hypermethylated genes

For the given data set, the GNG model identifies genes for which their promoter regions are hyper-as well as hypomethylated in the breast cancer cell lines. However, as our research focus is on genes down-regulated by DNA methylation, we will confine our discussion concerning the biological relevance of the target genes only to the subset of hypermethylated genes.

Cancer cell lines are known to have more methylation events than primary tumors due to the effect of long-term culturing of DNA methylation. This is evident in the number of methylated genes identified by the GNG model in the breast cancer cell lines. T47D has many methylated genes (1,051), followed closely by MCF7 (965), and MDA-MB-361 has significantly less methylated targets (375). This is equivalent to about 4–10% methylation as our microarray platform contains close to 11,000 unique genes. In primary breast tumors, we expect to detect around 1–2% methylation events. Of the genes flagged to be hypermethylated, 86 of them are common in all three cell lines. The Venn diagram in Figure 4 also shows that MCF7 and T47D share 212 additional hypermethylated genes, which more than doubles that shared by MDAMB361 and either of these two cell lines. This finding is consistent with our discussion of the characteristics of the cell lines earlier in Section 3.1.

Figure 4.

Venn diagram denoting the intersection between the sets of hypermethylated genes for each of the three experiments. The radius of each circle is relative to the number of genes in each set.

The Gene ontology tree machine (GOTM) Zhang et al. (2004) was used to conduct a statistical analysis of the GO terms enriched by the 86 commonly methylated genes. Using the entire genome as the reference population, the category “Sequence-specific DNA binding” is enriched in this set (p-value = 0.0019). Of the seven genes (HOXB13, FALZ, FOXO3A, GATA2, HOXB9, HOXD11, and TERT) present in this category, a couple of them deserve more detailed discussion. HOXB13 is a gene we have studied extensively before (Rodriguez et al. 2006). This gene encodes a transcription factor that belongs to the homeobox gene family. It was evaluated through quantitative methylation-specific PCR (qMSP) that higher level of methylation occurred in hormone-receptor positive primary breast tumors, corroborating our current microarray finding. Another sequence-specific DNA binding gene identified is FALZ. The protein coded by FALZ is highly similar to the largest subunit of the Drosophila NURF (nucleosome remodeling factor) complex which interacts with sequence-specific transcription factors to facilitate gene transcription. There is evidence that human NURF complex can act as a negative regulator of the JAK/STAT (cytokine signaling) pathway thereby decreasing growth rate and anchorage-independent growth of cancer cells. As such, hypermethylation of FALZ can confer growth advantage to these cell lines in many diverse ways. Another member of the sequence-specific DNA binding gene list, FOXO3A, belongs to the forkhead family of transcription factors. This gene likely functions as a trigger for apoptosis through expression of genes necessary for cell death. Therefore, for cell lines such as T47D that lacks caspase 3 (Radisavljevic, 2003), the FOXO3a kinase pathway provides a much need means to trigger cellular apoptosis.

In addition to the GO analysis, we decided to verify the methylation status of some of the other genes. Genes were first checked against the literature to establish that their suppression might play a role in cancer development and progression. Then a subset of genes with probes on BfaI fragment(s) that overlapped promoter CpG island were selected for validation by methylation-specific PCR in MCF7 cell lines (MSP) (Herman et al. 1996). Genes shown by MSP to be methylated are identified in Figure 3D as brown bullets. It is evident from the plot that the methylation levels of these gene targets varied from moderately methylated to highly methylated.

Yet another approach to evaluate hypermethylated genes identified by the GNG model is to seek out genes evaluated in other breast cancer studies. One such gene is LATS1. The GNG model identified this gene to be methylated in all three cell lines studied. Experimental evidence suggests that this gene is a tumor suppressor gene capable of modulating cell survival and negatively regulating cell proliferation by modulating CDC2/Cyclin A activity (Xia et al. 2002). Another gene found to be hypermethylated in MDA-MB-361 and T47D is ETV5, also known as ERM. ETV5 is a member of the superfamily of Ets-related transcription factors. Members of this family bind to related DNA sequence elements in the promoters of target genes to regulate gene transcription. Candidate target genes include ERBB2 (encodes Her2/neu) and genes encode matrix proteases (Bieche et al. 2004). Therefore, dysregulation of ETV5 (by DNA methylation) would lead to increase invasive potential of tumor cells. The observation of promoter hypermethylation in this gene in MDA-MB-361 is reasonable as this cell line overexpresses Her2/neu. Such observation in T47D but not in MCF7 is also reasonable as tumors occurring at a younger age (T47D originated from a 54 years old patient) tends to be more aggressive than tumors occurring at an older age (MCF7 originated from a 69 years old patient). Contrarily to ETV5, WT1 was found to be hypermethylated in MCF7 and T47D but not in MDA-MB-361. The WT1 gene was originally identified as a tumor suppressor gene responsible for Wilms’ tumor, a kidney neoplasm of childhood. Recent findings point to an oncogenic effect of WT1 in tumor types such as leukemia and breast cancer. In the breast cancer study, (Miyoshi et al. 2002) observed that a high expression of WT1 conferred a poor prognosis for breast cancer patients independent of other conventional prognostic factors. Therefore, it is reasonable to observe methylation of WT1 promoter in the less aggressive MCF7 and T47D but not in the more aggressive MDA-MB-361 cell line.

4. Discussion

In this paper, we proposed a single slide data analysis model (GNG mixture model) for identifying hypermethylated CpG islands loci based on methylation microarray data, and applied it to analyze data generated from three estrogen receptor-positive cell lines. From a statistical perspective, the model fits the data well, and outperforms the NU mixture model of Dean and Raftery (2005) that was proposed for gene expression analysis. To get a rough idea on how the GNG model would perform for gene expression data and to further compare with the NU model for such data, the GNG model was used to fit the gene expression data analyzed in Dean and Raftery (2005). In this preliminary study, the genes that were differentially expressed according to the NU mixture model were similarly classified by the GNG model (data not shown), although differences exist in some cases. These initial findings warrant further inquiry into the general applicability of the GNG model for classifying data obtained from various types of two-color arrays, but it is not within the scope of the current paper.

Through a bootstrap resampling technique, we evaluated the reliability of the loci that are flagged as hypermethylated. For all three cell lines, the call consistency was at least close to 95% based on B = 200 bootstrap samples, which provided us with high confidence in the loci labeled as hypermethylated. As discussed earlier regarding making inferences on classifying a probe as hypermethylated, we suggested setting the strength-of-evidence constant C to be inversely proportional to the ‘reliability’ of the microarray. Such a priori information may be learned through a standard known test set. However, if such information is not readily available, then the bootstrap procedure might be applied to various constants C's. The one that would lead to a high consistency rate could then be chosen for making inferences.

Four lines of evidence were sought to evaluate the biological relevance of the loci identified as hypermethylated among the three cell lines: the number of singularly and commonly (both two ways and three ways) methylated loci; the gene ontology analysis of the commonly methylated loci in all three cell lines; the experimental verification of a subset of the flagged loci using methylation-specific PCR; and the literature supports of key methylated loci identified. Since the three cell lines are all estrogen receptor-positive with key similarities, but there are also clear differences in their phenotypic characteristics, the results match well with such biological knowledge. The GNG model revealed that MCF7 and T47D cell lines had similar amount of hypermethylated genes whereas the cell line MDA-MB-361 had significantly less methylated targets. At this global level, the GNG model portrays a relationship between DNA methylation and aging. A more aggressive tumor increases in severity quickly (such as MDA-MB-361) thereby progresses with little accumulation of methylation marks. Tumors from older patients (such as MCF7) develop less aggressively thereby accumulate more methylation events along the way. With the advent of microarray technology (current microarray platform interrogates close to 11,000 gene promoters), the GNG model uncovers many methylation events previously not described in the literature. We chose to validate a small set of genes identified as hypermethylated in MCF7. These genes had varying degree of methylation in comparison to immortalized normal breast epithelial cell line MCF10A. This experimental outcome in part substantiates the usefulness of the GNG model. A survey of the literature uncovers gene targets methylated concurrently or in a subset of the three breast cancer cell lines studied as described by other researchers. Methylation profiles of some of these genes are highlighted (not exhaustively) to provide additional support for these gene targets (black bullets in Figure 2D).

It would also be of interest to evaluate the biological relevance of the methylated genes through its cross-correlation with gene expression data. To this end, we briefly examined a few of the genes that are commonly assayed in both the Affymetrix gene expression array based on which the gene expression data on the same cell lines were generated and the 44 K methylation array. It appears that some of the genes that are hypermethylated are indeed down-regulated in their gene expression, but given the small number of genes examined, it is difficult to draw any general conclusion at this point. However, it would be of great interest to find the genes that are both hypermethylated and down-regulated, as they might be potential biomarkers for cancer diagnosis, prognosis, or even targets of demethylation drugs.

Finally, we note that although the GNG model was proposed for analyzing single slide methylation data due to the need for such a procedure, the formulation can be easily ex tended to the analysis of multi-slide data. One possibility would simply be using the average normalized log ratios, following Dean and Raftery (2005). Another would use normalized (within and between slides) log ratios from each slide directly in the joint likelihood formulation to increase the sample size and hence accuracy of the parameter estimates. Be it as a single-, or a multi-slide analysis tool, the GNG model can be used for data from other microarray platforms, such as gene expression data.

Footnotes

Acknowledgments

This work is supported by the National Cancer Institute grant U54CA113001. DP is supported in part by the NCI grant T32-CA106196-03. This material is also based upon work partially supported by The National Science Foundation grants DMS-0306800 and DMS-0112050.

References

Ballestar

, Esteller

, and Richardson

B.C.

2006. The Epigenetic Face of Systemic Lupus Erythematosus. J. Immunol., 176(12): 7143–7147.

Bieche

, Tozlu

, Girault

, Onody

, Driouch

, Vidaud

, and Lidereau

2004. Expression of PEA3/E1AF/ETV4, an Ets-related transcription factor, in breast tumors: positive links to MMP2, NRG1 and CGB expression. Carcinogenesis, 25(3): 405–411.

Boyer

L.A.

, Lee

T.I.

, Cole

M.F.

, Johnstone

S.E.

, Levine

S.S.

, Zucker

J.P.

, Guenther

M.G.

, Kumar

R.M.

, Murray

H.L.

, Gifford

R.G.J.D.K.

, Melton

D.A.

, Jaenisch

, and Young

R.A.

2005. Core Transcriptional Regulatory Circuitry in Human Embryonic Stem Cells. Cell, 122(6): 947–956.

Carlin

, and Louis

1996. Bayes and Empirical Bayes Methods for Data Analysis. New York: Chapman and Hall.

Dean

, and Raftery

2005. Normal uniform mixture differential gene expression detection for cDNA microarrays. BMC Bioinformatics, 6(1): 173.

Efron

, and Morris

1973. Combining possibly related estimation problems (with discussion). Journal of the Royal Statistical Society, Series B, 35: 379421.

Herman

J.G.

, Graff

J.R.

, Myohanen

, Nelkin

B.D.

, and Baylin

S.B.

1996. Methylation-specific PCR: A novel PCR assay for methylation status of CpGislands. PNAS, 93(18): 9821–9826.

Jaenisch

, and Bird

2003. Epigenetic regulation of gene expression: how the genome integrates intrinsic and environmental signals. Nat. Genet., 33: 245–254.

Jones

P.A.

, and Martienssen

2005. A Blueprint for a Human Epigenome Project: The AACR Human Epigenome Workshop. Cancer Res., 65(24): 11241–11246.

10.

Karey

, and Sirbasku

1988. Differential responsiveness of human breast cancer cell lines MCF-7 and T47D to growth factors and 17 beta-estradiol. Cancer Res., 48(14): 4083–4092.

11.

Kriaucionis

, and Bird

2003. DNA methylation and Rett syndrome. Hum. Mol. Genet., 12(Review Issue 2): R221–227.

12.

Laird

P.W.

2005. Cancer epigenetics. Hum. Mol. Genet., 14(Review issue 1): R65–76.

13.

2006. Survival prediction of diffuse large-B-cell lymphoma based on both clinical and gene expression information. Bioinformatics, 22(4): 466–471.

14.

McLachlan

G.J.

, and Peel

2000. Finite Mixture Models. New York, NY: Wiley.

15.

Miyoshi

, Ando

, Egawa

, Taguchi

, Tamaki

, Sugiyama

, and Noguchi

2002. High Expression of Wilms’ Tumor Suppressor Gene Predicts Poor Prognosis in Breast Cancer Patients. Clin. Cancer Res., 8(5): 1167–1171.

16.

Nelder

J.A.

, and Mead

1965. A Simplex Method for function minimization. Computational Journal, 7: 308–313.

17.

Newton

, Kendziorski

, Richmond

, Blattner

, and Tsui

2001. On Differential Variability of Expression Ratios: Improving Statistical Inference about Gene Expression Changes from Microarray Data. Journal of Computational Biology, 8(1): 37–52.

18.

Patole

P.S.

, Zecher

, Pawar

R.D.

, Grone

H.J.

, Schlondorff

, and Anders

H.J.

2005. G-Rich DNA Suppresses Systemic Lupus. J. Am. Soc. Nephrol., 16(11): 3273–3280.

19.

Piotrowski

, Benetkiewicz

, Menzel

, de Stahl

T.D.

, Mantripragada

, Grigelionis

, Buckley

P.G.

, Jankowski

, Hoffman

, A D.B., Srutek

, Laskowski

, Zegarski

, and Du manski

J.P.

2006. Microarray-based survey of CpG islands identifies concurrent hyper- and hypomethylation patterns in tissues derived from patients with breast cancer. Genes, Chromosomes and Cancer, (7): 656–667.

20.

Pittman

, Huang

, Dressman

, Horng

C.F.

, Cheng

S.H.

, Tsou

M.H.

, Chen

C.M.

, Bild

, Iversen

E.S.

, Huang

A.T.

, Nevins

J.R.

, and West

2004. Integrated modeling of clinical and gene expression information for personalized prediction of disease outcomes. PNAS, 101(22): 8431–8436.

21.

Potti

, Mukherjee

, Petersen

, Dressman

H.K.

, Bild

, Koontz

, Kratzke

, Watson

M.A.

, Kelley

, Ginsburg

G.S.

, West

, Harpole

, David

, and Nevins

J.R.

2006. A Genomic Strategy to Refine Prognosis in Early-Stage Non-Small-Cell Lung Cancer. N. Engl. J. Med., 355(6): 570–580.

22.

Radisavljevic

2003. Nitric oxide suppression triggers apoptosis through the FKHRL1 (FOXO3A)/ROCK kinase pathway in human breast carcinoma cells. Cancer, 97(5): 1358–1363.

23.

Rakyan

V.K.

, Hildmann

, Novik

K.L.

, Lewin

, Tost

, Cox

A.V.

, Andrews

T.D.

, Howe

K.L.

, Otto

, Olek

, Fischer

, Gut

I.G.

, Berlin

, and Beck

2004. DNA Methylation Profiling of the Human Major Histocompatibility Complex: A Pilot Study for the Human Epigenome Project. J. PLoS Biology, 2(12): e405.

24.

Shen

, Chow

, Wang

, and Fan

2006. Abnormal CpG island methylation occurs during in vitro differentiation of human embryonic stem cells. Hum. Mol. Genet., 15(17): 2623–2635.

25.

Simonaro

C.M.

, Park

J.H.

, Eliyahu

, Shtraizent

, McGovern

M.M.

, and Schuchman

E.H.

2006. Im printing at the SMPD1 locus: implications for acid sphingomyelinase-deficient Niemann-Pick disease. Am. J. Hum. Genet., 78(5): 865–870.

26.

Speed

2003. Statistical Analysis of Gene Expression Microarray Data. Chapman & Hall/CRC.

27.

Verducci

J.S.

, Melfi

V.F.

, Lin

, Wang

, Roy

, and Sen

C.K.

2006. Microarray analysis of gene expression: considerations in data mining and statistical treatment. Physiol. Genomics, 25(3): 355–363.

28.

Wang

L.G.

, Beklemisheva

, Liu

X.M.

, Ferrari

A.C.

, Feng

, and Chiao

J.W.

2006. ual action on promoter demethylation and chromatin by an isothiocyanate restored GSTP1 silenced in prostate cancer. Molecular Carcinogenesis.

29.

West

, Ginsburg

G.S.

, Huang

A.T.

, and Nevins

J.R.

2006. Embracing the complexity of genomic data for personalized medicine. Genome Res., 16(5): 559–566.

30.

Xia

, Qi

, Li

, Pei

, Barton

, Blackstad

, Xu

, and Tao

2002. LATS1 tumor suppressor regulates G2/M transition and apoptosis. Oncogene, 21(8): 1233–1241.

31.

Yan

P.S.

, Chen

C.M.

, Shi

, Rahmatpanah

, Wei

S.H.

, and Huang

T.H.M.

2002. Applications of CpG Island Microarrays for High-Throughput Analysis of DNA Methylation. J. Nutr., 132(8): 2430S–2434.

32.

Yang

, Liu

, Wang

, Tsang

, and Ngan

2006. Differential DNA methylation profiles in gynecological cancers and correlation with clinico-pathological data. BMC Cancer, 6(12).

33.

Yang

Y.H.

, Dudoit

, Luu

, Lin

D.M.

, Peng

, Ngai

, and Speed

T.P.

2002. Normalization for cDNA microarray data: a robust composite method addressing single and multiple slide systematic variation. Nucl. Acids Res., 30(4): e15.

34.

Zhang

, Schmoyer

, Kirov

, and Snoddy

2004. GOTree Machine (GOTM): a web-based platform for interpreting sets of interesting genes using Gene Ontology hierarchies. BMC Bioinformatics, 5(1): 16.