The architecture of ChIA-PET Tool includes six modules: Linker filtering, PET mapping, PET classification, Binding site calling, Chromatin interaction calling, and ChIA-PET visualization (Figure ). First, in the Linker filtering module, the linkers from the input sequence reads are determined, and the PETs are sorted by the presence of readable linker barcodes. The PETs without readable linker barcodes are assigned as 'ambiguous PETs', whereas the PETs with readable barcodes are further assigned into chimeric PETs if they have heterogeneous linker barcode compositions (AB linker) or non-chimeric PETs if they have homogenous barcode compositions (AA or BB) (Figure ). Next, in the PET mapping module, the PET sequences are mapped to the corresponding reference genome. The mapped PETs are then classified based on the genomic spans of the two tag mapping locations into 'self-ligation PETs' (short genomic spans) and 'inter-ligation PETs' (long genomic spans or inter-chromosomal). The self-ligation PETs are used for calling putative binding sites, and the inter-ligation PETs are used for chromatin interaction analysis. The processed results are uploaded to a mySQL database for organization and visualization in ChIA-PET browser and the generic genome browser (G-browser) [10
Architecture of ChIA-PET Tool. The six main modules in ChIA-PET Tool (labeled) connect the input sequence data, the intermediate data, and the final results. The details of the processing are referred to in the main text.
To demonstrate the analysis procedure of ChIA-PET Tool, we used a real ChIA-PET library, IHH015A, for illustration. IHH015A is a part of the datasets of an ERα ChIA-PET study reported previously [9
]. This ChIA-PET library consists of 13,866,893 PETs generated by Illumina GAII paired-end sequencing, and was separated into chimeric PETs (IHH015C) and non-chimeric PETs (IHH015M) through the linker filtering procedure described below. The analysis results are summarized in Table and the remaining data analysis flow is mainly illustrated with the non-chimeric library IHH015M (Figure ).
Statistics of ChIA-PET data analysis
ChIA-PET data analysis flow of library IHH015A. A library, IHH015A, is used to demonstrate the ChIA-PET Tool data analysis flow, and the non-chimeric data IHH015M from IHH015A is used to show the analysis results.
A ChIA-PET input sequence has a pair of reads from the two opposite directions of the same ChIA-PET template. Both reads are 36 nucleotides long, and contain 20-nucleotide tag information and 16-nucleotide linker information (Figure ). To determine the linker type, we aligned the linker part of the reads (the last 16 nucleotide sequence of the 36 nucleotides) to the half-linker nucleotide sequence A (or B). If both the paired reads have good alignment with half-linker A (or B) and the specific nucleotide positions 9 and 10 (the nucleotide barcode) are CG (or AT), we classify the PET as having a homogenous full linker composition with AA (or BB). If one read in a PET has a good alignment with the half-linker A while the other one has a good alignment with the half-linker B, we classify this PET as having a heterogeneous full linker composition with AB (or BA). If there is no good alignment with any of the half-linkers (could also indicate low sequence quality), the PETs are classified as ambiguous PETs and will be discarded from further analyses. A PET sequence with a full linker AB indicates that this PET is derived from a ligation product formed between two different ChIP complexes from the two separate half-linker aliquots (Figure ). Therefore, the corresponding PETs with linker composition AB are most likely derived from random and non-specific ligations between two different ChIP complexes. Hence, we classified the PETs with linker AB as the chimeric PET dataset, and the PETs with linkers AA and BB as non-chimeric PET dataset (note that the PETs with linkers AA and BB may still have certain chimeric PETs). After the linker alignment, the linker parts in the short sequence reads will be trimmed and the tag sequences will proceed for further analyses.
With linker filtering, ChIA-PET library IHH015A data were divided into two pools: IHH015M (mix of PETs with AA and BB linkers) and IHH015C (chimeric PETs with AB linkers). The IHH015M dataset has 4,269,610 PETs (30.8% of total input PET sequences) and IHH015C has 6,157,038 PETs (44.4%). The remaining PETs were classified as PETs with ambiguous linkers.
PET mapping to reference genome
After linker trimming, the tags are mapped to the corresponding reference genome using the Batman package (C Tennakoon et al.
, manuscript submitted) with at most one mismatch. Batman is a Burrows-Wheeler-transform-based method [11
] that maps short sequences to a genome with very high speed. For each tag, Batman first considers the exact matches to the reference genome. If a single exact match is obtained, that location is taken as the mapping location of the tag, and the tag is classified as 'unique mapping'. If multiple exact matches are obtained, the tag is classified as 'multiple mappings'. If no exact match is obtained for the tag, a mapping is done with one mismatch allowed, and the result is similarly labeled as 'unique mapping' or 'multiple mappings'. If there is still no mapping for a tag with one mismatch, the tag is finally classified as 'non-mappable'.
After mapping the tags to the human reference genome (hg18) with Batman, only those PETs from IHH015M and IHH015C with both tags uniquely mapped to the reference genome were considered for further analyses. The remaining PETs with tags multiply mapped or unmapped to the reference genome were filtered out. We obtained 2,016,907 PETs in IHH015M and 2,707,860 PETs in IHH015C with unique mappings. To further avoid miscalling of clonal amplifications by PCR involved in sample preparation as enrichment of ChIP fragments, we merged all similarly mapped PETs (within ± 1 bp) into one unique PET. In this way, we reduced false positive calls resulting from the same PCR clonal amplification. Finally, we obtained 1,501,442 unique PETs with unique mapping from IHH015M and 1,710,335 unique PETs with unique mapping from IHH015C.
The ChIA-PET sequences can be classified into two categories: self-ligation and inter-ligation PETs (Figure ). 'Self-ligation PETs' are obtained from self-circularization ligation of the same chromatin fragments, and result in ChIA-PET sequences with both tags mapped within a short genomic distance of each other on the same chromosome in a head-to-tail orientation. 'Inter-ligation PETs' are derived from inter-ligation between two different DNA fragments, and can be partitioned into three different sub-categories: 'inter-chromosomal inter-ligation PETs', 'intra-chromosomal inter-ligation PETs', and 'different-orientation ligation PETs'. 'Inter-chromosomal inter-ligation PETs' are PETs with two tags mapped to two different chromosomes. 'Intra-chromosomal inter-ligation PETs' are PETs with both tags mapped to the same chromosome with a long genomic span, since PETs with long genomic span cannot arise from individual short chromatin fragments. 'Different-orientation ligation PETs' are PETs with both tags mapped to the same chromosome within a short genomic span, but with the wrong orientation or on different strands.
To determine the span cutoff between self-ligation and intra-chromosomal inter-ligation PETs, we plotted the genomic spans of the PETs mapped on the same chromosomes of the IHH015A dataset. The histogram shows that the vast majority of these PETs do not have genomic span over 2 kb (Figure ). Similarly, using log-log plot analysis of the same data (logarithm frequency against the logarithm span), we observed a mixture model with two straight distribution lines (Figure ), clearly representing two distinctive PET populations. The chromatin DNA was randomly chopped into short fragments (represented by self-ligation PETs) by sonication, which can be modeled by a 'stick-breaking process' [12
] for breaking long chromatin fibers. Our analysis and simulation suggest that the size distribution of the chromatin DNA fragments represented by self-ligation PETs follows a power-law distribution, which is a straight line in a log-log plot and corresponds to the left line in Figure . By contrast, the right line in Figure represents a PET population clearly different from the self-ligation PETs, which follows another power-law distribution and is an approximation for the intra-chromosomal chromatin interaction model as suggested by Dekker et al.
]. Therefore, the span cutoff between self-ligations and intra-chromosomal inter-ligations can be determined by the intersection of the two lines in the log-log plot. In our analysis for the IHH015M library data, the span cutoff called by our method is 4,595 bp. This estimation for DNA fragment size is consistent with agarose gel electrophoresis of the chromatin DNA sonication profile (Figure ) in the ChIA-PET protocol. The agarose gel result in Figure clearly shows that most DNA fragments are below the 1,650 bp mark and it is hard to see any DNA smear above 5,000 bp.
Figure 4 Span distribution of intra-chromosomal PETs and the cutoff between self-ligation PETs and inter-ligation PETs. (a) The distribution of the intra-chromosomal PET genomic spans. (b) The log-log plot of interaction frequencies against PET span, and the cutoff (more ...)
Accordingly, the PET datasets in IHH015M and IHH015C were classified into different categories. In IHH015M, 491,750 PETs were classified as self-ligation PETs. By contrast, IHH015C only yielded 12,155 self-ligation PETs. IHH015M contains 91,519 intra-chromosomal inter-ligation PETs and 893,213 inter-chromosomal inter-ligation PETs. IHH015C contains 94,743 intra-chromosomal inter-ligation PETs and 1,602,889 inter-chromosomal inter-ligation PETs. We note that the number of self-ligation PETs in IHH015M (with homogenous linker AA and BB; non-chimeric) is more than 40 times (491,750/12,155 = 40.5) that of the self-ligation PETs in IHH015C (with heterogeneous linker AB; chimeric). This is expected because self-ligations of the same DNA fragments are supposed to have the same linker types based on the experimental design (Figure ).
To check whether the inter-ligation PETs in IHH015M and IHH015C arise due to random ligation, the ratio of the intra-chromosomal inter-ligation PETs and the inter-chromosomal inter-ligation PETs was analyzed with a random model. From the PET datasets IHH015M and IHH015C, the numbers of DNA fragments from each chromosome were counted. Assuming that the ligation of the fragments occurs in a random manner and one fragment has an equal chance of ligating to any other fragments, the expected rate of finding an intra-chromosomal inter-ligation PET in a specific chromosome is proportional to the square of the number of DNA fragments in this chromosome. The total rate of intra-chromosomal inter-ligation PETs is proportional to the sum of the square of the numbers of DNA fragments over all individual chromosomes. Therefore, based on a random model, the expected rate of intra-chromosomal inter-ligation PETs is 0.0552 for IHH015C and 0.0546 for IHH015M. By contrast, the observed rate was 0.0558 for IHH015C and 0.0929 for IHH015M. The fold change between the observed rate and the expected rate for IHH015C (0.0558/0.0552 = 1.01; P-value 5.2E-4 from binomial test) was insignificant, validating the notion that the chimeric inter-ligation PETs are derived from random ligation. By contrast, the difference between the observed rate and the expected rate for IHH015M (0.0929/0.0546 = 1.70; P-value < 2.97E-323) was very significant, suggesting that the non-chimeric inter-ligation PETs are not random and probably enriched for specific chromatin interactions.
To visually illustrate the differences between the ChIA-PET libraries IHH015M and IHH015C, we represented every PET in these two datasets as a point (x, y) on a two-dimensional map where x and y axes are the genome loci of the two tags in a PET. Figure . shows the density map of all the PETs in the two ChIA-PET libraries, by normalizing the maximum density to 1. The darker the rectangle, the higher is the PET density between two corresponding chromosomes. For the chimeric PET dataset (IHH015C), there was no particular distribution pattern of PET density. By contrast, for the non-chimeric PET dataset (IHH015M), the PET density of the intra-chromosomal inter-ligations was much higher than that of the inter-chromosomal inter-ligations, suggesting that most potential chromatin interactions detected by ChIA-PET are intra-chromosomal, not inter-chromosomal.
Genome-wide ChIA-PET interaction density plots from the IHH015M and IHH015C PET datasets The square density plot in each graph is normalized to [0, 1]. Darker squares indicate higher interaction densities.
In summary, our analyses showed that non-chimeric PET dataset IHH015M is significantly different from the random model and the data can be used for further analyses.
Binding site analysis
Binding sites can be identified by looking for clusters of overlapping self-ligation PETs. As background noise is random, noisy PETs should distribute randomly throughout the genome. Only ChIP-enriched fragments would form clusters of overlapping self-ligation PETs, and are considered putative binding sites (Figure ). The probability that a cluster contains k
self-ligation PETs by random chance can be calculated by a Monte Carlo simulation, allowing false discovery rates to be assigned to different clusters. The clusters with false discovery rates below a particular threshold, such as 0.01, are putative binding sites. A similar method was previously applied in ChIP-PET data analysis [4
]. In considering the ChIP enrichment score of a binding site, the number of ChIP fragments at a binding site is counted. This is similar to the existing ChIP enrichment calculation protocols in ChIP-Seq data [5
After predicting binding sites, a post-processing step can be applied to remove suspicious binding sites, which can arise from different sources. For example, a library from female cells, such as MCF-7, is not expected to have any binding sites in chromosome Y. Binding sites in satellite repeat regions are also likely to be attributable to non-specific mapping and should be dropped. Further, certain binding sites in cancer genomes may be the results of genome amplifications [8
], and should be flagged, such that caution can be exercised in using them.
In our analysis, after calling putative binding sites and removing binding sites likely to be due to non-specific mapping, 4,035 binding sites were called from the IHH015M dataset at a false discovery rate <0.01, which covered 47,735 (9.7%) of the 491,750 self-ligation PETs in IHH015M. By contrast, only 24 binding sites were called from the chimeric PET dataset IHH015C at the same false discovery rate, which covered only 130 (1.1%) of 12,155 self-ligation PETs in IHH015C. This indicates that most of the binding sites of IHH015M are bona fide
. Indeed, many of the ERα binding sites identified in IHH015 library were validated using ChIP-qPCR [9
Chromatin interaction analysis
Chromatin interaction identification is done in two steps: first, define the ChIP enriched interaction anchor regions from inter-ligation PETs and then quantify the interaction frequency by counting the number of inter-ligation PETs between the two connected anchor regions. The identification of ChIP-enriched interaction anchors from inter-ligation PETs is performed by finding peaks and valleys from the overlapped tags from the inter-ligation PETs, in a similar manner employed by ChIP-Seq [15
]. The tag length of each inter-ligation PET is extended in a 5' to 3' manner by the 'tag extension length' to represent a ChIP DNA fragment. The 'tag extension length' is equivalent to the most frequently detected span of the self-ligation PETs, which is around 200 bp for the IHH015A library. Most of the interaction anchor regions identified from inter-ligation PETs should also be overlapped with the protein binding sites identified from self-ligation PETs. After defining the enriched anchor regions from inter-ligation PETs, the number of overlapping inter-ligation PETs between any two anchors is counted to reflect the relative interaction frequency. As each interaction involves two anchors and one loop, it is called a 'duplex interaction'. Similar to the binding site analysis, a real interaction is expected to involve multiple overlapping inter-ligation PETs connecting two anchors.
To determine if an interaction PET cluster between two anchors is a real chromatin interaction and not by random chance, a simple method is to count the number of inter-ligation PETs in the interaction cluster. If the cluster has a higher PET count, it has a higher probability to be a real chromatin interaction. This model, however, does not take into account the fact that, when the ChIP enrichments of two anchors are high, there is also a higher probability to form more inter-ligation PETs by random chance between these two anchors. To address such concerns, we formulated a statistical analysis framework to account for the random formation of any inter-ligation PETs between two anchors. The null hypothesis assumes that, in the ChIP-enriched chromatin fragment population, each chromatin fragment has an equal chance to ligate to any other fragments in a random manner, and the interactions between different anchors are independent of each other. Under this random model, the number of inter-ligation PETs that link two anchors follows a hyper-geometric distribution. The formula is provided in Equation 1:
The expected interaction frequencies between any two genomic loci and the false discovery rate for each interaction were calculated. Hence, both inter-ligation PET frequency and ChIP enrichment of the anchors are taken into account by this analysis.
Equation 1 considers a library with N inter-ligation PETs. RA and RB represent two anchor regions with cA and cB PETs, respectively, where cA, cB <<N. Equation 1 shows that, when cB ends are randomly chosen from 2N ends as ends in region RB, what is the probability of choosing IA, B ends from cA ends of region RA to form IA, B interactions between region RA and region RB. By this, we are able to compute a P-value to test if IA, B, the number of inter-ligation PETs between RA and RB, is over-represented. Given a cut-off threshold, T, of hypergeometric P-value, we are able to calculate the false discovery rate, which is the fraction of the clusters with P-value below T under the empirical random model generated by randomly permuting the ends of PETs. As cA and cB reflect the ChIP enrichment of two DNA anchors, any ChIP enrichment bias is accounted for by this model. An illustration of the random model is shown in Figure .
Figure 6 An illustration of the statistical model for probability analysis of ChIA-PET interactions. RA and RB represent two DNA regions ('anchors') with cA and cB PETs (here 9 and 7, respectively). IA, B is the number of inter-ligation PETs between RA and RB (more ...)
From the predicted interactions with three or more inter-ligation PETs between anchors and false discovery rate < 0.05, we filtered the interactions that could be a result of mis-mapping or noise, as we did with the binding sites. We filtered out any interactions wherein the anchors were present in chromosome Y or the mitochondria, and satellite repeat regions. We also filtered out a specific kind of noise in the interaction clusters from genome structural variations in cancer cell-lines. The cancer cell-lines like MCF-7 have intensive genome rearrangements: genome insertion, deletion, inversion, and translocation, which can be identified with DNA-PET data [8
]. Self-ligation PETs around the breakpoints of genome rearrangement from cancer cell-lines will be mapped as inter-ligation PETs in the reference genome, and interaction clusters related to such genome rearrangements should be removed.
Using the above analysis procedure and a threshold of 3 or more for the number of overlapping inter-ligation PETs, we identified 1,945 putative intra-chromosomal interactions and 12 putative inter-chromosomal interactions in the IHH015M dataset. Our validations, including 3C [13
], ChIP-3C [16
], and 4C [17
], as well as fluorescent in situ
hybridization (FISH) and small interfering RNA experiments, suggest that the majority of the intra-chromosomal chromatin interactions identified in this analysis are bona fide
loci bound by ERα as reported in Fullwood et al.
]. By contrast, the chimeric inter-ligation PETs in IHH015C yielded only 12 intra-chromosomal and zero inter-chromosomal inter-ligation PET clusters. Detailed manual curation verified that none of them constitute real interactions. Our comparison of a non-chimeric ChIA-PET dataset (IHH015M) with a similarly sized chimeric ChIA-PET dataset also indicates that the non-chimeric dataset generates statistically significant binding sites and interactions, whereas the chimeric dataset does not. This means that chimerism is not an issue in obtaining bona fide
binding sites and chromatin interactions in the ChIA-PET library. As an example, abundant non-chimeric inter-ligation PETs in the IHH015M dataset identified the ERα-bound chromatin interaction event at the KRT8/18 locus in the human genome (Figure .), but no chimeric PETs in the IHH015C dataset were found at the KRT8/18 locus. This interaction at the KRT8/18 locus has been validated by 4C experiments [9
Figure 7 An example of ChIA-PET mapping display in G-browser. A screenshot of the ChIA-PET G-browser is shown in human chromosome 12 at the keratin gene family region. The tracks are (from top to bottom): UCSC known genes track, which shows the KRT8 and KRT18 (more ...)
ChIA-PET data visualization
All ChIA-PET data, including the called binding sites and chromatin interaction clusters, are uploaded to a mySQL database. A centralized ChIA-PET browser is built to organize data reporting and visualization. The ChIA-PET browser consists of two components: a tabular browser and a graphic genome browser. The tabular browser provides an organization of all ChIA-PET experimental datasets (libraries) from a variety of projects. It reports the unique PETs with unique mapping, the binding sites and the interaction clusters in tabular forms, and the users can download the data for further analysis. Example screenshots of binding sites, interaction clusters and the whole genome interactions from IHH015M are provided in Figures , and . The graphical genome browser (G-browser) is created by adopting the 'generic genome browser system' [10
], which allows users to view and manually curate the data. A screenshot of the G-browser is shown in Figure More details of the browsers can be found from the ChIA-PET website [18
] (username, guest; password, guest).
Screen shot of binding site table view in ChIA-PET browser.
Screen shot of interaction cluster table view in ChIA-PET browser.
Screen shot of the whole genome interaction view in ChIA-PET browser.
Implementation and performance
ChIA-PET Tool is implemented with various software written in C, Java and scripting languages (PERL and Python), and has been tested with the Linux operating system. The components in the pipeline are linked up together to behave as a singular processing pipeline. The pipeline requires hardware support; for example, 32 Gbytes RAM and 1,024 Gigabytes hard disk are recommended on a duo-core server. mySQL (5.0.67 or higher) is used as the database engine. Other required software packages are Apache, Perl and Bioperl modules, and PHP (refer to the ChIA-PET Tool installation guide on the ChIA-PET website [18
] for details).
We have tested the ChIA-PET Tool pipeline with two hardware configurations: one with 32 Gbytes RAM and 8 CPUs, and another with 64 Gbytes RAM and 16 CPUs. The input to the pipeline is the original paired-end sequence data from the Illumina sequencing output file. A ChIA-PET library from a single lane requires approximately 1 Gigabyte space for storage, which includes intermediate files generated throughout the processing. On average, it takes 1 hour for ChIA-PET Tool to process a single lane ChIA-PET dataset and 1 hour to upload the GFF file to the G-browser database (the upload time is subject to existing database size and server load).