BACKGROUND
transcription factors  serve as the final molecules in signal transduction pathways that coordinate expression of target genes.
when activated in response to upstream signals, often encoded as chemical ligands and protein modification, tfs bind to their cis-regulatory sites to exert their regulatory effects on their target genes.
during the process, tfs often interact with other proteins, which further modulate the function and efficacy of tfs to achieve fine-tuned regulation of gene expression; studying such interactions and regulations is an increasingly important component of studying gene expression systems.
nuclear receptors , such as estrogen receptor α , are transcription factors that migrate to the nucleus  to regulate downstream target genes.
nrs play important biological roles in normal physiology and disease.
in particular erα plays an important role in both breast cancer and osteoporosis.
upon ligand binding, erα and other nrs are bound by proteins called co-regulators that recruit transcriptional machinery and chromatin modifying enzymes.
co-regulators are therefore critical in nr activity.
understanding the composition of functional nr/co-regulator complexes in specific signaling contexts could provide a basis for the development of novel nr- and co-regulator-targeted therapeutics.
the problem addressed in this paper arose from a study of the interaction between the major erα co-activator src- <dig> , also known as ncoa <dig>  with erα and the impact of such interactions gene expression  <cit> .
recently, chromatin immunoprecipitation coupled with high-throughput next-generation sequencing  has become the main technology for global characterization of the transcriptional impact of nrs and their co-regulators  <cit> .
chip-seq involves the short-read  sequencing of the chip-enriched dna fragments.
these short sequence reads  are then aligned to a reference genome.
then the actual binding loci from the positional tag distributions  are determined using 'peak calling' algorithms.
numerous peak calling algorithms have recently been developed for identifying chip-enriched genomic regions from chip-seq experiments  <cit>  but there is a wide range of discordance among the peak calls from different algorithms  <cit> .
therefore, there is a need for the methods that can integrate additional information besides chip-seq tags to identify functional tf binding sites.
furthermore, studying the interactions between tfs and their co-regulators through chip-seq technology poses an additional challenge since co-regulators do not directly bind dna.
co-regulator chip-seq measures the secondary protein-dna binding through primary tfs and leads to relatively weak sequencing signals--i.e.
relatively small number of sequence tags above noise.
as such, it remains a challenge for contemporary peak calling methods to detect weak secondary protein-dna-binding signals and simultaneously maintain a high specificity.
often, a well-designed experiment studying interaction between a tf and its co-regulator generates critical information in addition to the chip-seq data for the co-regulator binding.
for example, chip-seq data reflecting the binding of the primary tf of interest to its cis-regulatory sites are often collected; the genomic sequence surrounding chip-seq peaks are usually available, which can be used to reflect the intrinsic sequence characteristics of regulatory sites; transcriptomic data that reflect functional outcomes of the interaction of the tf and its co-regulators can also be monitored.
in this study, we investigated and compared different statistical and machine learning approaches to integrate multiple types of information to overcome the difficulty of identifying functional erα/src- <dig> interaction in presence of weak chip-seq signal.
RESULTS
the biological study underlying this paper aims to investigate the impact of erα/src- <dig> interaction on estrogen induced gene expression in a bone cell line transfected with erα , which may shed light on the effect of estrogen-related bone development, bone loss, and potentially bone metastasis.
we have generated chip-seq data using anti-erα and anti-src- <dig> antibodies in presence and absence of estradiol .
to further investigate the impact of interactions between this nr//co-regulator pair, we collected expression array data from the same cell lines with a combination of e <dig> treatment and src- <dig> knock down.
in general, the results of an src- <dig> chip-seq experiment would reflect secondary, indirect binding of src- <dig> to dna through multiple nrs.
however, in this study our experimental design aims to investigate specifically estrogen-induced interactions between erα, src- <dig>  and dna.
the detailed results of the experiments are being prepared for a separate publication .
in the current paper, we address the fundamental issue of identifying reliable and functional erα/src- <dig> dna binding sites.
identifying src- <dig> binding sites based on anti-src- <dig> chip-seq
we first set out to investigate the efficacy of studying erα/src- <dig> dna binding sites only based on the results of chip-seq experiments performed using an anti-src- <dig> antibody.
potential src- <dig> binding peaks were identified using three different algorithms: macs  <dig> .
<dig>  <cit> , bayespeak <cit> , and t-pic  <cit> .
table  <dig> shows the number of peaks identified by the above algorithms with different cut-off thresholds and the corresponding number of genes to which the peaks are mapped.
* union and intersection of the peaks by the three methods, as shown in figure  <dig>
the results of the table raise the following issues during interpretation: first, as expected, applying different cut-off thresholds to the results by a given algorithm leads to a different number of peaks being identified: there is inevitably a trade-off between the number of peaks recovered and the quality of the peaks.
second, different algorithms assess the quality of the peaks based on different assumptions and methodologies: there is no consensus on the "goodness" of quality scores of these algorithms.
we further noted that different versions of a same algorithm return different quality scores.
finally, different algorithms return disjointed sets of peaks, as shown in figure  <dig>  indicating that distinct assumptions and approaches enable an algorithm to discover some potential peaks that evade detection by other algorithms.
these issues force decisions that potentially impact the conclusions of the study such as: which algorithm performs better, what cut-off threshold for a given algorithm to pick, and how to consolidate the results from different algorithms so that one can maximize the number of high quality peaks.
making these choices remains challenging due to the lack of consensus in the field  <cit> .
in order to compare our results with a recently published study by lanz et al.
<cit> , which analyzed dna recruitment of the co-regulators src- <dig>  we studied the estrogen-induced src- <dig> peaks identified from our data using macs algorithm with a cut-off threshold of p value at 10- <dig>  a threshold based on their study.
our analysis yielded a total of  <dig>  peaks, which were further mapped to  <dig> genes.
the number of peaks identified by us with the above condition is far fewer compared to their study.
the discrepancy is likely to be, at least in part, due to amplification of src- <dig> in mcf- <dig> cells used in their study and possible differences in antibody affinities.
however, these results also raised the hypothesis that chip-seq signals of secondary binding at physiologic levels are usually weaker.
this suggests that conventional cutoff thresholds for peak calling algorithms may be too stringent, neglecting weak peaks  potentially resulting from real erα/src-1/dna interactions.
therefore, additional information besides src- <dig> chip-seq tags should be capitalized to enhance identification of functional binding sites.
integrating multiple sources of biological information for identifying src- <dig> binding sites
to corroborate the results of src- <dig> chip-seq, we also studied the erα chip-seq data  and investigated peaks overlapping between the erα and the src- <dig> chip-seq results.
by varying the cut-off threshold of macs, we identified different numbers of overlapping peaks between erα and src- <dig>  with the number of overlapping peaks increasing as the cut-off threshold relaxes .
the results again indicated that the conventional cut-off thresholds are failing to identify putative erα/src- <dig> dna binding sites .
on the other hand, simply relaxing the cut-off threshold is likely leading to increased false positive peak calls.
thus a principled method is needed to further identify functional erα/src- <dig> dna binding sites.
to elucidate functional erα/src- <dig> dna binding events, i.e., the binding events that influence gene expression, we generated and analyzed expression array data from cells that were treated with vehicle and e <dig> in the presence and absence of src- <dig> sirna .
the microarray data enabled us to identify  <dig> genes whose expression response to estrogen treatment required src- <dig>  which are hereto referred to as src-1-sensitive genes.
when we compared the list of src-1-sensitive genes and the list of genes with src- <dig> binding sites derived by macs at p value = 10- <dig>  we noted that only  <dig> genes overlapped among the lists.
while the discrepancy between the number of src-1-sensitive genes and genes with src- <dig> peaks could be explained by other biological factors, such as that many src- <dig> peaks were not functional or secondary expression effects, it also supported our general hypothesis that erα/src- <dig> interaction chip-seq signal is relatively weak and potentially true functional erα/src- <dig> dna binding sites were missed by the stringent setting of the peak calling algorithm.
while it may be tempting to directly combine the information from src- <dig> chip-seq, erα chip-seq and microarray data by identifying the intersections of overlapping genes and peaks, such an approach is overly simplistic and ignores other potentially informative data, e.g., the genome-sequence characteristics of erα/src- <dig> interaction sites and the prior information of known erα/src- <dig> interactions.
these considerations motivated us to investigate and compare different principled machine learning approaches in order to improve the sensitivity and specificity of detecting erα/src- <dig> dna binding sites by integrating multiple types of information.
an integrative approach to detect erα/src- <dig> dna binding sites
the overall framework and rationale of our information integration approaches are as follows.
we formulated the task of identifying functional erα/src- <dig> dna binding sites as a classification task, in which we performed a binary classification to label a potential src- <dig> binding site derived from chip-seq analysis as either functional or nonfunctional.
we investigated both supervised learning, which allows us to take advantage of existing knowledge of erα/src- <dig> interactions, and unsupervised learning, which allowed us to take an unbiased approach.
the classification formulation allowed us to pool more candidate peaks identified by different peak calling algorithms at relaxed cutoff thresholds so that we did not have to rely on a single "best" algorithms and "optimal" parameterization but resorted to our classification to identify functional erα/src- <dig> dna binding.
in this study, we collected the union of the peaks returned by all three algorithms at the cutoff threshold as follows, macs: p value cutoff 10- <dig>  bayespeak: posterior probability  ≥  <dig>  and t-pic: p value cutoff 10- <dig>  this led to a pool of  <dig>  candidate peaks.
another important advantage of the classification approach is that it allows us to integrate multiple types of biological information collected from our experiments and public databases by representing them as features for a classifier.
in this way, multiple types of information contribute to the classification of potential peaks and their impact can be determined by learning algorithms.
for each candidate src- <dig> peak, we constructed the following features: a vector of binary features representing presence/absence of nucleotide trigrams , which reflects the intrinsic characteristics of genome sequence surrounding the summit of a peak region; an average of predicted nucleosome-occupancy scores, which represents the chromatin structure characteristics around the peak summit; a binary feature reflecting if a primary binding peak, i.e., the erα chip-seq peak, overlaps with the src- <dig> peak; and a binary feature representing the functional outcome of erα/src- <dig> interaction, i.e., whether the peak is mapped to an src- <dig> sensitive gene.
detailed descriptions of features are presented in methods section.
we evaluated the results of predictions from classification algorithms by determining if conserved erα binding motif can be found in the classified peaks, as an indication that a peak is the result of erα/src- <dig> dna binding.
searching for instances of conserved tf binding motifs at the predicted binding loci is considered the most prominent verification method for validating peaks  <cit> .
unsupervised classification
first, we explored if the candidate peaks could be divided into two distinct groups by unsupervised learning in an unbiased manner.
we applied a k-means clustering procedure to the data and the results are listed in table  <dig>  we inspected the genome sequences of the peaks to assess if a conserved motif for estrogen response element  was detected by a motif classification algorithm referred to as clover <cit> .
in cluster 46% of  <dig>  peaks contained the ere motif, and, in cluster 27% of  <dig>  peaks contained the ere motif.
we believe this is not a good separation of the peaks in that, even though cluster  <dig> has more ere-containing peaks, it is a bigger cluster and only 46% of peaks contain eres.
thus the results would likely lead to a high false positive rate with respect to src- <dig> binding.
supervised classification
supervised learning requires labeled data as training cases.
obtaining a training set through large-scale experimental validation of erα/src- <dig> dna binding is costly and difficult to perform.
therefore, we investigated whether a relatively small amount of labeled data in the supervised learning task would result in better separation compared to unsupervised clustering with our feature set.
we experimentally validated  <dig> src- <dig> peaks as functional peaks by quantitative pcr  experiment .
we used these peaks as positive training cases, together with a set of randomly drawn control  chip-seq peaks as negative training cases, to train supervised classifiers.
we investigated the performance of three state-of-the-art classifiers: naive bayes , support vector machines  and random forest .
for nb and random forest classifiers, we set the classification thresholds at  <dig>  and  <dig>  respectively.
since the ratio of the positive and negative training cases may have impact on classification algorithms, e.g.
nb and svm, we explored using different ratios for training, between 1: <dig>  1: <dig>  and 1: <dig>  and classifiers were built from these training sets.
our test set consisted of all  <dig>  candidate peaks.
table  <dig> lists the total number of peaks in each class, the number ere-containing peaks in each class and the ratio reflecting ere enrichment.
results for classifiers with 1: <dig>  1: <dig>  and 1: <dig> training ratio  were very similar to each other .
therefore, just results for the 1: <dig> ratio were shown.
we noted that the supervised classification approaches have significantly increased the number of positive peaks when compared to those derived by peak calling algorithms based on recommended cutoff thresholds.
for example, nb returned  <dig>  positive peaks in comparison to  <dig>  and  <dig>  peaks returned by macs with cutoff p value set at 1e- <dig> and 1e- <dig>  which reflected a 6-fold and  <dig> -fold increase, respectively.
through further evaluation enrichment of ere in the genome sequences surrounding the peaks, we found that a similar percentage of peaks contained ere element: 69% for nb, and 72% and 66% macs at 1e- <dig> and 1e- <dig> respectively.
thus, the results indicate that the qualities of the positive peaks returned by nb were as good as those returned by the stringent peak calling in terms of ere enrichment.
we further inspected if the classification approaches retrieved additional functional peaks, i.e., the peaks that were mapped to src-1-sensitive genes derived from microarray experiment.
table  <dig> shows the results of the src- <dig> peaks returned by different peak calling approaches that were mapped to src-1-sensitive genes.
we noted that, by setting macs cutoff p values at 1e- <dig>  1e- <dig>  and 1e- <dig>  a total of  <dig>   <dig>  and  <dig> peaks were mapped to src-1-sensitive genes.
on the other hand, nb has identified  <dig> peaks that overlap with src-1-sensitive genes.
we also note that bayespeak exclusively identified some of the newly "discovered" functional peaks.
similarly macs also exclusively discovered some new peaks.
these results indicate that, based on different assumptions and criteria, different peak calling algorithms are capable of identifying potential peaks to complement other peak calling algorithms.
thus it is more sensible to consider candidate peaks from more than one peak-calling algorithm as long as an objective approach can be further applied to consolidate the results.
we also performed a 9-fold cross-validation experiment to assess if the algorithm can correctly identify the experimentally validated erα/src- <dig> training cases from the candidate peak pool.
results were listed in table  <dig>  nb classifier showed 89% precision, 100% recall and 96% accuracy, see methods section for the descriptions of the metrics.
this result increased our confidence that positive calls from our algorithm are likely to reflect real erα/src- <dig> dna binding.the results in the table indicate that the nb classifier performed better than the svm and rf classifiers, judging from relative enrichment of ere containing binding sites in the predicted positive peaks.
among the three classifiers tested in this study, the rf classifier performed worst.
we noted that the number of features that were used by rf during the learning was much smaller than the number of features utilized by other classifiers, which may partially explain the inferior performance of this algorithm in this experiment.
the svm method also performed worse on this task than probabilistic nb.
we conjecture that the reason might be that svm is complex algorithm with many parameters to adjust and therefore finding optimal parameters for decision boundary might be challenging for this task.
we therefore concentrated on the nb classifier because it could be readily used in both supervised and semi-supervised learning environment.
semi-supervised classification
our number of training cases is relatively sparse compared to a conventional machine learning setting.
semi-supervised approaches have been applied to conquer limitations of supervised and unsupervised methods when labeled data is scarce and obtaining large amounts of labeled data is expensive and time consuming  <cit> .
this is done by incrementally assigning instances, which are called with high confidence by a classifier, from unlabelled data into training cases in order to increase the number of training cases and thus enhance the generalizability of classification.
therefore, we investigated semi-supervised classification to see whether we could further increase the performance in identifying erα/src- <dig> dna binding.
in this study, we applied a self-training algorithm  <cit>  using nb as the base classifier because of its probabilistic outputs.
we iteratively assigned the most confident positive instances called by our classifier into training cases and found that performance of the self-training became stable after  <dig> iterations and stopped further training.
figure  <dig> shows the trends of percentage of ere-containing peaks in the positive calls in semi-supervised learning.
it is interesting to note that initially as a few pseudo-positive cases were imputed into the training cases the precision of the called positive peaks decreased but later became stable after  <dig> iterations.
a similar total number of peaks and the percentage of ere-containing peaks were identified by our semi-supervised learning algorithm when compared to other supervised learning experiments, see table  <dig>  thus the results do not show obvious advantage of semi-supervised learning over supervised learning algorithms in our experiment.
identification of informative features
biologically, it is of interest to identify the features that significantly contribute to the classification in that it will reveal the relationships between input features and outcome.
we rank key features by roc class reparability criteria using matlab   <cit> , using a training dataset containing the  <dig> true positive peaks and  <dig> random non-binding sites from igg peak calls.
following were the  <dig> top ranked features: "aac", peaks-mapped-to-src-1-dependent-genes, "gcg", "cgt", overlapping-with-erα-peak, "aag", "aca", "acc", "cgc", "aga", "agc", "agg", "cga", "act", "atc", see methods section for detailed descriptions of the features.
among the top-ranking features, we noted that the features reflecting the function outcome  and the interaction between er and src- <dig>  were ranked high, indicating the learning algorithm correctly recognized their importance in classification.
it is interesting to note that many nucleotide trigrams, which reflect the characteristics of sequences of peaks, were among the high-ranking features.
we aligned the top-ranking trigrams to the ere motif, as shown in figure  <dig>  indeed, the trigrams correspond well with the important components of the ere motif.
these results indicate that the classification learning algorithm, like the motif searching algorithms , is able to identify highly conserved "words" that constitute one of the important motifs of the training sequences.
we noted that the feature reflecting nucleosome occupancy at peak regions was ranked as 48th.
this may indicate either that nucleosome occupancy is a dynamic process and our static feature does not reflect the true occupancy status during the experiments or that the era/src- <dig> dna binding is not heavily dependent on nucleosome occupancy.
biological insights from improved peak calling
we further examined the impact of the improved src- <dig> peak calling on biological insights drawn from the dataset.
we conducted gene ontology analysis using the database for annotation, visualization, and integrated discovery   <cit>  on genes with an src- <dig> peak within  <dig> kb of the tss as determined by macs  or by our method.
we observed a dramatic difference in the identification of genes enriched in specific biological processes.
specifically, our method yielded in the calling of peaks in gene sets which were highly enriched for genes involved in blood vessel development  and actin filament-based processes .
indeed blood vessel development has previous been implicated in bone generation  <cit> .
further, within the genes enriched in these biological processes, we identified a number of genes with known functions in bone development.
since src- <dig> has already been implicated in e <dig> mediate bone maintenance, this observation provides evidence for the mechanism underlying this phenotype.
in contrast, genes with src- <dig> peaks determined by macs were not significantly associated with any biological processes.
CONCLUSIONS
we believe that the ability to improve chip-seq peak calling by utilizing available sources of biological information for indirect co-regulator binding in the presence of weak chip-seq signal is an important research area.
due to the intrinsic variability in the affinity of interactions between a tf and its co-regulators, it is inevitable that the chip-seq signal of these types of studies would span a broad spectrum and that the weak signal scenario, as in this study, would be likely to occur often.
the need for methods to address this problem is acute considering the increasing number of studies using chip-seq to study nr and their co-regulators due to their importance in normal development and in many diseases such as breast cancers.
our work strives to explore whether the peak calling can be improved through the integration of available diverse biological sources via machine learning approaches.
our results demonstrate that it is informative to generate, collect, and integrate the following information: chip-seq data reflecting location of the primary interaction of the tf of interest to its cis-regulatory sites, gene expression data reflecting functional outcomes of interaction of the tf and its co-regulators, and finally genomic sequence data of the identified regions.
other types of data which is highly likely to be useful include histone modification marks, recruitment of rna polymerase ii, and relative location of the insulator protein ctcf  <cit> .
in summary, our results indicate that a supervised classification approach enables one to utilize even limited amounts of existing knowledge together with multiple types of biological data to enhance the sensitivity and specificity of identifying dna binding sites for co-regulators proteins.
our feature selection experiments indicate that experimental inputs complementary to chip-seq are critical in identifying biological significant signals from chip studies with weak signals due to indirect dna binding.

 transcription factors serve as the final molecules in signal transduction pathways that coordinate expression of target genes. when activated in response to upstream signals, often encoded as chemical ligands and protein modification, tfs bind to their cis-regulatory sites to exert their regulatory effects on their target genes. during the process, tfs often interact with other proteins, which further modulate the function and efficacy of tfs to achieve fine-tuned regulation of gene expression; studying such interactions and regulations is an increasingly important component of studying gene expression systems. nuclear receptors, such as estrogen receptor α, are transcription factors that migrate to the nucleus to regulate downstream target genes. nrs play important biological roles in normal physiology and disease. in particular erα plays an important role in both breast cancer and osteoporosis. upon ligand binding, erα and other nrs are bound by proteins called co-regulators that recruit transcriptional machinery and chromatin modifying enzymes. co-regulators are therefore critical in nr activity. understanding the composition of functional nr/co-regulator complexes in specific signaling contexts could provide a basis for the development of novel nr- and co-regulator-targeted therapeutics. the problem addressed in this paper arose from a study of the interaction between the major erα co-activator src-, also known as ncoa with erα and the impact of such interactions gene expression. recently, chromatin immunoprecipitation coupled with high-throughput next-generation sequencing has become the main technology for global characterization of the transcriptional impact of nrs and their co-regulators. chip-seq involves the short-read sequencing of the chip-enriched dna fragments. these short sequence reads are then aligned to a reference genome. then the actual binding loci from the positional tag distributions are determined using'peak calling' algorithms. numerous peak calling algorithms have recently been developed for identifying chip-enriched genomic regions from chip-seq experiments but there is a wide range of discordance among the peak calls from different algorithms. therefore, there is a need for the methods that can integrate additional information besides chip-seq tags to identify functional tf binding sites. furthermore, studying the interactions between tfs and their co-regulators through chip-seq technology poses an additional challenge since co-regulators do not directly bind dna. co-regulator chip-seq measures the secondary protein-dna binding through primary tfs and leads to relatively weak sequencing signals--i.e. relatively small number of sequence tags above noise. as such, it remains a challenge for contemporary peak calling methods to detect weak secondary protein-dna-binding signals and simultaneously maintain a high specificity. often, a well-designed experiment studying interaction between a tf and its co-regulator generates critical information in addition to the chip-seq data for the co-regulator binding. for example, chip-seq data reflecting the binding of the primary tf of interest to its cis-regulatory sites are often collected; the genomic sequence surrounding chip-seq peaks are usually available, which can be used to reflect the intrinsic sequence characteristics of regulatory sites; transcriptomic data that reflect functional outcomes of the interaction of the tf and its co-regulators can also be monitored. in this study, we investigated and compared different statistical and machine learning approaches to integrate multiple types of information to overcome the difficulty of identifying functional erα/src- interaction in presence of weak chip-seq signal.  the biological study underlying this paper aims to investigate the impact of erα/src- interaction on estrogen induced gene expression in a bone cell line transfected with erα, which may shed light on the effect of estrogen-related bone development, bone loss, and potentially bone metastasis. we have generated chip-seq data using anti-erα and anti-src- antibodies in presence and absence of estradiol. to further investigate the impact of interactions between this nr//co-regulator pair, we collected expression array data from the same cell lines with a combination of e treatment and src- knock down. in general, the results of an src- chip-seq experiment would reflect secondary, indirect binding of src- to dna through multiple nrs. however, in this study our experimental design aims to investigate specifically estrogen-induced interactions between erα, src- and dna. the detailed results of the experiments are being prepared for a separate publication. in the current paper, we address the fundamental issue of identifying reliable and functional erα/src- dna binding sites. identifying src- binding sites based on anti-src- chip-seq we first set out to investigate the efficacy of studying erα/src- dna binding sites only based on the results of chip-seq experiments performed using an anti-src- antibody. potential src- binding peaks were identified using three different algorithms: macs., bayespeak, and t-pic. table shows the number of peaks identified by the above algorithms with different cut-off thresholds and the corresponding number of genes to which the peaks are mapped. * union and intersection of the peaks by the three methods, as shown in figure the results of the table raise the following issues during interpretation: first, as expected, applying different cut-off thresholds to the results by a given algorithm leads to a different number of peaks being identified: there is inevitably a trade-off between the number of peaks recovered and the quality of the peaks. second, different algorithms assess the quality of the peaks based on different assumptions and methodologies: there is no consensus on the "goodness" of quality scores of these algorithms. we further noted that different versions of a same algorithm return different quality scores. finally, different algorithms return disjointed sets of peaks, as shown in figure indicating that distinct assumptions and approaches enable an algorithm to discover some potential peaks that evade detection by other algorithms. these issues force decisions that potentially impact the conclusions of the study such as: which algorithm performs better, what cut-off threshold for a given algorithm to pick, and how to consolidate the results from different algorithms so that one can maximize the number of high quality peaks. making these choices remains challenging due to the lack of consensus in the field. in order to compare our results with a recently published study by lanz et al., which analyzed dna recruitment of the co-regulators src- we studied the estrogen-induced src- peaks identified from our data using macs algorithm with a cut-off threshold of p value at 10- a threshold based on their study. our analysis yielded a total of peaks, which were further mapped to genes. the number of peaks identified by us with the above condition is far fewer compared to their study. the discrepancy is likely to be, at least in part, due to amplification of src- in mcf- cells used in their study and possible differences in antibody affinities. however, these results also raised the hypothesis that chip-seq signals of secondary binding at physiologic levels are usually weaker. this suggests that conventional cutoff thresholds for peak calling algorithms may be too stringent, neglecting weak peaks potentially resulting from real erα/src-1/dna interactions. therefore, additional information besides src- chip-seq tags should be capitalized to enhance identification of functional binding sites. integrating multiple sources of biological information for identifying src- binding sites to corroborate the results of src- chip-seq, we also studied the erα chip-seq data and investigated peaks overlapping between the erα and the src- chip-seq results. by varying the cut-off threshold of macs, we identified different numbers of overlapping peaks between erα and src- with the number of overlapping peaks increasing as the cut-off threshold relaxes. the results again indicated that the conventional cut-off thresholds are failing to identify putative erα/src- dna binding sites. on the other hand, simply relaxing the cut-off threshold is likely leading to increased false positive peak calls. thus a principled method is needed to further identify functional erα/src- dna binding sites. to elucidate functional erα/src- dna binding events, i.e., the binding events that influence gene expression, we generated and analyzed expression array data from cells that were treated with vehicle and e in the presence and absence of src- sirna. the microarray data enabled us to identify genes whose expression response to estrogen treatment required src- which are hereto referred to as src-1-sensitive genes. when we compared the list of src-1-sensitive genes and the list of genes with src- binding sites derived by macs at p value = 10- we noted that only genes overlapped among the lists. while the discrepancy between the number of src-1-sensitive genes and genes with src- peaks could be explained by other biological factors, such as that many src- peaks were not functional or secondary expression effects, it also supported our general hypothesis that erα/src- interaction chip-seq signal is relatively weak and potentially true functional erα/src- dna binding sites were missed by the stringent setting of the peak calling algorithm. while it may be tempting to directly combine the information from src- chip-seq, erα chip-seq and microarray data by identifying the intersections of overlapping genes and peaks, such an approach is overly simplistic and ignores other potentially informative data, e.g., the genome-sequence characteristics of erα/src- interaction sites and the prior information of known erα/src- interactions. these considerations motivated us to investigate and compare different principled machine learning approaches in order to improve the sensitivity and specificity of detecting erα/src- dna binding sites by integrating multiple types of information. an integrative approach to detect erα/src- dna binding sites the overall framework and rationale of our information integration approaches are as follows. we formulated the task of identifying functional erα/src- dna binding sites as a classification task, in which we performed a binary classification to label a potential src- binding site derived from chip-seq analysis as either functional or nonfunctional. we investigated both supervised learning, which allows us to take advantage of existing knowledge of erα/src- interactions, and unsupervised learning, which allowed us to take an unbiased approach. the classification formulation allowed us to pool more candidate peaks identified by different peak calling algorithms at relaxed cutoff thresholds so that we did not have to rely on a single "best" algorithms and "optimal" parameterization but resorted to our classification to identify functional erα/src- dna binding. in this study, we collected the union of the peaks returned by all three algorithms at the cutoff threshold as follows, macs: p value cutoff 10- bayespeak: posterior probability ≥ and t-pic: p value cutoff 10- this led to a pool of candidate peaks. another important advantage of the classification approach is that it allows us to integrate multiple types of biological information collected from our experiments and public databases by representing them as features for a classifier. in this way, multiple types of information contribute to the classification of potential peaks and their impact can be determined by learning algorithms. for each candidate src- peak, we constructed the following features: a vector of binary features representing presence/absence of nucleotide trigrams, which reflects the intrinsic characteristics of genome sequence surrounding the summit of a peak region; an average of predicted nucleosome-occupancy scores, which represents the chromatin structure characteristics around the peak summit; a binary feature reflecting if a primary binding peak, i.e., the erα chip-seq peak, overlaps with the src- peak; and a binary feature representing the functional outcome of erα/src- interaction, i.e., whether the peak is mapped to an src- sensitive gene. detailed descriptions of features are presented in methods section. we evaluated the results of predictions from classification algorithms by determining if conserved erα binding motif can be found in the classified peaks, as an indication that a peak is the result of erα/src- dna binding. searching for instances of conserved tf binding motifs at the predicted binding loci is considered the most prominent verification method for validating peaks. unsupervised classification first, we explored if the candidate peaks could be divided into two distinct groups by unsupervised learning in an unbiased manner. we applied a k-means clustering procedure to the data and the results are listed in table we inspected the genome sequences of the peaks to assess if a conserved motif for estrogen response element was detected by a motif classification algorithm referred to as clover. in cluster 46% of peaks contained the ere motif, and, in cluster 27% of peaks contained the ere motif. we believe this is not a good separation of the peaks in that, even though cluster has more ere-containing peaks, it is a bigger cluster and only 46% of peaks contain eres. thus the results would likely lead to a high false positive rate with respect to src- binding. supervised classification supervised learning requires labeled data as training cases. obtaining a training set through large-scale experimental validation of erα/src- dna binding is costly and difficult to perform. therefore, we investigated whether a relatively small amount of labeled data in the supervised learning task would result in better separation compared to unsupervised clustering with our feature set. we experimentally validated src- peaks as functional peaks by quantitative pcr experiment. we used these peaks as positive training cases, together with a set of randomly drawn control chip-seq peaks as negative training cases, to train supervised classifiers. we investigated the performance of three state-of-the-art classifiers: naive bayes, support vector machines and random forest. for nb and random forest classifiers, we set the classification thresholds at and respectively. since the ratio of the positive and negative training cases may have impact on classification algorithms, e.g. nb and svm, we explored using different ratios for training, between 1: 1: and 1: and classifiers were built from these training sets. our test set consisted of all candidate peaks. table lists the total number of peaks in each class, the number ere-containing peaks in each class and the ratio reflecting ere enrichment. results for classifiers with 1: 1: and 1: training ratio were very similar to each other. therefore, just results for the 1: ratio were shown. we noted that the supervised classification approaches have significantly increased the number of positive peaks when compared to those derived by peak calling algorithms based on recommended cutoff thresholds. for example, nb returned positive peaks in comparison to and peaks returned by macs with cutoff p value set at 1e- and 1e- which reflected a 6-fold and -fold increase, respectively. through further evaluation enrichment of ere in the genome sequences surrounding the peaks, we found that a similar percentage of peaks contained ere element: 69% for nb, and 72% and 66% macs at 1e- and 1e- respectively. thus, the results indicate that the qualities of the positive peaks returned by nb were as good as those returned by the stringent peak calling in terms of ere enrichment. we further inspected if the classification approaches retrieved additional functional peaks, i.e., the peaks that were mapped to src-1-sensitive genes derived from microarray experiment. table shows the results of the src- peaks returned by different peak calling approaches that were mapped to src-1-sensitive genes. we noted that, by setting macs cutoff p values at 1e- 1e- and 1e- a total of and peaks were mapped to src-1-sensitive genes. on the other hand, nb has identified peaks that overlap with src-1-sensitive genes. we also note that bayespeak exclusively identified some of the newly "discovered" functional peaks. similarly macs also exclusively discovered some new peaks. these results indicate that, based on different assumptions and criteria, different peak calling algorithms are capable of identifying potential peaks to complement other peak calling algorithms. thus it is more sensible to consider candidate peaks from more than one peak-calling algorithm as long as an objective approach can be further applied to consolidate the results. we also performed a 9-fold cross-validation experiment to assess if the algorithm can correctly identify the experimentally validated erα/src- training cases from the candidate peak pool. results were listed in table nb classifier showed 89% precision, 100% recall and 96% accuracy, see methods section for the descriptions of the metrics. this result increased our confidence that positive calls from our algorithm are likely to reflect real erα/src- dna binding.the results in the table indicate that the nb classifier performed better than the svm and rf classifiers, judging from relative enrichment of ere containing binding sites in the predicted positive peaks. among the three classifiers tested in this study, the rf classifier performed worst. we noted that the number of features that were used by rf during the learning was much smaller than the number of features utilized by other classifiers, which may partially explain the inferior performance of this algorithm in this experiment. the svm method also performed worse on this task than probabilistic nb. we conjecture that the reason might be that svm is complex algorithm with many parameters to adjust and therefore finding optimal parameters for decision boundary might be challenging for this task. we therefore concentrated on the nb classifier because it could be readily used in both supervised and semi-supervised learning environment. semi-supervised classification our number of training cases is relatively sparse compared to a conventional machine learning setting. semi-supervised approaches have been applied to conquer limitations of supervised and unsupervised methods when labeled data is scarce and obtaining large amounts of labeled data is expensive and time consuming. this is done by incrementally assigning instances, which are called with high confidence by a classifier, from unlabelled data into training cases in order to increase the number of training cases and thus enhance the generalizability of classification. therefore, we investigated semi-supervised classification to see whether we could further increase the performance in identifying erα/src- dna binding. in this study, we applied a self-training algorithm using nb as the base classifier because of its probabilistic outputs. we iteratively assigned the most confident positive instances called by our classifier into training cases and found that performance of the self-training became stable after iterations and stopped further training. figure shows the trends of percentage of ere-containing peaks in the positive calls in semi-supervised learning. it is interesting to note that initially as a few pseudo-positive cases were imputed into the training cases the precision of the called positive peaks decreased but later became stable after iterations. a similar total number of peaks and the percentage of ere-containing peaks were identified by our semi-supervised learning algorithm when compared to other supervised learning experiments, see table thus the results do not show obvious advantage of semi-supervised learning over supervised learning algorithms in our experiment. identification of informative features biologically, it is of interest to identify the features that significantly contribute to the classification in that it will reveal the relationships between input features and outcome. we rank key features by roc class reparability criteria using matlab, using a training dataset containing the true positive peaks and random non-binding sites from igg peak calls. following were the top ranked features: "aac", peaks-mapped-to-src-1-dependent-genes, "gcg", "cgt", overlapping-with-erα-peak, "aag", "aca", "acc", "cgc", "aga", "agc", "agg", "cga", "act", "atc", see methods section for detailed descriptions of the features. among the top-ranking features, we noted that the features reflecting the function outcome and the interaction between er and src- were ranked high, indicating the learning algorithm correctly recognized their importance in classification. it is interesting to note that many nucleotide trigrams, which reflect the characteristics of sequences of peaks, were among the high-ranking features. we aligned the top-ranking trigrams to the ere motif, as shown in figure indeed, the trigrams correspond well with the important components of the ere motif. these results indicate that the classification learning algorithm, like the motif searching algorithms, is able to identify highly conserved "words" that constitute one of the important motifs of the training sequences. we noted that the feature reflecting nucleosome occupancy at peak regions was ranked as 48th. this may indicate either that nucleosome occupancy is a dynamic process and our static feature does not reflect the true occupancy status during the experiments or that the era/src- dna binding is not heavily dependent on nucleosome occupancy. biological insights from improved peak calling we further examined the impact of the improved src- peak calling on biological insights drawn from the dataset. we conducted gene ontology analysis using the database for annotation, visualization, and integrated discovery on genes with an src- peak within kb of the tss as determined by macs or by our method. we observed a dramatic difference in the identification of genes enriched in specific biological processes. specifically, our method yielded in the calling of peaks in gene sets which were highly enriched for genes involved in blood vessel development and actin filament-based processes. indeed blood vessel development has previous been implicated in bone generation. further, within the genes enriched in these biological processes, we identified a number of genes with known functions in bone development. since src- has already been implicated in e mediate bone maintenance, this observation provides evidence for the mechanism underlying this phenotype. in contrast, genes with src- peaks determined by macs were not significantly associated with any biological processes.  we believe that the ability to improve chip-seq peak calling by utilizing available sources of biological information for indirect co-regulator binding in the presence of weak chip-seq signal is an important research area. due to the intrinsic variability in the affinity of interactions between a tf and its co-regulators, it is inevitable that the chip-seq signal of these types of studies would span a broad spectrum and that the weak signal scenario, as in this study, would be likely to occur often. the need for methods to address this problem is acute considering the increasing number of studies using chip-seq to study nr and their co-regulators due to their importance in normal development and in many diseases such as breast cancers. our work strives to explore whether the peak calling can be improved through the integration of available diverse biological sources via machine learning approaches. our results demonstrate that it is informative to generate, collect, and integrate the following information: chip-seq data reflecting location of the primary interaction of the tf of interest to its cis-regulatory sites, gene expression data reflecting functional outcomes of interaction of the tf and its co-regulators, and finally genomic sequence data of the identified regions. other types of data which is highly likely to be useful include histone modification marks, recruitment of rna polymerase ii, and relative location of the insulator protein ctcf. in summary, our results indicate that a supervised classification approach enables one to utilize even limited amounts of existing knowledge together with multiple types of biological data to enhance the sensitivity and specificity of identifying dna binding sites for co-regulators proteins. our feature selection experiments indicate that experimental inputs complementary to chip-seq are critical in identifying biological significant signals from chip studies with weak signals due to indirect dna binding. 
