Annotation of gene promoters by integrative data-mining of ChIP-seq Pol-II enrichment data

Background Use of alternative gene promoters that drive widespread cell-type, tissue-type or developmental gene regulation in mammalian genomes is a common phenomenon. Chromatin immunoprecipitation methods coupled with DNA microarray (ChIP-chip) or massive parallel sequencing (ChIP-seq) are enabling genome-wide identification of active promoters in different cellular conditions using antibodies against Pol-II. However, these methods produce enrichment not only near the gene promoters but also inside the genes and other genomic regions due to the non-specificity of the antibodies used in ChIP. Further, the use of these methods is limited by their high cost and strong dependence on cellular type and context. Methods We trained and tested different state-of-art ensemble and meta classification methods for identification of Pol-II enriched promoter and Pol-II enriched non-promoter sequences, each of length 500 bp. The classification models were trained and tested on a bench-mark dataset, using a set of 39 different feature variables that are based on chromatin modification signatures and various DNA sequence features. The best performing model was applied on seven published ChIP-seq Pol-II datasets to provide genome wide annotation of mouse gene promoters. Results We present a novel algorithm based on supervised learning methods to discriminate promoter associated Pol-II enrichment from enrichment elsewhere in the genome in ChIP-chip/seq profiles. We accumulated a dataset of 11,773 promoter and 46,167 non-promoter sequences, each of length 500 bp, generated from RNA Pol-II ChIP-seq data of five tissues (Brain, Kidney, Liver, Lung and Spleen). We evaluated the classification models in building the best predictor and found that Bagging and Random Forest based approaches give the best accuracy. We implemented the algorithm on seven different published ChIP-seq datasets to provide a comprehensive set of promoter annotations for both protein-coding and non-coding genes in the mouse genome. The resulting annotations contain 13,413 (4,747) protein-coding (non-coding) genes with single promoters and 9,929 (1,858) protein-coding (non-coding) genes with two or more alternative promoters, and a significant number of unassigned novel promoters. Conclusion Our new algorithm can successfully predict the promoters from the genome wide profile of Pol-II bound regions. In addition, our algorithm performs significantly better than existing promoter prediction methods and can be applied for genome-wide predictions of Pol-II promoters.


Background
Alternative promoter usage is known to affect more than half of all human and mouse genes, and has been proposed as a primary driver of varied transcriptional regulation in different cellular conditions or developmental stages [1][2][3][4]. Numerous genes displaying complex transcriptional regulation, because of the use of alternative promoters, have been studied thoroughly [5]. Recent annotations of the human and mouse genomes suggest that many differentiation and disease-associated genes contain alternative promoters. Annotation of all human and mouse gene promoters that are differentially used during development, in different cell/tissue types or aberrantly activated in disease conditions is still incomplete and is essential for defining the transcriptome and proteome of human and mouse genomes.
The development of chromatin immunoprecipitation methods coupled with DNA microarray (ChIP-chip) technology and massively parallel sequencing (ChIPseq) has enabled genome-wide identification of promoters, using antibody against RNA polymerase II (Pol-II) in different cells or tissues [6,7]. The combined signatures of RNA Pol-II binding and histone modification marks like H3K4me3, H3K9Ac obtained by these high throughput technologies are being used to identify human and mouse transcriptional units [8]. However, there are some challenges in predicting promoter usage based on the enrichment regions/peaks observed in these ChIP-chip/seq experiments. The ChIP-chip/seq technology requires antibodies with extremely high affinity and specificity for the target transcription factors. Unfortunately, such antibodies are not available for most human transcription factors, including Pol-II, producing nonspecific enrichment in the ChIP-chip/seq profiles. The non-specific enrichment regions could be eliminated from the analysis by performing a ChIP-chip/seq experiment on the same cell or tissue lacking the specific factor. However, in most cases this is not feasible and we have to rely on other methodologies like the use of nonspecific IgG ChIP-chip/seq to decrease the non-specific enrichment background. The major challenge in annotating promoters based on RNA Pol-II enriched regions/ peaks is the spread of the transcribing polymerase throughout the gene and as a result all genomic regions bound by RNA Pol-II are enriched in these experiments, producing significantly large number of enriched peaks after the initial statistical analysis [9]. Though the initiator form of RNA polymerase II (phosphorylated CTD at Ser5) is enriched at a higher level in promoter region of actively transcribed genes, it is not restricted to the promoter region. Moreover, the promoters with stalled RNA Pol-II do not show an enrichment for the Ser5 phosphorylated form of RNA Pol-II [10]. Similarly, the histone marks namely H3K4me3 and H3K9Ac, which are highly enriched in promoter regions, are not exclusively present in promoter regions [8,11]. Currently it is not possible to identify promoters with high confidence based on RNA Pol-II ChIP-chip/seq enrichment data alone, thus warranting development of better classification algorithms for accurate identification of promoter related Pol-II enriched regions.
Here, we developed a computational method to discriminate promoter associated RNA Pol-II enriched regions of length 500 bp from the enrichment at other genomic regions, using the rich source of existing promoter data and associated chromatin modification signatures and various DNA sequence features. We prepared a data-set consisting of 11,773 Pol-II enriched promoters and 46,167 Pol-II enriched non-promoter regions from our recent ChIP-seq experiments, using antibody against Pol-II on five mouse tissues. We systematically trained and evaluated recent ensemble classifiers on this data set, using both 10-fold cross validation and testing on independent test set, and selected Bagging and Random Forest classifier for the final algorithm.

