Prediction of CpG-island function: CpG clustering vs. sliding-window methods
© Hackenberg et al. 2010
Received: 26 January 2010
Accepted: 26 May 2010
Published: 26 May 2010
Skip to main content
© Hackenberg et al. 2010
Received: 26 January 2010
Accepted: 26 May 2010
Published: 26 May 2010
Unmethylated stretches of CpG dinucleotides (CpG islands) are an outstanding property of mammal genomes. Conventionally, these regions are detected by sliding window approaches using %G + C, CpG observed/expected ratio and length thresholds as main parameters. Recently, clustering methods directly detect clusters of CpG dinucleotides as a statistical property of the genome sequence.
We compare sliding-window to clustering (i.e. CpGcluster) predictions by applying new ways to detect putative functionality of CpG islands. Analyzing the co-localization with several genomic regions as a function of window size vs. statistical significance (p-value), CpGcluster shows a higher overlap with promoter regions and highly conserved elements, at the same time showing less overlap with Alu retrotransposons. The major difference in the prediction was found for short islands (CpG islets), often exclusively predicted by CpGcluster. Many of these islets seem to be functional, as they are unmethylated, highly conserved and/or located within the promoter region. Finally, we show that window-based islands can spuriously overlap several, differentially regulated promoters as well as different methylation domains, which might indicate a wrong merge of several CpG islands into a single, very long island. The shorter CpGcluster islands seem to be much more specific when concerning the overlap with alternative transcription start sites or the detection of homogenous methylation domains.
The main difference between sliding-window approaches and clustering methods is the length of the predicted islands. Short islands, often differentially methylated, are almost exclusively predicted by CpGcluster. This suggests that CpGcluster may be the algorithm of choice to explore the function of these short, but putatively functional CpG islands.
The methylation of CpG dinucleotides is an important epigenetic modification of DNA, required in mammals for embryonic development, genomic imprinting and X-chromosome inactivation [1–3]. Around 80% of all CpG dinucleotides are methylated in mammal genomes. The exceptions are short stretches of CpG dinucleotides (CpG islands or CGIs), which are predominantly hypomethylated in healthy tissues [4, 5]. CGIs are thought to be predominantly located in the promoter region of genes; around 70% of all genes have a CGI overlapping its promoter region. Moreover, virtually all housekeeping genes are associated to CGIs, while only half of the tissue specific genes show such association . Given its location in the promoters, CGIs may play important roles in the regulation of gene expression. An example is the aberrant methylation of CGIs observed in many cancer types [7–11]. Moreover, evidence exist that the differential or tissue specific methylation of CpG islands may be involved in the regulation of tissue specific genes .
Accurate prediction tools are therefore needed and a considerable effort has been carried out over the last decade to detect CGIs in mammal genomes. Many different algorithms have been proposed, most of them based on the criteria of Gardiner-Frommer . These authors proposed in 1987 thresholds for the detection of CGIs: GC-content (50%), CpG observed/expected (O/E) ratio (0.6) and length (200 bp). Many of the published methods simply readjust these thresholds. However, it has been shown that filtering criteria-based definitions of CpG islands are mathematically incomplete and non-operational, as the sliding window methods frequently fail to identify a large percentage of subsequences that meet the filtering criteria .
Recently, methods based on the clustering of CpGs along the genome sequence detect CGIs as a statistical property, thereby not relying on thresholds of GC-content, O/E ratio and length. The first algorithm published in this category was the CpGcluster method , which detects the CGIs by means of the distances between CpGs, then assigning a statistical significance to each cluster of CpG dinucleotides. Subsequently, CpGcluster was followed by other methods detecting CGIs by means of the CpG densities [15–18]. In the same way, many other features could also contribute to determine the boundary of individual CpG-islands, such as transcription factors and nucleosome location. The nucleosome code could be an important ingredient of future CGI models, although sequence features will probably remain as the principal component (see, for example, ). Epigenetic information may be also of help in detecting CGIs by making use of contextual information .
Given the conceptual differences between sliding window algorithms (SWA) using a high parameter space and those detecting CGIs as a statistical property of the CpG clustering in DNA sequences, disagreement exist on the way CGIs should be predicted. Recently, a comparison between islands detected by the window-based Takai-Jones (TJ) program  and those detected by CpGcluster was published . The comparison evaluated mainly the co-localization of CGIs and known promoters and concludes an overall advantage for the TJ approach over CpGcluster.
We present here new ways to detect putative function of CGIs, emphasizing the basic difference between CpGcluster and SWA predictions: the statistical significance introduced by CpGcluster instead of the conventional length threshold. We show that the statistical significance assigned to each CpGcluster island is a key criterion to control the overlap with promoter regions, evolutionarily conserved elements and spurious Alu elements. Finally, we show that many short (<200 bp) islands (CpG islets) may be also functional, given its overlap with either promoter or evolutionary conserved regions and the absence of methylation in at least one tissue. As many of these islets are exclusively predicted by CpGcluster, this may be the algorithm of choice for experimental essays aimed to verify the function of these short islands.
Therefore, the first part of this work is basically aimed to clarify: 1) the differences between the length threshold used by SWA and the statistically significance used by CpGcluster; and 2) the consequences that the differences in the number of predicted islands and the mean length might have on the prediction quality.
Prediction quality has been assigned conventionally by the percentage of overlap with promoter regions and spurious Alu elements. In the original publication of CpGcluster  we added the overlap with evolutionarily conserved elements or PhastCons  as an indicator of putative functionality. Here, we add several new types of analysis to assess the prediction quality, namely the capability to distinguish between different methylation domains or different alternative Transcription Start Sites (TSSs).
Co-localization of CpG islands and the promoter region.
Number of predicted islands
Genome coverage (%)
Promoter overlap (R13)
Number of islands
Correspondence between the number of predicted islands, log (p-value) and window length.
No. of predicted islands
The results in Figures 3, 4, 5 are straightforward in comparing the relative strengths of the two main parameters involved in CGI quality (length and p-value). The increased stringency in the conventional parameters used by the TJ program excluded contaminating Alu elements, but it also reduced the number of gene promoter associated islands, suggesting that bona fide CGIs were also being discarded . However, raising the statistical significance (i.e. decreasing the p-value) of CpGcluster leads to an exponential increase in the overlap with promoters or PhastCons, simultaneously decreasing the overlap with Alu elements. CpGcluster algorithm is, therefore, a more rational and powerful way to increase CGI prediction quality. An additional advantage is that CpGcluster p-value would be particularly useful in comparative genomics of CGIs, making possible the comparison of CGIs with the same statistical significance, but obtained from different species, despite variations in G + C content or CpG density.
The most important criterion to assess putative functionality of a CpG island is the absence of methylation. Therefore, the comparison to experimentally verified, unmethylated regions is another important analysis type to establish prediction quality.
Recently, the methylation status of 697 hypermethylated and 6,987 hypomethylated promoter regions in WI38 primary lung fibroblast  have been used to compare the prediction quality of TJ and CpGcluster algorithms . In this study, the prediction quality was measured in the following way: i) true positives (TP): hypomethylated promoters containing a predicted island, ii) false positives (FP): hypermethylated promoters containing a predicted island, iii) true negatives (TN): hypermethylated promoters not containing a predicted island, and iv) false negatives (FN): hypomethylated promoters not containing a predicted island.
However, in our opinion, there is an important pitfall in such an approach. It is known that the methylation state of a given region can change among different tissues; therefore, assigning a "false positive" label to a predicted island which has been shown to be methylated in a single tissue may be misleading, as the same prediction could be perfectly "true positive" if measured in a different tissue.
Fortunately, Weber et al.  also determined the methylation states in sperm. Analyzing fibroblast and sperm data together, we observed that 11,260 regions are unmethylated in both tissues but 1,550 are unmethylated in one tissue but methylated in the other one. This means that around 12% of the regions are differentially methylated; therefore, a substantial number of FPs were actually TPs. Given these data, in our opinion, without the knowledge of the methylation state in a vast number of different tissues, the number of "false positive predictions" cannot be assessed in this way.
Prediction of unmethylated regions (Bird's islands, N = 17,383).
Number of predicted islands
Number of islands overlapping a Bird's island
Number of Bird's islands 'touched' by the prediction
Prediction of unmethylated regions (Weber's regions, N = 13,277).
Number of predicted islands
Number of islands overlapping a Weber's region
Number of Weber's regions 'touched' by the prediction
Overlap of different CGIs with 3,465 domains bound by the polycomb repressive complex 2 (PRC2).
Number of predicted islands
Number of islands overlapping PRC2 domains
Number of PRC2 domains 'touched' by the prediction
Given these facts, it might be that either conventional SWA predictions erroneously merge smaller islands into longer ones, or that CpGcluster erroneously fragments longer islands into many smaller ones. Next, we use alternative TSSs and single CpG resolution methylation data to shed light on these questions.
Frequently, CpGcluster predicts more than one island within the promoter region. It has been shown  that 37.8% of all RefSeq genes have more than one CpGcluster island, while only 3.2% have more than one TJ island. Following the premise "one promoter one CpG island", this observation was interpreted as a disadvantage of CpGcluster . However, in recent years, new insights into the regulation of gene expression became available, showing among other things a frequent use of alternative TSSs. The existence of alternative TSSs opens the possibility that more than one island per gene might exist. Therefore, the high percentage of genes with more than one CpGcluster island might instead indicate a more specific relation of CpGcluster islands to alternative promoters or TSSs. To check this possibility, we used the DBTSS database . Out of 15,194 RefSeq genes annotated in the latest DBTSS release, 7,895 (52%) have at least one alternative TSS. With such scenario, one might expect up to 52% of all promoters having more than one island in its promoter (one for each TSS). Given these numbers, the reported 37.8% of genes with more than one CpGcluster island might look not so inadequate.
Co-localization of CpG islands and alternative promoters.
Numbers of overlapping islands
All the TSSs
In the human genome, there are a total of 166 bi-directional promoter pairs which share one long SWA CGI but two separated CpGcluster CGIs. The gene-pair shown in Figure 7 may be just an example of extreme differentiation in gene-expression: while the first member of the gene-pair is a housekeeping gene, the second one is a tissue-specific gene. However, one cannot reasonably expect that this may be the rule for all the bidirectional gene-pairs. In fact, after analyzing the expression profiles in a sample of 73 healthy tissues, only 16 (or 9.64%) gene-pairs show a completely divergent pattern of gene-expression (coexpression value ≤ 0.2, see Methods), while 13 (or 7.83%) exhibit complete coexpression (coexpression value = 1). The remaining gene-pairs show intermediate values of coexpression.
On the other hand, by using single base resolution methylation data , we also analyzed methylation differences between the CGIs overlapping bi-directional promoters. We found that 10 (or 11.24%) of these island-pairs in H1 stem cells, and 15 (or 16.85%) in the IMR90 fetal lung fibroblasts, show significant differences (Mann-Whitney non-parametric test) in their methylation average (p ≤ 0.05).
A functional CpG island should show a rather homogenous methylation profile among the different CpGs and over the different tissues. For example, the existence of more than one methylation domain within a predicted island might indicate an erroneous merging of two small islands into a single longer island.
CpG "islets", genomic regions not conventionally classified as CGIs because of their short length (<200 bp), but having a GC content and observed-to-expected CpG ratio characteristic of a CGI, have recently been identified in a 6.76-Mbp chromosomal region (10q25) containing a neocentromere . Some of these islets remain unmethylated, corresponding to sites of active transcription and/or boundaries that separate major chromatin sub-domains. This suggests that, as conventional islands, the islets can also participate in the maintenance of a particular genomic pattern of methylated/unmethylated CpGs, thus contributing to the differential regulation of gene expression [3–5].
Overlap of CpG islets (N = 88,137) with different sets of promoters and evolutionarily conserved elements.
Number of overlapping CpG islets
Number of overlapping CpG islets exclusively predicted by CpG cluster
Promoters from RefSeq database
TSSs from DBTSS database
Promoter regions from DBTSS database
Number of unmethylated and differentially methylated CpG 'islets'.
Number of CpG islets
CpG 'islets' exclusively predicted by CpGcluster
HEP (12 tissues)**
Lister et al. 2009 (2 cell lines)***
We systematically compared conventional SWA for detecting CGIs to a clustering method, namely the CpGcluster algorithm. We showed than both approaches perform very similar when predicting long, unmethylated regions or polycomb sites. However, we found three scenarios where the CpGcluster algorithm seems to have advantages. First, the statistical significance assigned to each CpGcluster island seems to be a better quality parameter than the window size of conventional finders, as it reduces more efficiently false positive predictions. Second, we have shown that CpGcluster islands co-localize in a more specific way to alternative TSSs and methylation domains. Third, we have shown that many of the small islands predicted by CpGcluster might be functional, given the overlap with conserved elements or promoter regions. Moreover, 30% of the differentially methylated islets are exclusively predicted by CpGcluster, which suggests this method as the option of choice for the experimental verification of islet functionality.
We used human genome assembly NCBI 36.1 (hg18), downloaded from the UCSC genome browser http://hgdownload.cse.ucsc.edu/downloads.html#human.
To quantify the co-localization of the predictions with promoter regions and principal transcription start sites we used the RefSeq gene annotation . We furthermore used the DBTSS database version 6.0 , as it annotates also alternative transcription start sites, as well as start sites which cannot be assigned to a known RefSeq transcript. From both, the RefSeq and DBTSS annotation, we extracted the coordinates of two regions; the transcription start site (TSS) and the promoter regions, defined as TSS-1500 bp to TSS+500 bp.
We determined the overlap of CGIs with conserved elements (PhastCons) and spurious Alu elements. The evolutionarily conserved elements  and the RepeatMasker  annotation of repeated elements where downloaded from the UCSC table browser . In general, we consider two measures to quantify the overlap between CGIs and genomic elements. First, we define the mean coverage of a CGI prediction as the mean value of all coverage fractions. The coverage fraction can be calculated as the number of bases of an island corresponding to a given genomic element divided by the island length. Furthermore, we calculate the overlap fraction as the number of islands which overlap in at least one base with a given genomic element divided by the total number of predicted islands.
For SWA CGI finders, a CpG island was at least 200 bp long, which excluded the detection of any shorter tracts. To detect CpG-rich regions, disregarding its length, we used a recently published CpG island finder algorithm (CpGcluster, ) which does not rely on any length threshold but directly predicts statistically significant CpG clusters. Briefly, the CpGcluster algorithm can be divided into two steps. First, based on a distance threshold, the individual CpGs which are below this threshold are clustered along the DNA sequence. Second, by means of the negative binomial distribution a p-value is assigned to each CpG cluster, which allows the prediction of highly significant clusters such as CpG islands.
We considered five computational predictions of CpG islands. For the CpGcluster algorithm , we generated two prediction sets by setting the assigned p-value to two different thresholds. We generated a relaxed set with p-value <= 1E-5 and a strict set by setting the threshold to 1E-20. We implemented the TJ algorithm, as explained in , by setting the thresholds to: length ≥ 500 bp, GC content ≥ 55%, ObsCpG/ExpCpG ≥ 0.65 and minCpG >= 0.6*Lisland/16 (to avoid "mathematical" islands). We generated the CpGproD prediction  running the program http://pbil.univ-lyon1.fr/software/cpgprod.html with default parameters. Finally, we downloaded the UCSC CpG island predictions from the UCSC table browser .
We used the GeneAtlas2 expression data  to determine the co-expression of gene pairs sharing a bi-directional promoter. The "coexpression value" for a couple of genes is the ratio of the number of tissues in which both genes are simultaneously expressed (signal levels > 200) or simultaneously not expressed (signal levels <= 200), and the number of healthy tissues with expression data.
Since the lack of methylation of a CpG island is a very good indicator of function , we used several different sources of experimental methylation data. Weber et al.  detected methylation states in two different tissues, fibroblast and sperm. We extracted 13,277 non-overlapping regions which are unmethylated in at least one of the two tissues (scaled 5mC log2 ratio < 0.3). Next, we used 17,383 CpG island recently detected in blood cells by means of a new technique .
Finally, we assigned methylation states (unmethylated, methylated and differentially methylated) to our CpGcluster predictions by means of the data from the HEP-human epigenome project . The data comprises about 1.9 million CpG methylation values, obtained from the analysis of 2,524 amplicons across chromosomes 6, 20 and 22 in 43 samples (derived from 12 different tissues). We first calculated the mean methylation of each CpG dinucleotide over the different clones, then deleting all CpGs which have been detected in less than 2 clones or in less than 6 different tissues. Subsequently, the individual CpGs were labeled as methylated (mean methylation >= 80), intermediate methylated (80-20) and unmethylated (under 20) for each of the different tissues. Next, we define the methylation states of the CpGs over the different tissues in the following way: i) methylated CpG: methylated in more than 50% of tissues and never unmethylated, ii) unmethylated CpG: unmethylated in more than 50% of tissues and never methylated, iii) differentially methylated CpG: both, methylated and unmethylated in different tissues, the number of intermediate methylation states being smaller than 50%. Finally, we assign a methylation label to the CpG islands which have methylation data for more than 50% of its CpGs: i) methylated: more than 50% of the CpGs are methylated and no unmethylated CpG exist, ii) unmethylated: more than 50% of the CpGs are unmethylated and no methylated CpG exist, iii) differentially methylated: more than 50% of all CpGs need to be differentially methylated.
We also used the sequence reads obtained by MethylC-Seq for two human cell lines , H1 human embryonic stem cells and IMR90 fetal lung fibroblasts, to get the average methylation level of single cytosines at both DNA strands for these two methylomes. All islands need more than 50% of its CpGs covered. Only cytosines covered by at least 10 reads were counted.
When comparing the prediction of CpG islands to a gold standard (e.g. experimentally verified islands), we define:
True Positives (TP): An island overlapping in at least 1 bp with the gold standard
False Positives (FP): An island not overlapping with the gold standard
False Negative (FN): An island in the gold standard that has not been predicted.
Note that we consider all islands not overlapping with the gold standard as false positive predictions. However, no complete gold standard exists, and therefore an unknown number of these islands will be actually true positive predictions. This assumption does not affect the sensitivity, as FP does not occur in the equation, but it affects the PPV. Consequently, and since the PPV can only increase when some FPs turn out to be TPs, the value used in this work is the lower boundary PPV of the prediction, e.g. the worst case scenario when all islands which do not overlap with the gold standard are indeed false positives.
Ratio between observed and expected CpG frequencies
%G + C: Molecular fraction of guanine and cytosine
Phylogenetic Conserved Elements
The sensitivity of the prediction
Positive Predictive Value of the prediction
Takai/Jones program or island
Transcription Start Site
We thank Andy Choo and Nicholas Wong, from the University of Melbourne (Parkville, Victoria, Australia), by sharing with us their data on CpG islets. We acknowledge the Spanish Government (Grant No. BIO2008-01353) and the Spanish Junta de Andalucía (Grant Nos. P06-FQM1858 and P07-FQM3163) financial support.
MH acknowledges financial support from the 'Juan de la Cierva' grant from the Spanish Government. GB acknowledges financial support from the 'Programa de formación de investigadores del Departamento de Educación, Universidades e Investigación' grant from the Basque Country Government.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.