Main

Understanding the human genome will require comprehensive delineation of functional elements within the 98% of genomic terrain that does not encode protein. In vivo, cis-regulatory modalities colocalize with focal alterations in chromatin structure1,2,3,4, and this governs the accessibility of genomic sequences to critical regulatory factors. Exploitation of the close connection between functional elements and chromatin structure should offer a powerful and generic approach for de novo identification of cis-regulatory sequences in the context of complex gene domains.

Active regulatory elements within complex genomes are distinguished by pronounced sensitivity to the nonspecific endonuclease DNase I3,4,5 when exposed in the context of intact nuclei. DNase I hypersensitive sites are the sine qua non of a diverse spectrum of classical transcriptional and chromosomal regulatory activities including enhancers, promoters, silencers, insulators, boundary elements and locus control regions1,3,6. Indeed, in the human genome, many functional elements were first identified as major DNase I hypersensitive sites and only later were found to have specific regulatory roles. Analysis of chromatin structure may enable generic delineation of functional elements across the genome, provided it exhibits direct sequence specificity, quantitative data output that permits automated analysis, and adaptability to a high-throughput format.

DNase I hypersensitive sites in native genomic domains have traditionally been localized by an approach relying on Southern transfer followed by indirect end-labeling5. Although widely applied, this technique is not quantitative and has numerous technical and resource-related limitations that prevent its application on a genome-wide scale. The major limitations of conventional hypersensitivity assays are the low throughput and the lack of sequence specificity. Conventional Southern blot–based hypersensitivity assays are tremendously time-consuming, and comprehensive mapping of DNase I hypersensitive sites over an extended gene locus (50–200 kilobase, kb) requires months, or in some cases years, of dedicated effort. Because it is an indirect assay, another critical drawback of this conventional approach is its lack of sequence specificity. Identification of the core 50–200 base pair (bp) elements over which cooperative trans factor binding occurs often requires extensive additional functional analysis, such as deletion studies7.

We therefore sought to develop a high-throughput, sequence-specific approach to mapping hypersensitive sites. We showed that the continuous distribution of DNase I sensitivity across a gene locus can be accurately sampled in a high-throughput format. These data were analyzed using a rigorous algorithm to construct a quantitative 'chromatin profile' and to identify and score hypersensitive outliers via a signal-to-noise ratio (SNR). Application of this quantitative chromatin profiling (QCP) approach to 300 kb of diverse human genomic terrain demonstrated that peaks in hypersensitivity SNR are highly sensitive and specific for DNase I hypersensitive site core elements and the associated cis-regulatory sequences. This approach may also reveal new elements even in the context of well-explored genomic loci.

Results

QCP

Sensitivity of chromatin to DNase I in vivo is expected to vary as a continuous function of genomic position. We reasoned that if this variation could be captured quantitatively, it would be possible to distinguish DNase I hypersensitive sites by performing outlier detection using rigorous statistical methods. We hypothesized that continuous quantitative profiles of DNase I sensitivity could be obtained by (i) tiling PCR amplicons across a candidate gene locus, (ii) quantifying the relative DNase I sensitivity of genomic DNA corresponding to each amplicon, and (iii) plotting replicate determinations as a function of genomic position (Fig. 1). We further expected that such profiles could be analyzed to reveal a DNase I sensitivity 'baseline' (together with 95% confidence bounds), relative to which statistically significant outliers could be reliably detected using a generic algorithm. DNase I hypersensitive sites typically comprise a core element of roughly 250 bp, where regulatory factors bind and exclude a canonical nucleosome. However, the chromatin structural alterations produced by such complexes vary considerably and may extend to flanking sequences7. We therefore predicted that core DNase I hypersensitive site elements would be situated at the local maxima of hypersensitivity, which could be identified rigorously by computing its SNR as a function of genomic position (Fig. 1).

Figure 1: Overview of quantitative chromatin profiling.
figure 1