Methods
Dataset of Pol-II enriched promoters and non-promoters for training the classification models The training set was generated from RNA Pol-II ChIP-seq data of five mouse tissues (Brain, Kidney, Liver, Lung and Spleen) generated by our lab. The RNA Pol-II ChIPseq data was first processed and Pol-II enrichment peaks were identified at an FDR of 0.001, by assuming that peaks in the random background would follow Poisson statistics. The identified peaks were compared with TSS of non-redundant gene list generated from four different sources: RefSeq, Vega, Ensembl and UCSC. The gene lists were downloaded from UCSC genome browser [12] for mm9. Any peak that falls within -300 bp to +200 bp of known TSS (from compiled non-redundant TSS) were taken as promoter peak. The Pol-II peaks which are inside a gene but do not overlap with any of the known TSS are candidate non-promoter peaks. For our negative set, we consider only those peaks which fall within transcripts that possess a Pol-II enriched peak at the corresponding promoter (-300 bp to +200 bp). Also, any peak which falls within the promoter region of homologous gene known in other organism and within the 5' end of compiled set of non-redundant expression sequence tags (ESTs) was removed from the negative set, because some of those could be undiscovered novel promoters in the mouse genome. The homologous gene track (xenoRefGene track) was downloaded from UCSC genome browser. After identification of promoter and non-promoter peaks, the Pol-II peak was annotated as actual TSS for the transcript in that particular tissue. For each annotated peak (in both promoter and nonpromoter sets), sequences were generated each of length 500 bp (-300 bp to +200 bp around the peak). After performing this procedure for all the tissues, we again compared the records in promoter and non-promoter sets with each other and removed the overlapping records from the non-promoter set. Finally, our complete dataset has 11,773 records in the promoter set and 46,167 records in the non-promoter set. We used 8,793 and 34,686 records from promoter and non-promoter sets respectively for training the classification models. In addition to 10-fold cross-validation, we used the remaining records (2,980 promoter and 11,481 nonpromoter cases) to test the performance of the fitted models. The data sets are available as supplementary information at [13].

Classification models
We tried different state-of-art ensemble and meta classifiers for identification of promoter and nonpromoter classes. The WEKA data-mining toolbox [14] was used for building the classification models. The different classifiers tested on 39-dimensional feature vector are: LogitBoost [15], Bagging [16], Rotational Forest [17] and Random Forest [18]. The detailed description of the classification methods is provided in Supplementary Method.

