Efficient yeast ChIP-Seq using multiplex short-read DNA sequencing

Similar barcode behavior with the same DNA sample

To first test whether barcodes generate similar results when added to the same population of DNA fragments, we divided an input DNA sample into five aliquots. The input sample consists of Saccharomyces cerevisiae DNA prepared from crosslinked and sonicated chromatin without immunoprecipation and it represents breaks at open chromatin regions (R.K. Auerbach, G.M. Euskirchen, unpublished). Four of the samples were differentially barcoded and the libraries were mixed together in the same flowcell lane. As a control, standard non-barcoded Illumina genomic DNA adapters were added to the fifth aliquot and this library was sequenced in another lane. After sorting reads by barcodes, the barcode sequences were stripped and 26 base pairs were aligned back to the yeast genome. Similar numbers of reads were obtained for each barcoded library. As shown in Figure 2, signal tracks are very similar for the five aliquots. The normalized average tag counts in 500 bp windows are also highly comparable among aliquots (Table 2, R² = 0.93–0.97). The results for each barcoded library were scored against the control without a barcode to find significant differences; only 36,992 out of 12,156,677 unique bps (0.304%) from the S. cerevisiae genome differ significantly for all four barcoded libraries (Figure 2; See Methods for analysis). In addition, comparison of peaks among the different libraries revealed high degrees of similarity. These characteristics highly resemble those of a non-barcoded input library analyzed in two different sequencing reactions. Thus, overall biases are not observed in efficiencies of various barcoded libraries and barcoding does not appear to induce differences relative to non-barcoded libraries.

Experimental design for multiplex yeast ChIP-Seq