To study a candidate gene of interest in the context of a given cell type, intact nuclei are isolated, and DNase I–treated and untreated samples are prepared. Primers are designed to amplify tiled contiguous 250 bp amplicons spanning the candidate gene locus. The relative number of intact copies of the genomic DNA sequence corresponding to each amplicon is then determined (using a suitable modality, such as real-time PCR) in DNase I–treated and untreated samples in replicate (6–9 × ). The resulting DNase I sensitivity ratios (relative copies in treated / relative copies in untreated samples) are then plotted on a genomic axis. Mean values for 6–9 × replicates from each amplicon are computed and define a moving DNase I sensitivity 'baseline,' relative to which 95% confidence bounds are determined empirically. A robust statistical algorithm is then applied to identify outliers, which signify hypersensitive amplicons. A hypersensitivity SNR is then computed for each outlier amplicon and plotted as a function of genomic position. For a given outlier feature, the SNR provides a quantification of the number of standard deviations by which the feature exceeds the expected background variability, given its empirical distribution. Note that SNR plots contain only values for significant outliers. In a final step, local maxima in SNR are determined, defining the peaks in hypersensitivity that overlie core functional elements.

Under standard buffer conditions, DNase I introduces both single-stranded (ss) nicks and double-stranded (ds) cuts into chromatin-assembled genomic DNA templates. Given a population of templates in which some undergo modification whereas others do not, the sensitivity of a given amplicon to DNase I will be reflected by a reduction in its apparent amplification efficiency from a DNase I–treated as compared to an untreated template population. Such differences can be readily detected using replicate quantitative real-time PCR measurements, which may be performed in a high-throughput format. The feasibility of using this protocol to measure the kinetics of DNase I digestion at known hypersensitive sites has been demonstrated previously8. However, accurate de novo localization of DNase I hypersensitive sites against a moving baseline of in vivo DNase I sensitivity presents a considerably greater challenge requiring the development of a quantitative methodology as outlined above.

QCP accurately delineates DNase I hypersensitive sites

To test these concepts, we produced a quantitative chromatin profile over a 21 kb interval spanning the human HBB (β-globin) locus control region (LCR), a paradigmatic complex cis-regulatory element9. The HBB LCR comprises an array of functional elements that in vivo give rise to five major DNase I hypersensitive sites (designated HS1–HS5) upstream of the HBE1 (ε-globin) gene on chromosome 11 (refs. 10, 11, 12). The core elements of these DNase I hypersensitive sites have been localized at the sequence level through extensive functional and chromatin structural studies13.

To produce a quantitative chromatin profile of the human HBB LCR, we prepared genomic DNA from DNase I–treated and untreated nuclei of K562 cells, an erythroid cell line in which the lineage-specific DNase I hypersensitive sites of the LCR are known to be active10. We designed a series of 89 contiguous 225 bp amplicons spanning the LCR, and then obtained nine replicate DNase I sensitivity measurements (ratios of the copy number in treated versus untreated samples) for each amplicon (801 total individual measurements). To arrive at a common DNase sensitivity scale that could be applied to other loci, relative copy ratios were normalized to a standardized reference amplicon from the rhodopsin gene locus on chromosome 3, which is transcriptionally inactive and resistant to DNase I in K562 and other cell types used herein (data not shown). We performed nine replicate measurements for each amplicon and determined both the trend and 'baseline' behavior of DNase I sensitivity across the locus. The measurement errors for DNase I sensitivity values clustered around the baseline, and hence confidence bounds on outliers and extreme values for this distribution. We then identified amplicons yielding replicate measurements having low variance (with respect to the mean measurement error) and a mean sensitivity ratio outside of the 95% confidence intervals about the baseline. To visualize the chromatin profile, we computed SNRs for each outlier relative to variability about the baseline and plotted these values on a genomic axis. The SNR is a broadly applied instrument in quantitative assays, in which numerical values reflect the number of standard deviations by which the observed measurement exceeds the expected background variability. Computed hypersensitivity SNRs across the HBB LCR are shown in Figure 2. Note that because SNRs are only shown for amplicons with significant (outlier) sensitivity ratios, the numerical range on the y-axis is arbitrary; in practice this range is scaled according to the distribution of SNR scores in the plotted segment.

Figure 2: Peaks in hypersensitivity SNR identify core functional elements.
figure 2

(a) Relative DNase I sensitivity measurements (DNase I–treated versus untreated; y-axis) over 25.6 kb spanning the HBB LCR (x-axis: chromosome 11 coordinates). (b) DNase I hypersensitivity in K562 cells expressed as computed SNR. HS1–HS5 are clearly evident; HS1′ identifies a previously described non-erythroid-specific HS31. Peaks in SNR coincide with the core functional elements of HS1–HS5 (red boxes) as localized through extensive in vivo studies. Notably, the amplicons corresponding to SNR peaks localize precisely to the core regulatory factor–binding regions of HS1–HS5 as determined by in vivo footprinting and other functional studies (see Supplementary Fig. 1 online).