Feature variables used in classification model
Each sequence record (500 bp window) of the promoter and non-promoter set is represented by a 39-dimensional feature vector. The feature values were calculated based on (a) DNA sequence composition; (b) DNA physico-chemical-structural properties, and (c) experimental data. Most of the conversion tables that are based on DNA physical-chemical-structural properties were downloaded from [19]. The feature variables used in the classification models are briefly described below.
Let a given DNA sequence be: (a) Properties based on DNA sequence composition We calculate 10 different features in this category. The first 7 features are calculated from single nucleotide composition. Let n x represents the total number of times symbol x appeared in S and x Δ 1 for single nucleotide features.
1. A_Fraction: n A /L 2. C_Fraction: n C /L 3. G_Fraction: n G /L 4. T_Fraction: n T /L 5. PurPyr_Fraction: (n A + n G -n C -n T )/L 6. AmKe_Fraction: (n A + n C -n G -n T )/L 7. WeSt_Fraction: (n A + n T -n C -n G )/L The remaining 3 features are related to CpG island. One of the features is based on di-nucleotide composition, where x Δ 2 . And remaining 2 features are based on tri-nucleotide composition, where x Δ 3 . Similar CpG features were used in [19] for promoter prediction. 8. CpG1: (2*n CG + 2*n GC )/(L-1) 9. CpG2: (n ACG + n AGC + n CAG + n CCG + n CGA + n CGC + 2* n CGG + n CGT + n CTG + n GAC + n GCA + 2* n GCC + n GCG + n GCT + 2* n GGC + n GTC + n TCG + n TGC )/(L-2) 10. CpG3: (4* n CAG + n CCG + n CGG + 4* n CTG + 4* n GAC + n GCC + n GGC + 4* n GTC )/(L-2) (b) Properties based on physico-chemical-structural property of DNA sequences In this category we calculate 22 features. Let p (x) represent a mapping function for a property 'P', where x Δ 1 or x Δ 2 or x Δ 3 or x Δ 4 depending upon given property. The feature value for a given sequence 'S' based on property 'P' is given by f P , where n x represents total number of times symbol x has appeared in S. And Δ ≡ Δ 1 or Δ ≡ Δ 2 or Δ ≡ Δ 3 or Δ ≡ Δ 4 depending up on the property 'P'. The performance of the fitted models and other best performing promoter prediction programs, which are publicly available for download to run on whole chromosomes on a local computer, was tested on an independent un-seen data set. When a predicted promoter overlaps with Pol-II enriched promoter, then the respective record is counted as true positive (TP). And when such case is missed it is termed as false negative (FN). When a predicted promoter overlaps with Pol-II enriched non-promoter, then such case is counted as false positive (FP). And if such case is predicted as nonpromoter then it is termed as true negative (TN). The performance of classifier is evaluated based on the promoter prediction metrics suggested by Bajic et. al.

Results
Classification models to predict promoters using chromatin modification signatures and DNA sequence features For selecting the best performing classifier, we trained four different ensemble and meta classification models on 3/4 th of the dataset and tested on the remaining 1/4 th of the dataset. The performance measures obtained by using 10-fold cross-validation and independent test set are presented in Tables 1 and 2. The performance measures, in terms of sensitivity and positive predictive accuracy, among the classifiers do not vary much over the four different models. Bagging and Random Forest models are slightly better than the other two models, showing lower error rates. Figure 1 allows for a more informative discussion on the relative predictive performance of the models. It is clear that Bagging, LogitBoost and Random Forest perform more or less similar and slightly better than Rotational Forest, with overall positive predictive value greaten than 95 and correlation coefficient greaten than 0.9. We then implemented the classification models given by Bagging and Random Forest methods in our algorithm which, will be applied to scan all the Pol-II enriched peaks in the mouse genome.
While the classification methods used here are considered "black box" methods, with no interpretable classification model, the methods still provide useful BMC Bioinformatics 2010, 11(Suppl 1):S65 http://www.biomedcentral.com/1471-2105/11/S1/S65 information, such as variable importance. One of the measures of variable importance in Random Forest method is the mean decrease in accuracy, calculated using the out-of-bag sample. The difference between the prediction accuracy on the untouched out-of-bag sample and that on the out-of-bag sample permuted on one predictor variable is averaged over all trees in the forest and normalized by the standard error. This gives the mean decrease in accuracy of that particular predictor variable which has been permuted. Thus, the importance of the predictor variables can be ranked by their mean decrease in accuracy. Figure 2 shows the list of feature variables ranked according to mean decrease in accuracy of classification. It is interesting to note that feature variables based on experimental data such as CAGE tags, Pol-II enrichment, and H3K4me3 enrichment rank among the most discriminative variables from mean decrease in accuracy graph (Figure 2).