We next analyzed the distributions of three diverse DNA-binding proteins across the yeast S. cerevisiae genome using three to four biological replicates per factor. This study allowed us to 1) further examine the reproducibility of a specified barcode using a different biological replicate 2) determine whether there was crossover signal between barcodes and 3) use the technique to map the binding distributions of three DNA-binding proteins. The proteins selected were the centromeric histone H3 variant Cse4 [16], RNA PolII (each grown in vegetative growth) and the transcription factor Ste12 [31] (incubated in pseudohyphal growth conditions). Cse4 is expected to have a very distinct binding distribution from that of Ste12 which binds near promoters. We also included input DNA as a reference for ChIP-Seq as a fourth sample. Initially we compared three biological replicates for each factor using barcoded sequences for two of the replicates and non-barcoded sequences for the third one. Each biological replicate was made into its own DNA library. The barcoded libraries were pooled and sequenced with one pool per lane (Table 3) whereas the non-barcoded libraries were sequenced in individual lanes. Two lanes were run using the barcoded adapters: Pol II-ACGT, Input-CATT, Cse4-GTAT and Ste12-TGCT (Pool #1, Table 3) and Input-ACGT, Cse4-CATT, Ste12-GTAT and PolII-TGCT (Pool #2, Table 3). To compare reproducibility of alternative barcodes relative to reproducibility between biological replicates, we sequenced another replicate of Pol II ChIP DNA in a separate pool using the same barcode as one of the other replicates (Pol II-ACGT) (Pool #3, Table 3). After sorting reads by barcodes, the barcodes were stripped and 26 base pairs were aligned to the yeast genome. We scored ChIP samples against 8 million input reads derived equally from barcoded and non-barcoded input DNA. Rank-order comparisons between the top targets for two different biological replicates were performed to assess similarity. We obtained 148 targets for Cse4, 823 targets for Ste12 and 2508 targets for PolII (FDR = 0.05; Additional Files 1, 2, 3).

Comparison of PolII biological replicates

By comparing two biological replicates of PolII ChIP DNA with the same barcode (ACGT) we found a close correlation between the target lists (Figure 3a and 3c). A similar close correlation was observed by comparing the targets of the individual biological replicates using different barcodes among each of the different factors (Figure 3c). Furthermore, comparison of normalized average tag counts in 500 bp windows between biological replicates for PolII and Ste12 revealed strong correlations (Table 2, R² = 0.89–0.97). These values are similar to the correlation coefficients between biological replicates found in a recent yeast RNA-Seq study (R² = 0.87–0.90) [1]. Thus, consistent with barcoded input DNA, the different barcodes give similar results when analyzing similar libraries. Our results strongly suggest that the barcode has negligible impact on binding site detection.

Comparison of ChIP-Seq data to previously published ChIP-chip data

We also compared the binding site distributions for Ste12, Cse4 and PolII to previously published data [16,36,37,54,56]. Our Ste12 ChIP-Seq targets overlap with 67% (232/346) of the targets found in Borneman et al. [37] while our PolII ChIP-Seq targets overlap with 72% (667/929) of the PolII ChIP-chip targets from the Steinmetz et al. data [54]. These overlaps between ChIP-Seq and ChIP-chip do not differ significantly from those determined in other studies [6], typically in the 64–71% range. However, we observe 2- to 4-fold more targets using ChIP-Seq as compared to ChIP-chip, consistent with the observation that ChIP-Seq is more sensitive than ChIP-chip [6].

Non-overlapping binding profiles among samples sequenced simultaneously

By comparing binding profiles obtained from a single lane containing four different barcoded libraries, we observed that many targets belong uniquely to a single signal track and their distributions are consistent with the biological roles of the factors studied (Figure 3b and 3c). This is particularly the case for Ste12 and Cse4 which have highly dissimilar binding patterns and almost no target overlap (Figure 3c). Note that some peaks are present in both ChIP and input samples and are removed after scoring. Thus, analysis of different ChIP samples with distinct barcodes indicates that the signal from one barcoded library does not crossover to those of other barcoded libraries (Figure 3b). We define signal crossover as the erroneous assignment of a barcoded read from a particular ChIP sample to the wrong ChIP sample during multiplex ChIP-Seq, generating a signal profile for one barcoded sample with reads from other ChIPs sequenced simultaneously.

Ste12, PolII and Cse4 binding distributions

Analysis of the binding targets of Ste12, Pol II and Cse2 reveals expected and novel results. Ste12p binds upstream of genes implicated in pseudohyphal growth such as Tec1, Mga1, Phd1 and Flo8 (Figure 4 and [31,36]). As expected, Ste12 targets are enriched for GO categories including pseudohyphal growth, cell wall organization, and stress response (Table 4 and Additional File 5), which is consistent with the role of Ste12p in regulating the pseudohyphal pathway. The PolII distribution in mid-log phase intersects ~2500 genes, which is fewer than that reported by recent RNA-Seq studies [1]. This may be due to differences in the analysis of ChIP-Seq and RNA-Seq data. Enriched regions in ChIP-Seq data are determined using a threshold for a minimal read count difference at a particular region between the ChIP sample and input DNA while RNA-Seq data is analyzed using a global threshold for a minimal number of mapped reads. Overall, we recapitulated all the binding profiles displayed in the genome-wide PolII ChIP-chip study by Steinmetz et al. (Figure 5) [54]. As expected, PolII targets are overrepresented in GO categories such as translation, transcription, RNA processing and primary metabolism, which are crucial for exponential growth in rich media (Table 4 and Additional File 5). Cse4p binds strongly to all yeast 16 centromeres as shown previously (Figure 6 and[16]). Our data indicate that Cse4 tightly occupies a narrow region around the centromere. In addition, we also detected significant binding relative to background at 132 novel non-centromeric sites; Cse4p is less abundant at these sites relative to centromeric regions (Additional File 1). These non-centromeric sites are enriched for metabolic genes, in particular those involved in glucose metabolism (e.g. glycolysis, hexose and glucose catabolism), which are highly expressed under our conditions (Table 4 and Additional File 5). In addition, nearly all non-centromeric Cse4 target genes are present among the 100 highest ranked PolII targets. To explain this novel finding, we hypothesize that Cse4 transiently integrates at regions of high histone turnover which could now be detected using ChIP-Seq. Cse4p might temporarily integrate at these non-centromeric sites before being outcompeted by histone H3 since it has been demonstrated by ChIP-qPCR that a Cse4 mutant protein can localize to non-centromeric locations when the stoichiometry between histone H3 and Cse4p is perturbed [27].

Artifacts absent in barcoded libraries but present in non-barcoded libraries

Finally, we did observe one set of differences between barcoded and non-barcoded libraries. A curious rare artifact was found in non-barcoded libraries and excluded in barcoded libraries. ~31 abnormally-shaped peaks were present in non-barcoded ChIP libraries in which the same sequences was obtained many times resulting in block-shaped peaks (Additional File 4). This artifact was not observed in barcoded libraries.

Although the source of this artifact is not clear, it might be due to specific priming from the non-barcoded adaptor on yeast DNA. Regardless, these artifacts are readily eliminated by removing non-unique reads. We did not find additional areas of dissimilarity between barcoded and non-barcoded libraries.

Data accession and supplementary methods

The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [58] and are accessible through GEO Series accession number GSE13322 [59]. Sgr files and other data can be accessed at this URL address: http://archive.gersteinlab.org/proj/Lefrancois_et_al_2008/

Cell growth

Yeast strains are described in Table 8. Yeast strains CMY288-1B, CMY8058-3-4 and CMY8082-20-3 were grown in YPAD rich media to exponential mid-log phase (OD₆₀₀ = 1.0, 500 mL culture). For yeast strains YJM339 and CMY291, cells were grown overnight to mid-log phase to OD₆₀₀ = 0.6 and then pseudohyphal growth was conducted in nitrogen-depleted SLAD for 4 h as described [36].

Chromatin immunoprecipation (ChIP)

Chromatin immunoprecipitations were performed as described [60]. All ChIP experiments were completed as biological triplicates for Cse4 and Ste12 and as biological quadruplicates for PolII. For the immunoprecipitations of Ste12 and Cse4, Ste12-13X Myc (160 uL of a 50% anti-Myc EZiew affinity gel; Sigma) and Cse4-3HA (160 uL of a 50% anti-HA EZiew affinity gel; Sigma) were pre-washed three times in lysis/IP buffer containing protease inhibitors and PMSF and added to the lysates from their respective epitope-tagged strains and untagged controls. For the immunoprecipitation of native RNA polymerase II from strain CMY288-1B, 20 μL of mouse ascites containing RNA polymerase II 8WG16 mouse monoclonal antibody (Covance, Cat. #MMS-126R) was added to one set of triplicates. The other set was left without antibody as a control.

Immunoprecipitations were carried out with inverting at 4°C for 14–16 h. For the PolII samples, 250 uL of a 50% Protein G agarose slurry (washed twice in lysis/IP buffer) was added to each of the PolII samples as well as to the no antibody controls, and incubated for 1 h at 4°C with gentle inverting.

Input DNA preparation

To isolate input DNA, 250 uL of clarified CMY288-1B lysate was reserved before immunoprecipitation and combined with 250 μL of 1× TE [pH 8.0]/1% SDS. Crosslinks were reversed by an overnight incubation at 65°C and Proteinase K treatment was carried out as described above. Instead of direct DNA precipitation as for ChIP samples, input DNA was extracted three times using phenol:chloroform:isoamyl alcohol (25:24:1) (Fluka) and once with chloroform alone. To precipitate DNA, 50 μL of 5 M LiCl and 1 mL of 100% ethanol were mixed with the upper phase isolated from the chloroform extraction and precipitation occurred for 1 h at -20°C. Samples were centrifuged for 20 min at 14,000 RPM at 4 °C and were air-dried for 10 min. DNA was resuspended in 30 μL 1× TE [pH 8.0] and input DNA was incubated at 37°C for 30 min with 2 μL of DNase-free RNase A. Finally, input DNA was further purified using a Qiagen MinElute column.

Real-time quantitative polymerase chain reaction

Real time quantitative PCR (qPCR) with SYBR green dye was performed using a Roche LightCycler 480 qPCR machine to confirm enrichment of control regions. A two-fold enrichment between experimental samples (epitope-tagged strains for Ste12 and Cse4, primary antibody for PolII) and reference samples (untagged strains for Ste12 and Cse4, protein G beads only for PolII) for known binding sites was set as the threshold for samples to continue forward towards Illumina library preparation. For each DNA-binding protein studied, four primer pairs were designed and these included at least two primer pairs amplifying known binding sites and at least one primer pair amplifying a random region as a negative control. Primers were generated using Primer3 http://primer3.sourceforge.net/ with the following criteria: fragment size between 200 and 250 bp, primer length 20 and melting temperature between 59°C and 61°C. Other settings were kept as default. Primer sequences are given in Table 9. A standard curve with a dilution series of genomic DNA was generated for each primer pair to determinate primer pair efficiency [61]. Each 10 μL reaction contained 2 μL of nuclease-free water (Gibco), 5 μL of LightCycler 480 SYBR Green I Master mix (Roche, Cat. #04 707 516 001), 0.5 μL of each primer (10 μM stocks) and 2 μL of either water (negative control), diluted genomic DNA (standard curve) or ChIP DNA (1/21 dilution). For ChIP DNA, reactions were done in duplicates on the 384-well plate. Phenol-chloroform extracted genomic DNA was diluted to 8^-1, 8^-2, 8^-3, 8^-4, 8^-5, 8^-6 and 8^-7 to generate a dilution series for determination of the standard curve. Each qPCR reaction was run using the same program: pre-incubation (95°C for 5 min), 45 cycles of amplification (95°C for 20 s, 54°C for 30 s, 72°C for s), melting curves using a heat ramp and cool down. Crossing point values (Cp values) were obtained using the second derivative maximal analysis tool included in the Roche LightCycler480 software. To determine enrichment, the 2^-ΔΔCp method was used by comparing enrichment values for positive primer pairs to a negative primer pair between experimental and reference ChIP experiments. The negative primer pair should not give differences between experimental samples and control samples. Average enrichment for all biological replicates was determined for each primer pair and standard deviation to the mean was indicated as error bars. Mean enrichments over negative primer pair (+/- SD) were plotted in Excel (Figure 7).

Adapter design and annealing

Each barcode has a 3 base index followed by a 'T' (position 4) to promote pairing with the 'A'-overhang that is added to the samples as part of the Illumina genomic DNA library preparation. Barcodes were created in such a way that no barcode had the same base at positions 1, 2 and 3 for two reasons. First, we wanted to have a balanced base composition. Second, since ELAND alignment software allows 2 mismatches to map any read, we wanted to make sure that sequencing errors would not influence mapping. These four barcodes were used for this study: GTAT, CATT, ACGT and TGCT.

Oligos were synthesized by MWG at a 0.05 μmol scale with HPLC purification. Oligo sequences are given in Table 1. Each oligo was resuspended in annealing buffer (10 mM Tris [pH 7.5], 50 mM NaCl, 1 mM EDTA) to 200 μM. The forward and reverse oligos for each pair were mixed in equal volumes to 100 μM and denatured for 5 min at 95°C in a wet heat block. The heat block was then removed to room temperature and allowed to cool slowly over 45 min to promote annealing. For input DNA, Illumina genomic DNA adapters were diluted 1:20 using RNase-free DNase-free water (Gibco) and our adapters were diluted 1:30. For ChIP DNA libraries, the DNA concentrations of the barcoded-adapters were measured by a Nanodrop spectrophotometer and adjusted to the same working concentration as Illumina's adapters. The following dilutions were therefore determined: 1:500 for barcode GTAT, 1:450 for barcode CATT, 1:750 for barcode ACGT and 1:330 for barcode TGCT.

Illumina DNA library generation

We followed the manufacturer's protocol for creating genomic DNA libraries with a few modifications based on experience gained from ChIP-Sequencing [6]. Here we present our modifications for generating barcoded libraries. ChIP DNA and input DNA were first band-isolated on a 2% agarose to obtain fragments between 150 and 350 base pairs and DNA was extracted using the QIAquick gel extraction kit (Qiagen) and eluted in 34 μL. For input DNA, because DNA amounts are higher than DNA recovered by ChIP, input DNA after gel extraction was diluted 1:5. After end-repair and addition of a single adenosine ("A") nucleotide, adapters were ligated to samples for 15 min at room temperature in the following fashion: the samples eluted from the MinElute column in 10 μL were ligated to 1 μL of adapters using 1.3 μL of LigaFast T4 DNA Ligase (3 Units/μL; Promega) and 12.3 μL Rapid Ligation Buffer (Promega). For input DNA libraries, the Illumina genomic DNA adapters were diluted 1:20 and all barcoded adapters were diluted 1:30. For ChIP DNA libraries, the Illumina genomic DNA adapters were diluted 1:40 and barcoded adapter dilutions are described in the previous section. After 15 min, samples were purified with the MinElute PCR purification kit (Qiagen).

We found that size-selection and purification of the ligation products by agarose gel before amplification of the library by PCR resulted in better libraries and decreased the incidence and intensity of an adapter-adapter band at ~120 bp. To eliminate those adapters that lack a fragment insert, samples were run on a 2% agarose E-Gel (Invitrogen) for 20 min, together with Track-It 50 bp DNA ladder (Invitrogen). DNA fragments ranging from 150 base pairs to 500 base pairs were extracted and recovered in 28 μL EB with a QIAquick gel extraction kit (Qiagen). To amplify the library, PCR was performed using Illumina genomic DNA primer "1.1" and Illumina genomic DNA primer "2.1" as described [6,62] with 15 cycles (Input DNA) or 17 cycles (ChIP DNA) of amplification. A final size selection was performed using a 2% agarose E-Gel to obtain a library with a median length of ~230 bp which is within the recommended size range for cluster generation on Illumina's flowcell. The library was recovered in 20 μL EB using MinElute Gel Extraction kit (Qiagen). Finally, DNA concentrations and purities (A_260/280 nm ratios) were measured on a Nanodrop spectrophotometer and are given in Table 10.

Mapping values and scoring

Raw data from the Illumina Genome Analyzer I and II were analyzed with Illumina's Firecrest, Bustard and GERALD modules for image analysis, basecalling and run metrics respectively, and a PhiX174 control lane was used for matrix and phasing estimations, as per the manufacturer's instructions. At this stage, a Perl script was used to partition the reads and remove the barcodes, i.e. the first four bases of each read. The next 26 bases of each read were aligned against the reference genome S288c Saccharomyces cerevisiae, using Illumina's ELAND program in standalone mode. For each barcode from each flowcell lane of barcoded libraries, the numbers of total and mapped reads were determined. Reads lacking a fully-intact barcode were discarded in a fifth bin called unclassified (not shown) and were not used in individual barcode mapping analysis, although they are calculated in the global lane mapping analysis. Mapping values and thresholds for scoring are given in Table 10. Uniquely-mapped reads with ChIP factors and barcoding schemes in common were pooled to produce 19 different data sets (Table 10). A full lane of non-barcoded input DNA comprising 2,455,181 mapped reads was used as a control when scoring barcoded input DNA. For scoring ChIP DNA relative to input DNA, a reference set consisting of two full lanes of input DNA and two barcoded input DNA data sets were combined, adding up to 13,198,172 total reads. Signal files were created using a 200 bp sliding window and scoring was performed with the PeakSeq program [63] using a mappability fraction of 1.0. Significant "ChIP hits" were determined relative to the corresponding input DNA by PeakSeq and further filtered by requiring a hit length of at least 100 bp, a p-value of < 0.05, a ratio of at least 2.0 between ChIP DNA and input DNA read counts and a difference of at least 10 between the ChIP DNA and input DNA read counts.

In the section "Similar barcode behavior with the same DNA sample", each barcoded input DNA (Table 10, experiments #15–18) was scored against the non-barcoded input DNA (Table 10, experiment #14) using PeakSeq. Similar criteria were applied to filter barcoded input DNA hits. In these significant "input hits" from all four barcoded input DNA samples, the number of overlapping nucleotides was determined. Finally, the total number of nucleotide positions from all significant "input hits" was calculated and compared to the S. cerevisiae genome size.

Signal tracks and browser

Signal tracks were produced from uniquely-mapped reads and the S. cerevisiae genome using a 200 base pair sliding window. The Integrated Genome Browser (IGB, Affymetrix) was used to view images of signal tracks and to overlay them against the October 2003 version of the Saccharomyces cerevisiae genome. For ORFs and other annotations, the Saccharomyces cerevisiae genome was imported from SGD, the Saccharomyces Genome Database http://www.yeastgenome.org. After barcode parsing, each barcoded sample was compared over input DNA signal. To build signal tracks for comparing samples with different number of sequencing reads, the y-axis was normalized for each sample according to the total number of mapped reads.

Comparisons of the average tag counts in 500 bp bins between two experiments

We divided the yeast genome into 500 bp bins. For each experiment (Table 10), we calculated the average tag count for each bin by summing the tag counts at each nucleotide position within a bin and dividing the total by 500. We then normalized each experiment by dividing the average tag count for each bin by the total tag count for the experiment. Normalized average tag count for each 500 bp bin were compared between two experiments among barcoded and non-barcoded samples for a particular factor. Linear regression was performed using R and the R² (R-squared) and p-values recorded.

Target list annotation

For RNA PolII, target lists were overlapped with SGD ORFs and then annotated using SQL operation chains. Targets that did not overlap with any ORF were manually annotated to the closest gene by comparing PolII bed files in IGB. If we could not discriminate between two genes, then both were selected.

For Cse4, target lists were manually curated. CENs and ORFs were identified manually by looking at Cse4 signal tracks and Cse4 bed files in IGB and associating them to overlapping features. If one target overlapped two ORFs, both ORFs were included in subsequent analyses. If one target was located between two ORFs, the closest ORF was selected if the target's signal was noticeably higher towards this ORF and if the target was positioned closer to this ORF by at least 200 bp. Otherwise, both ORFs were considered.

For Ste12, target lists were manually curated. Targets located within 1.5 kb upstream of an ORF's transcription start site were assigned uniquely to that ORF. In the case of two ORFs with divergent promoters, both ORFs were considered if they were located less than 1 kb away from the target. For all other targets, the closest ORF was selected except when distances between a target and its two closest ORFs were less than 200 bp. In this case, both ORFs were included. Targets whose signals were marginally above background and that were distal to annotated regions of the genome were discarded.

Comparisons between ChIP-Seq and published ChIP-chip data

Comparisons between PolII ChIP-Seq signal tracks and PolII ChIP-chip signal profiles around genes were made using data from Steinmetz et al. [54] and were performed visually in IGB. Comparisons between Ste12 ChIP-Seq targets and Ste12 ChIP-chip targets were done manually using IGB. If a ChIP-chip target was located less than 150 bp away from a ChIP-Seq target in at least one of the biological replicates, they were considered as overlapping for analysis purposes. ChIP-chip targets located more than 2 kb from a gene were excluded from comparisons. For Ste12, the Ste12 ChIP-chip target list from Borneman et al. [37] was used. For PolII ChIP-chip data from Steinmetz et al. [54], readily available on SGD and GEO (Series GSE6293, Sample GSM144667) [64], similar analyses and criteria as in Borneman et al. [37] were used to determine PolII ChIP-chip targets for consistency. 929 PolII ChIP-chip targets were visually inspected in IGB to ensure that they matched the PolII ChIP-chip signal profile [54]. Comparison with PolII ChIP-Seq data was done manually using IGB with similar criteria as those described for the Ste12 comparison between ChIP-Seq and ChIP-chip targets.

Rank-rank correlations

In order to demonstrate that the rank-order lists of target binding sites for different samples are correlated, we performed the following comparisons. Lists of target binding sites are ranked first by p-value from most significant to least significant and in second by the difference between the ChIP DNA and input DNA read counts from the highest to the lowest. For a pairwise comparison of two target lists, we can vary the fraction of each list that is compared. The number of targets that agree (i.e. overlap by at least one bp) between the two lists can be plotted as a function of the fraction of the two lists compared. Ideally for targets that are identical the rank-rank plot will be linear having a value of 100% agreement when the entire rank lists are compared. Deviation for this linear correlation shows where the rank lists being compared begin to differ.

Gene Ontology (GO) analysis

GO analyses were performed on SGD http://db.yeastgenome.org/cgi-bin/GO/goTermFinder.pl. We chose the GO Process Ontology with a p-value cutoff at p < 0.01. Other settings were set at default. For all three factors, ORFs were required to be present in at least one biological replicate. For Cse4 and PolII, all ORFs were used for the analyses. For Ste12, only those ORFs that were downstream of a target were analyzed.

Simulation of depth of sequencing

To investigate the sequencing depth necessary to saturate the number of detectable binding sites, we performed the following computational simulation. We simulated a genome of size S_G, with N_sites binding sites of average fold enrichment f, using N_reads mapped sequence reads. In order to perform the simulation efficiently, we scaled the genome size to 1 Mb, and correspondingly scaled the number of binding sites and mapped reads. We randomly generated sequence reads using a uniform distribution over the 1 Mb interval for both a simulated input DNA sample as well as for a ChIP DNA sample with regions with enriched numbers of sequence reads corresponding to binding sites. These simulated datasets were scored in the same manner as ChIP DNA data was scored against input using a false discovery rate of 5%. The sensitivity (i.e. the fraction of simulated binding sites identified) as well as the positive predictive value (i.e. the accuracy of the target sites identified) of the scored results was computed using the simulated data. The simulation was repeated 50 times in order to accurately estimate the sensitivity and positive predictive values. By varying the number of simulated mapped sequence-reads N_reads, we can compute the minimum number of reads required in order to achieve a sensitivity of 95% (i.e. to be able to identify 95% of the simulated binding sites).