All of the DNase I hypersensitive sites of the human HBB LCR were correctly and rigorously identified using this method. A weaker hypersensitive site described previously13 in K562 cells was also found (designated HS1′, Fig. 2). We then examined the correlation between the peaks in hypersensitivity SNR and the DNase I hypersensitive site core sequence. In each case, we found that the amplicon with the highest SNR precisely corresponded with sequences that had been determined previously to encompass the core regulatory factor binding domains of HS1–HS5 (Supplementary Fig. 1 online). These findings demonstrated the ability of QCP to rapidly recognize DNase I hypersensitive sites and to localize their core elements in a sequence-specific fashion.

High-throughout mapping of human gene regulatory regions

We then asked whether QCP could similarly expose the cis-regulatory architecture of other complex regulatory regions. Chromatin profiles were produced for the HBA@ (α-globin) locus upstream regulatory region on chromosome 16, the ADA locus (adenosine deaminase) on chromosome 20, CD2 locus (CD2 antigen, p50) on chromosome 1, the MYC locus (also known as c-myc) on chromosome 8, and the TCRA (T-cell receptor-α, TCR-α) locus downstream regulatory region on chromosome 14 (Figs. 3, 4, 5, 6).

We produced a chromatin profile spanning 66.4 kb upstream of the HBZ (ξ-globin) gene in the context of K562 cells (Fig. 3a). This region contains the HBA@ major regulatory region, which is situated in an intron of an unrelated upstream gene14. The hypersensitive sites identified by QCP encompass all of the major transcriptional control elements of the HBA@ upstream regulatory domain and all major hypersensitive sites previously reported following laborious studies of this region with conventional hypersensitivity assays14.

Figure 3: Identification of functionally diverse cis-regulatory sequences with QCP.
figure 3

Peaks in SNR (vertical axes) define functional sequences relative to genomic positions (horizontal axes). Genes (blue) are shown below genomic positions. (a) Quantitative chromatin profile over 66 kb spanning the HBA upstream regulatory region in erythroid cells (K562) using 2,439 DNase I sensitivity measurements from 271 amplicons. Analysis of the resulting profile revealed eight hypersensitive sites at −8 kb, −10 kb, −14 kb, −33 kb, −40 kb, −46 kb, −48 kb and −56 kb relative to the HBZ cap site. (b) Analysis of the newly identified −48 kb element (arrow) with conventional hypersensitive site assays in erythroid (K562) and lymphoid (Jurkat) cells. The parental band is a 10 kb XbaI fragment, and the probe (508 bp fragment, chromosome 16 94,343–94,922) is situated at the 5′ end. The new element shares the property of erythroid specificity with the two major elements of the locus (HS-40 and HS-33) and may therefore be important in α-globin regulation.

One of the major questions surrounding HBA@ locus regulation has been the location of additional erythroid-specific regulatory sequences that complement the activity of the major element at HS-40 and the auxiliary element at HS-33. Using QCP, we identified a site at −48 kb (Fig. 3a) that had not been detected in previous extensive surveys14. We confirmed that this site could be visualized with conventional assays (Fig. 3b), suggesting that it had been overlooked or was not detectable for technical reasons (such as a lack of appropriately positioned restriction sites). Like HS-40 and HS-33, this new element is erythroid specific, suggesting that it may contain the sought-after regulatory role that complements the other erythroid sites.

ADA is under control of a regulatory region positioned centrally within the 18 kb first intron15. The active core of this region is distinguished by the presence of a strong central tissue-specific DNase I hypersensitive site, which encodes a powerful transcriptional enhancer that also acts as a locus control region and is the dominant regulatory element of ADA. In the thymus, and in cell lines with a T-cell phenotype, this site is the only strongly evident hypersensitive site in the region. We performed QCP over the 58 kb ADA locus in Jurkat T cells and readily delineated this site, as well as the ADA promoter, which had not been previously analyzed at the chromatin level (Supplementary Fig. 2a online). Using the same preparation of Jurkat T cells, we applied QCP to a 27 kb region encompassing the CD2 locus. The resulting profile localized the CD2 promoter and the 3′ enhancer–LCR, as well as disclosing a previously unknown intronic element (Supplementary Fig. 2b online).

To demonstrate the ability of QCP to resolve complex, densely co-located elements, we studied the MYC locus on chromosome 8. MYC belongs to a class of highly regulated genes for which multiple distinct promoter elements have been described16, which are contained within hypersensitive sites in vivo17. A 30.4 kb chromatin profile performed in hepatocellular carcinoma (HepG2) cells revealed three clustered central sites, the SNR peaks of which corresponded precisely to the three core elements of the MYC promoter complex (Fig. 4).