Comparison with other promoter prediction programs
We compared our algorithm with some of the existing best performing promoter prediction methods [44]:   EP3 [19], Eponine [45], FirstEF [46], ProSOM 2.5 [47]. Table 2 and Figure 3 show that our classification model out-performs these existing programs based on independent (unseen) test set of Pol-II enriched promoters and Pol-II enriched non-promoters.
Annotation of promoters in mouse genome using Pol-II ChIP-seq data Although extensive promoter annotations are available from the EBI and UCSC genome servers, most annotations do not contain information about tissue or  cell-type information from experimental data. To demonstrate the efficacy of our new algorithm in finding promoters and to provide the annotation of potential novel promoters, we used the Pol-II enrichment peaks obtained from seven ChIP-seq datasets available on in vitro adipocyte differentiation of mouse 3T3-L1 cells and mouse ES cells, and the results are presented in Table 3. We applied the Random forest classification model (built using the training set) on a sequence of length 500 bp around each Pol-II enriched peak (-300 to +200 bp) in both strands. If 500 bp region from both strands are predicted as peak then they are merged and counted as one promoter. For each predicted promoter region, we annotated it to a nearby gene, if the Pol-II enriched peak is located within -2 Kbp to +500 bp around the corresponding gene TSS. We generated a non-redundant TSS file for coding genes from Refseq, Vega, ensembl and UCSC genes. For non-coding genes we used information available at RefSeq, Vega, ensembl, UCSC, miRBase [48] and recently discovered large non-coding RNAs (lincRNA) [49]. Table 3 presents the total number of peaks predicted by the model as promoters and also the number of annotated coding and non-coding genes in each sample. We then combined the results from all seven samples of predicted promoters in order to identify alternative promoters for each gene. For protein coding gene set, we found that there are 13413, 5064, and 4865 genes with one promoter, two promoters, and three or more promoters respectively. In other words, based on these annotations, 42.5% of the protein coding genes in mouse genome have two or more alternative promoters. For non-coding genes, we found that 4757, 1181, and 677 genes with one promoter, two promoters, and three or more promoters respectively.

Future directions
We will use this program to annotate human gene Pol-II promoters by running on all the publicly available ChIPseq Pol-II enrichment profiles. Our method successfully predicts 500 bp promoter regions (-300 bp to + 200 bp) and to better localize the core-promoter regions within the predicted promoters, we will apply recently developed CoreBoost_HM program published by Zhang Laboratory at CSHL [50].

Discussion
Chromatin modification and transcription factor binding profiles in the mammalian genomes is rapidly accumulating with the advent of next generation sequencing approaches. However, computational methods to effectively integrate these profiles to identify and annotate the promoter usage in specific cell/tissue types or developmental stages, are still limited. Recently, machine learning strategies have been applied to combine some of the wealth of published ChIP-seq data sets, such as chromatin modification signatures, to predict core promoter regions [50]. A logical step in analyzing the Pol-II enriched genomic regions is to scan those regions by existing promoter prediction methods to predict whether the enriched region is a promoter or non-promoter. However, we found that the performance of the existing methods is not satisfactory, and we speculate that the training set used in building the classifier was mostly responsible for their poor performance. We, therefore, build a bench-mark data set of Pol-II enriched promoters and Pol-II enriched nonpromoters to train the classifiers, which shows significant improvement over the existing programs. Theoretical and empirical works using classification, regression trees, variable selection in linear and non-linear regression have shown that bagging and ensemble based methods can generate substantial prediction gain. In fact, based on the evaluation of 10-fold cross validation and testing on an independent data set, we found that both Bagging and Random Forest methods performed with highest accuracy (better than 95% prediction accuracy).

Conclusion
In conclusion, we have developed a novel algorithm based on Bagging and Random Forest based classification methods to predict Pol-II bound promoters from ChIP-seq profiles. The present algorithm will help the discovery of novel promoters and ongoing annotation of alternative promoters of human and mouse genes from different ChIP-seq experiments.