however, due to the nature of detecting indirect protein-dna interaction, chip-seq signals from co-regulators can be relatively weak and thus biologically meaningful interactions remain difficult to identify. our results indicate that supervised learning with naïve bayes algorithm significantly enhances peak calling of weak chip-seq signals and outperforms other machine learning algorithms. there is an increasing interest in understanding the mechanism of action of co-regulator proteins, which do not bind dna directly, but exert their effects by binding to transcription factors such as the estrogen receptor. our results indicate that a supervised classification approach enables one to utilize limited amounts of prior knowledge together with multiple types of biological data to enhance the sensitivity and specificity of the identification of dna binding sites from co-regulator proteins. our integrative approach revealed many potential erα/src- dna binding sites that would otherwise be missed by conventional peak calling algorithms with default settings. in this study, we investigated and compared different statistical and machine learning approaches including unsupervised, supervised, and semi-supervised classification approaches to integrate multiple types of genomic and transcriptomic information derived from our experiments and public database to overcome difficulty of identifying functional dna binding sites of the co-regulator src- in the context of estrogen response. chromatin immunoprecipitation coupled with massively parallel sequencing is increasingly being applied to study genome-wide binding sites of transcription factors.
 chromatin immunoprecipitation coupled with massively parallel sequencing is increasingly being applied to study genome-wide binding sites of transcription factors. there is an increasing interest in understanding the mechanism of action of co-regulator proteins, which do not bind dna directly, but exert their effects by binding to transcription factors such as the estrogen receptor. however, due to the nature of detecting indirect protein-dna interaction, chip-seq signals from co-regulators can be relatively weak and thus biologically meaningful interactions remain difficult to identify.  in this study, we investigated and compared different statistical and machine learning approaches including unsupervised, supervised, and semi-supervised classification approaches to integrate multiple types of genomic and transcriptomic information derived from our experiments and public database to overcome difficulty of identifying functional dna binding sites of the co-regulator src- in the context of estrogen response. our results indicate that supervised learning with naïve bayes algorithm significantly enhances peak calling of weak chip-seq signals and outperforms other machine learning algorithms. our integrative approach revealed many potential erα/src- dna binding sites that would otherwise be missed by conventional peak calling algorithms with default settings.  our results indicate that a supervised classification approach enables one to utilize limited amounts of prior knowledge together with multiple types of biological data to enhance the sensitivity and specificity of the identification of dna binding sites from co-regulator proteins. 17- january the tenth asia pacific bioinformatics conference melbourne, australia 