Figure 4: QCP of the MYC locus (30.4 kb) on chromosome 8.
figure 4

Shown is a 30.4 kb chromatin profile derived from 855 measurements over 95 amplicons spanning the MYC locus, assayed in HepG2 hepatocellular carcinoma cells, in which MYC is transcriptionally active. Chromatin structural studies have identified three functional elements designated HS1–HS317. The MYC P1/2 promoter is contained in HSIII and the P0 promoter in HS2. We observed three clustered central sites, the SNR peaks of which corresponded precisely with the core elements of HS1–HS3. Two additional sites are evident, one at approximately −11 kb upstream and a second about 1.5 kb downstream of the gene; both correspond to previously documented elements32.

Delineation of cis elements across a multigene locus

We next asked whether the approach could resolve correctly, in a single experiment, the cis-active elements over a large multigene locus. To test this, we profiled the entire human β-like globin gene domain on chromosome 11, which comprises five genes (HBE, HBG2, HBG1, HBD and HBB, encoding ε-, Gγ-, Aγ-, δ- and β-globin, respectively) organized in a 5′-to-3′ fashion that corresponds to the timing of their expression during development and differentiation13. In addition to the LCR, known cis-active elements of the locus comprise the promoters of HBE, HBG2, HBG1, HBD and HBB, the HBG2 enhancer, and the 3′ HBB enhancer. We performed QCP over a 90 kb region encompassing the locus in K562 cells and localized all of the aforementioned elements (Fig. 5), as well as new elements. In particular, we observed an element distal to the HBB 3′ enhancer whose presence we also confirmed with conventional hypersensitivity assays (data not shown). This element lies within a 650 bp gap between two large repetitive tracts and coincides with a reported focus of intergenic transcription in erythroid cells18, suggesting a new regulatory role.

Figure 5: The 90.4 kb quantitative chromatin profile of the human β-globin locus in K562 cells.
figure 5

Shown (3′ to 5′) on the horizontal axis are the genomic positions of the HBE (ε), HBG1 (Gγ), HBG2 (Aγ), HBD (δ) and HBD (β) genes, as well as an olfactory receptor–like gene (OR5814) located 5′ of the LCR. The 90.4 kb (3,393 measurements over 377 amplicons) HBB quantitative chromatin profile from K562 cells revealed, in a single experiment, all of the major cis-regulatory elements of the globin locus, including the LCR (HS1–HS5); the HBE promoter together with an upstream element; the HBG1 promoter; the HBG2 promoter; the HBG2 3′ enhancer; the HBD promoter; the HBB promoter; and the HBB 3′ enhancer. The profile identified novel features (unlabeled peaks), including an element downstream of the HBB enhancer (arrow, far left) that coincides with a focus of intergenic transcription in erythroid cells18. Note the prominence of the HBD promoter, which is consistent with active δ-globin transcription in K562 cells33.

Analysis of the human TCRA regulatory domain

The human TCRA locus is embedded in a highly conserved segment of chromosome 14 and has not been analyzed systematically at the chromatin level. Knowledge of the cis-regulatory features of this locus derives largely from functional studies of the homologous mouse domain19, including the identification of a locus control region20 and an insulator activity21. In human cells, functional studies have disclosed an enhancer located 4.5 kb 3′ to the TCRA22.

We applied QCP to examine whether the spatial organization of the human regulatory region parallels that of the mouse. Comparative genomic analyses provide little insight into this question because of the high level of background noncoding sequence conservation (Fig. 6). We detected four prominent hypersensitive site core elements, which cleanly delineated the TCRA enhancer22, a juxtaposed insulator element, and the human homolog of the TCRA locus control region. A previously unindentified human element was evident within an intron of the downstream widely expressed DAD1 gene (Fig. 6). This example illustrates the power of a functionally based approach to illuminate and inform phylogenetic analyses.

Figure 6: Human TCRA (TCR-α) downstream regulatory region.
figure 6

(a) Quantitative chromatin profile (26.3 kb; 864 measurements over 96 amplicons assayed in Jurkat T cells) shows the spatial organization of the human TCRA regulatory region to be similar to that of the mouse domain, including the previously described human TCRA enhancer22, the human homolog of the mouse T cell α LCR, and an insulator element. A prominent hypersensitive site is also detected within an intron of the widely expressed DAD1 gene. Alignment of orthologous human and mouse sequences reveals extensive conservation across the TCRA regulatory region (below, % human-mouse conservation). Specific in vivo elements defined by QCP show varying degrees of conservation, but discrimination from other sequences in the locus showing similar levels of conservation is difficult. (b) Detection of DNase I hypersensitive sites (arrows) within the TCRA regulatory region in Jurkat cells. Parental band is a 9.2 kb BspHI fragment. (c) Detection of intronic DNase I hypersensitive site in DAD1 (arrow). The parental band is the 9.4 kb NdeI fragment. An additional hypersensitive site evident in the figure ('*') lies further downstream of the region profiled in a.

Accuracy of QCP

Of the 23 previously described classical cis-regulatory sequences expected to be active in the study cell types, 100% were successfully detected by QCP (Supplementary Table 1 online). Moreover, 29 of 29 previously described major human DNase I hypersensitive sites were also detected. We surveyed a total 1,167 amplicons in the best explored human genomic domains and identified a total of 41 SNR peaks, of which 33 coincided with previously defined sequence elements. Therefore, on the basis of this result alone (that is, not considering any additional validation studies, such as those presented in Figs. 4 and 6), the specificity of QCP (defined as the ratio of true-negative results to true-negative plus false-positive results) is at least 99.3% (= [1,167 − 41] / {[1,167 − 41] + [41 − 33]}).

Discussion

Sequence-specific localization of functional elements in the noncoding 98% of the human genome is a basic requirement for studies of gene regulation, for computational analysis of genomic sequences, for the interpretation of comparative genomic data and for systematic identification of a major class of functional genetic variation that has heretofore been largely obscure. We have shown that quantitative chromatin profiling may be applied to achieve precise localization of DNase I hypersensitive sites and, thereby, a diverse array of functional noncoding elements. QCP also revealed previously unidentified elements, which, surprisingly, were detected even in the context of the most heavily explored human genomic domains, the HBB and HBA@ loci.

Functional noncoding elements are rare in the human genome, and the critical test of any methodology applied to their detection over genomic distances is the specificity of the result. The specificity of QCP exceeds 99% and thereby substantially eclipses the performance of any described molecular, phylogenetic or computational approach. In addition to the accuracy and sequence specificity, a major strength of QCP is its applicability to diverse functional elements. This contrasts sharply with methodologies such as chromatin immunoprecipitation, which can only detect the presence of individual factors for which suitable antibodies are available.

The vast majority of human genetic variation resides in noncoding regions23. Genetic alterations in functional noncoding elements are expected to play a central role in the modulation of quantitative traits, including those that characterize the phenotypes of major common diseases. The search for such variation, however, has been complicated by the lack of methodologies for systematic exposition of functional noncoding elements within candidate gene loci. The ability of QCP to discriminate rapidly and accurately functional regions of the noncoding genome should have a substantial impact on our ability to discern functionally important alleles, identify heretofore elusive regulatory polymorphisms and dissect cis-acting contributors to phenotypic variation.

Methods

Cell culture and DNase I digestion.

See Supplementary Methods online.

Primer selection.

We designed primers to amplify contiguous or minimally overlapping 250 bp amplicons across target genomic regions. Primers were designed using Primer3 (ref. 24) restricting several parameters, including target amplimer size (250 bp ± 50 bases), primer Tm (optimal, 60 °C ± 2 °C), %GC (50% optimal, range 40–80%) and length (optimal 24, range 19–27), and the poly X (maximum 4). Primers were then scanned for repetitive sequences by alignment with a database of repetitive sequences.

Conventional DNase I assays.

Conventional DNase I hypersensitivity studies were performed using the indirect end-label technique5 according to a standard protocol as described previously25.

Measurement of relative DNase I sensitivity by quantitative PCR.

We assembled 15 μl real-time quantitative PCR reactions using 0.9 μM forward and reverse primers, 30 ng template DNA (untreated or DNase I treated) and master mix composed of 1 × FastStart buffer (Roche), 200 μM of each dATP, dCTP, dGTP, dTTP, 3 mM MgCl2 and FastStart Taq DNA polymerase (0.033 U/μl). The reaction mixture was supplemented with 0.33 × SYBR green I stain and 300 nM 6-ROX (Molecular Probes) to detect the accumulation of PCR product during amplification and to normalize fluorescence intensity, respectively. All quantitative PCR reactions were set up robotically with a Biomek FX (Beckman). Samples were run in triplicate on individual 384-well plates, and thermocycled with an ABI 7900HT Sequence Detection System (Applied Biosystems).

Primary data analysis.

We analyzed primary data in four phases: (i) determination of cycle threshold (Ct) values, (ii) amplification efficiency correction, (iii) melting curve analysis, and (iv) calculation of DNase I sensitivity ratios. Normalized fluorescence data were obtained and an amplification curve and nth-order polynomial fit were computed for each reaction. The Ct values were then determined for each curve. The amplification efficiency of a reference amplicon selected from the inactive and DNase-insensitive RHO (rhodopsin) locus (3q21–q24) was determined empirically for every reaction plate using a standard dilution series of DNA and the equation E = 10−1 / slope. We derived the efficiency of each test amplicon from the amplification curve26. Efficiency corrections were performed on all test amplicons with respect to the reference amplicon, and then we calculated relative copy number differences using the comparative Ct method27. Melting curve analysis was done for each amplicon to discard those yielding multiple products. Efficiency-corrected Ct values were then used to compute a relative copy number ratio by applying the formula 2−ΔΔCt or 2−[treated (target − reference) − calibrator (target − reference)]. Relative DNase I sensitivity ratios (relative copy ratios) were thus obtained. Ratios <1 were indicative of relative copy loss resulting from preferential cleavage of chromatin by DNase I.

Statistical analysis of DNase I sensitivity and hypersensitivity.

For a given genomic locus we expect, subject to inherent measurement error, a baseline response indicating general trends in DNase sensitivity across the locus. The deviation of repeated DNase I sensitivity measurements (at a given coordinate) from this baseline signals the presence of a hypersensitive site. The presence and significance of hypersensitive sites were quantified by using the following steps:

1. Identification of the trend or baseline behavior of a locus over an extended region.

2. Determination of the measurement error for DNase I sensitivity measurement values clustered around the baseline, and hence confidence bounds on outliers and extreme values for this distribution.

3. Identification of outliers that under repeated measurement have a clustering behavior or low variance with respect to the mean measurement error.

4. Assignment of a SNR to quantify the significance of this observation from the baseline.

An important observation that recurs throughout the analysis of DNase I sensitivity measurements (as well as other genomic data types such as DNA microarray data) is the non-Gaussian behavior of measurements. As the standard error increases, the ratio of measurements from Gaussian random variates approaches the Cauchy or Lorentz distribution. This has been shown to be the case for DNA microarray data28, where more robust methods for treating outliers are often necessary.

To determine the baseline of DNase I sensitivity measurements across a locus, a linear pass through the locus dataset is first performed using a 20% trimmed mean to remove glaring outliers. The remaining data set is smoothed using a locally weighted least squares (LOWESS) smoother that is adapted for robust locally weighted time series and scatter plot smoothing, as described29. A parameter (values 0–1) controls the size of the smoothing window and the standard value of 0.2 (which corresponding to using 20% of the locus in the fit) provides good results across a wide range of genomic loci. An important consideration in our choice of this method is its robustness in handing non-Gaussian time series.

To quantify the noise about the smooth baseline, data are mean centered about the moving baseline and about the outliers of this distribution using the median average deviation approach30. According to this method, a value X is declared to be an outlier if 0.6745 |XM| / MAD > 2.24 where M is the median and MAD is the average median deviation. Trimming the data in this way removes both the lower and upper extremes of the distribution, thereby addressing the problem of masking resulting from low sample number breakdown.

Clusters of hypersensitive site scores whose trimmed means occur outside the lower noise threshold from the baseline were then identified. To achieve this while keeping the analysis robust with respect to measurement error, another linear pass of the data is performed identifying groups at a common genomic position whose 20% trimmed mean lies strictly below the interpolated value at the lower shifted baseline. A small correction factor eliminates from consideration groups with very high variance or those consisting of a single point (zero variance).

Outliers corresponding to hypersensitive sites are further scored using SNR. To compute the SNR we compute S / Ni = |HSiBi| / MADB (1 + σC), where the SNR at site i is measured as the absolute deviation of the trimmed mean of the hypersensitive site cluster (HSi) from the interpolated baseline (Bi), divided by the median average deviation of the centered baseline. The remaining term in the denominator is a small correction factor that penalizes larger variances in hypersensitive site clusters and rewards highly compact clusters. The scoring method described above has merit in that it makes few distribution assumptions about the data and is robust under a relatively broad set of profiling scenarios.

Note: Supplementary information is available on the Nature Methods website.