ChIPulate: A comprehensive ChIP-seq simulation pipeline
[['ChIP-seq ( Chromatin Immunoprecipitation and sequencing ) is a popular high-throughput experimental technique to find locations that are bound in vivo by a single transcription factor ( TF ) [1] .', 'Upon mapping of the DNA fragments bound by the TF to the reference genome , the genomic loci bound by the TF are identified as high density mapped regions or peaks , where each peak is associated with an intensity based on the number of sequenced fragments arising from it .', 'The intensity reflects the in vivo occupancy of the TF at that locus .', 'Several studies of ChIP-seq data have focussed on the biological factors distinguishing the loci bound by the TF .', 'It has been shown that in addition to the affinities of binding sites present at a locus , nucleosome positioning is a strong determinant of TF binding in vivo [2–5] .', 'Other studies have shown that the concentration of the target TF [6 , 7] , short-range cooperative interactions between the target TF and other TFs [8] , and variation in chromatin accessibility [5 , 7] explain the variation in intensities across peaks .', 'Some of the variation can arise due to indirect binding , where the target TF binds DNA indirectly via a second DNA-bound TF [9–11] .', 'The intensity of such peaks is then no longer directly dependent on the affinity of the target TF to sequence at the bound locus .', 'Since the distribution of ChIP-seq peaks and their intensities depend on factors other than the affinity of the target TF towards sequence at a locus , it impacts two kinds of inferences often made from ChIP-seq data .', 'First , the highest intensity peaks are used to infer position weight matrix ( PWM ) or motif models of the target TF [12] .', 'While it is known that changes in the concentration of the target TF can distort the inferred PWM [13 , 14] , the extent to which cooperative and indirect interactions distort the PWM is unclear .', 'Second , along with variations in chromatin accessibility , these interactions weaken the statistical dependence of peak intensity on target TF binding site affinity , which means that a peak with a higher intensity need not necessarily contain a higher affinity target TF binding site .', 'In addition to these biological factors , ChIP-seq peaks are affected by purely experimental sources of noise .', 'The ChIP-seq protocol broadly consists of three key steps [1 , 15 , 16]—', '( i ) the extraction of fragments that are bound by the target TF', '( ii ) PCR amplification of these extracted fragments , and', '( iii ) sequencing of these fragments .', 'It is known that fragments can be extracted more easily at some genomic regions than others due to differences in chromatin structure and cross-linking efficiency across the genome [17–19] .', 'Similarly , certain fragments are more efficiently amplified by PCR than others due to differences in GC content or the presence of nucleotide repeats [20–22] .', 'The extent to which these factors distort peak intensities and make otherwise identical genomic regions appear differentially occupied has not been quantified .', 'Multiple biological replicates of ChIP-seq are recommended [23] to overcome these issues , but the quantitative improvement in accuracy as more replicates are performed has not been ascertained .', 'To evaluate the influence of the aforementioned biological and experimental sources of variation on ChIP-seq peak intensities , we have developed a comprehensive pipeline to simulate a ChIP-seq experiment called ChIPulate .', 'In particular , we simulate genome-wide TF-DNA binding based on established biophysical models of binding and also simulate the extraction , amplification and sequencing steps of the ChIP-seq protocol .', 'Our simulation thus associates a set of read counts at each genomic locus with its occupancy by the target TF .', 'Whereas existing statistical models [14 , 24–26] assume that read counts follow a prescribed probability distribution and treat the ChIP-seq protocol as a black box , we explicitly model key steps of the ChIP-seq protocol .', 'Our approach allows us to individually evaluate the impact of variation in extraction and PCR amplification efficiencies , as well as chromatin accessibility , indirect binding , and cooperative binding , on peak intensities and PWM inference .', 'We are also able to quantify the extent to which additional replicates of ChIP-seq improve its ability to robustly measure the occupancy of a genomic region .', 'We find that biological factors such as indirectly and cooperatively bound sequences distort inferred PWMs more than experimental sources of variation .', 'Variations in extraction efficiency across the genome distort peak intensities and lower their ability to distinguish between occupied loci .', 'Poor extraction efficiency also increases the number of false positive peaks , which are peak calls that do not contain a binding site for the target TF .', 'In contrast to the effect of variations in extraction efficiency , even drastic variations in PCR amplification efficiency have little impact on peak intensities , and hence do not affect the inferred PWM or increase the number of false positive peaks .', 'Finally , we found that at least two biological replicates of ChIP-seq read counts are necessary to reliably infer the binding energy of a genomic region .', 'Our work provides a general framework and a software tool for simulating ChIP-seq read counts through a realistic model of TF-DNA binding and the steps of the ChIP-seq protocol .', 'Improvements in the protocol to lower variation in extraction efficiency , rather than PCR amplification efficiency , are more likely to improve our ability to distinguish between genomic regions of differing occupancy .', 'Further , changes in the protocol that allow multiple biological ChIP-seq replicates to be performed , or computational approaches that can reliably combine read counts from ChIP-seq experiments , would allow more accurate inferences of TF-DNA occupancy , especially in regions containing multiple binding sites .'], ['A schematic of our ChIP-seq simulation framework is shown in Fig 1A , with details in Methods .', 'Broadly , the goal of our simulation is to take as an input the occupancy of a TF at multiple locations across a genome and output a set of read counts in both ChIP and input experiments at each location .', 'A single TF’s occupancy at a genomic location is determined by its site-specific binding energy , and its chemical potential , which depends on the logarithm of the concentration of the TF [27] .', 'Each genomic location is associated with two experimental parameters—its extraction efficiency and PCR efficiency .', 'The extraction efficiency , which determines how many bound fragments are successfully extracted in Step 1 in Fig 1A , can vary between different genomic locations .', 'The extracted fragments are then subjected to a number of PCR amplification cycles in Step 2 .', 'The average number of amplified fragments obtained at each location depends on the PCR efficiency associated with the locus , which can also vary between different locations .', 'Step 3 represents the sequencing step where the amplified fragments from all locations are sequenced .', 'The number of sequenced fragments is the total read count of the experiment .', 'Since unique and duplicate fragments ( red and black in Fig 1A ) are tracked through the simulation , the exact number of unique reads at each location is known .', 'In addition to the ChIP experiment , we simulate the sequencing of a genomic control input sample .', 'This simulation differs from the ChIP experiment in two ways .', 'First , the background binding energy ( ϵbg ) at each occupied genomic locus is set to a fixed value of 1kB T for all locations ( see Methods for justification ) .', 'Second , we set the number of cells to be only 10% of the number of cells employed in the ChIP simulation .', 'A larger number of cells are used in the ChIP simulation since , in practice , DNA is less efficiently extracted in the ChIP sample due to the use of an antibody whereas no antibody is used in the input sample .', 'Thus , to extract a given amount of DNA , fewer cells are employed in the input sample than in the ChIP sample .', 'In the event that the number of extracted fragments in the input is more than the number of extracted fragments in the ChIP sample , we down-sample the input fragments such that the total number of extracted fragments is identical in both ChIP and input samples before the PCR amplification step .', 'Similarly , we down-sample the ChIP fragments if the number of extracted fragments in the ChIP sample is more than than that in the input sample .', 'At each locus , the extraction and PCR efficiencies in the input are kept identical to the values used in the ChIP simulation .', 'After simulating both ChIP and input experiments , we compute the ratio of the read counts obtained from the ChIP sample to that of the input sample at each genomic locus .', 'The read count ratio at a location is equivalent to peak intensity computed by commonly used peak callers from ChIP-seq reads [24] .', 'The read counts obtained from our simulations are shown in Fig 1B with the values of the simulation parameters set to the default values described in Methods .', 'When simulating genome-wide TF-DNA binding before Step 1 in Fig 1A , we assume that there is a single TF and a single binding site in each genomic region .', 'We relax this assumption in later sections when simulating indirect and cooperative binding .', 'We assume that binding in different genomic locations are independent of each other .', 'The extraction efficiency parameter in Step 1 aggregates the effect of the many steps of fragmentation , cross-linking , pull-down and size selection that occur in a ChIP-seq experimental protocol .', 'We show in Methods that this parameter takes all these steps into account when modeling the extraction of DNA fragments from cells .', 'In the PCR amplification step , we assume the amplification process to be in the exponential phase of amplification [28] where increasing the number of PCR cycles exponentially increases the number of amplified fragments .', 'We assume that any potential PCR mutational errors do not change the amplification ratio at any genomic location .', 'In the sequencing step ( Step 3 ) , we assume that the sequencing and mapping are error-free , i . e . , reads are uniquely and correctly mapped , and they can be perfectly de-duplicated at each location .', 'We evaluated the impact of heterogeneity in extraction and PCR efficiencies on two outcomes of ChIP-seq .', 'One is the recoverability of TF motif based on the top genomic locations in terms of their read count ratios .', 'The second is fidelity , defined as the probability ( frequency ) with which a site X with a read count ratio that is at least 10% higher than site Y in fact contains a sequence with a higher affinity ( or lower binding energy ) .', 'We simulated the process of inferring the motif of an arbitrarily chosen S . cerevisiae transcription factor ( Tye7 ) from ChIP-seq ( Fig 2A ) .', 'We assigned binding energies for 1000 genomic regions from our default binding energy distribution , which is a power law distribution ( with an exponent of 0 . 5 ) between 0 and 10kB T ( see Methods ) .', 'By using the Tye7 binding energy matrix estimated by BEEML from protein binding microarray measurements [13 , 29] , we found the 10 bp sequence whose binding energy was closest to this assigned value , and virtually planted it at each location .', 'We simulated a ChIP-seq experiment with a fixed extraction and PCR efficiency across the genome , following which we selected sites having the top 100 read count ratios and constructed the PWM of Tye7 .', 'This PWM , which we refer to as the baseline PWM , is used to compute the extent to which heterogeneity in extraction and PCR efficiency changes the derived PWM .', 'Ideally , one would use a PWM of Tye7 published in a database as the baseline motif against which to make these comparisons .', 'However , these motifs can differ between databases due to differences not only in the experimental assay used to measure TF-DNA binding but also the underlying binding energy distribution of the bound sequences .', 'This means that the baseline Tye7 PWM derived from the simulation may differ from that present in the ScerTF database [30] .', 'However , our goal is to specifically assess the effect of extraction and PCR amplification heterogeneity on motif recoverability .', 'For this purpose , it suffices to compare the motif derived in the presence of heterogeneity with the baseline motif derived in the absence of heterogeneity .', 'Fig 2B shows the difference between the baseline PWM and the Tye7 PWM that is derived in the presence of extraction or PCR amplification heterogeneity , respectively , when the efficiencies follow a truncated normal distribution .', 'This difference is measured in terms of the mean Kullback-Leibler ( K-L ) distance , measured in bits ( Methods ) , between the two PWMs .', 'The impact of heterogeneity in these parameters on the inferred motif is small as the highest K-L distance from the baseline motif is about 0 . 08 bits .', 'The effect of PCR amplification and extraction heterogeneity on motif inference is more pronounced when the heterogeneity follows a power law distribution ( S1A Fig ) than when it follows a truncated normal distribution .', 'The distortion of the PWM is still relatively low , however , with the highest K-L distance from the baseline motif being less than 0 . 1 bits per position .', 'We also evaluated the impact of varying the chemical potential between 1kB T and 6kB T on the inferred motif when the baseline motif is derived at μ = 3kB T . Examples of the effect of a different chemical potential on the read count ratios obtained across the genome are shown in S2A Fig . We found that the inferred motif deviates up to 0 . 08 bits per position from the baseline motif ( S2B Fig ) .', 'This is within the same order of magnitude as the effect of a low mean extraction efficiency on the inferred motif .', 'We checked the extent to which the trends in Fig 2B depended on our choice of TF .', 'We selected 21 TFs that possessed different DNA binding domains [31] ( S1 Table ) , and for each TF , we simulated motif inference when there is heterogeneity in extraction efficiency ( S3 Fig ) or the PCR amplification ratio ( S4 Fig ) across the genome .', 'The binding energy distribution of each TF , along with all other parameters of ChIPulate , was kept at their default values .', 'We found that motif inference in each of these TFs was affected in a manner similar to Tye7 , where a low mean extraction efficiency distorted the inferred motif more than a low mean amplification ratio .', 'However , we found that the extent of motif distortion due to extraction heterogeneity was significantly correlated ( R2 = −0 . 67 , p < 10−3 ) with the baseline motif information content of each TF ( S3A Fig ) , while there was no such correlation in the presence of amplification heterogeneity ( S4A Fig , R2 = −0 . 36 , p = 0 . 1 ) .', 'Finally , we found a weak dependence of the extent of motif distortion on the length of the motif , where a low mean extraction efficiency distorted a 6 bp motif upto nearly 0 . 1 bits per position with a lower magnitude of distortion for longer motifs ( S5A Fig ) .', 'There was no such dependence of motif distortion on motif length when the amplification ratio varied across the genome ( S5B Fig ) .', 'Results shown in Fig 2C show that ChIP-seq fidelity is generally higher for sites with intermediate read count ratios ( 50th percentile to 75th percentile ) and deteriorates at lower and higher read count ratios .', 'This is partly due to the value of the chemical potential of 3kB T employed in our simulation .', 'When the chemical potential is varied between 1kB T and 6kB T , the fidelity amongst the top 25th percentile of read count ratios noticeably decreases and reaches close to 0 . 5 at μ = 6kB T ( S2C Fig ) .', 'The fidelity amongst read count ratios of the bottom 25% , however , increases over the same range of the chemical potential .', 'This is due to the fact that the occupancy of locations whose binding energies are much lower than the chemical potential are close to 1 while locations with binding energies above the chemical potential are close to zero .', 'As the chemical potential increases , the increase in occupancy at locations with high binding energies is sufficiently large for their read count ratios to differ from one another .', 'On the other hand , the occupancies of locations with low binding energies rapidly tend to one , which results in similar read count ratios between sites that possess different binding energies .', 'Fig 2D shows the impact of heterogeneity in extraction efficiency and PCR amplification ratio on fidelity .', 'The fidelity is highest when computed across all binding site pairs , than amongst pairs in the top 25% of read count ratios .', 'This is because the expected difference between the binding energies of sites in each pair is larger when the pairs are selected from the entire range .', 'Changes in the mean and variance in the amplification ratio have little impact on fidelity , both overall as well as in each read count ratio bin .', 'The fidelity across all regions is also not substantially affected by the mean and variance in extraction efficiency , but the fidelity is markedly lowered for sites with the lowest read count ratios .', 'When PCR amplification heterogeneity is power-law distributed , its impact on fidelity is still quite low ( S1B Fig ) .', 'In contrast , when the extraction efficiency is power-law distributed across the genome , its impact on fidelity is more drastic than when it is normally distributed .', 'We tested the effect of heterogeneity in extraction and PCR efficiency on the ability of ChIP-seq to distinguish peak calls that contain a binding site for the target TF ( true positive ) from those peak calls that do not contain a binding site for the target TF ( false positive ) .', 'Fig 3A shows a simulation where genomic locations that do not harbor a binding site for the target TF give rise to sequence reads .', 'We refer to such locations as false positive peaks .', 'This differs from indirectly bound peaks , where the target TF is part of the complex of TFs that binds DNA .', 'In practice , a false positive peak can arise due to poor antibody quality , where the antibody binds epitopes on TFs other than the target TF .', 'We assign binding energies to these false positive locations from a truncated power law such that the mean occupancy ratio between true and false positive sites are fixed while the variance in binding energies amongst false positive sites and true positive sites are equal .', 'Fig 3B shows the receiver operating characteristic ( ROC ) curve for the simulation shown in Fig 3A .', 'The ROC curve shows the change in the number of true positives and false positives when the read count ratio threshold at which a location is declared as a peak is changed .', 'The area under the ROC curve ( auROC ) provides an overall measure of the accuracy of using ChIP-seq read count ratios in distinguishing between locations with target TF binding sites from those without these sites .', 'Fig 3C and 3D show the impact of heterogeneity in extraction efficiency and PCR efficiency , respectively , on the auROC of ChIP-seq read count ratios .', 'We found that increases in the mean PCR efficiency did not appreciably increase the auROC .', 'Similarly , an increase in the heterogeneity of PCR efficiency , in terms of its coefficient of variation across the genome , did not decrease the auROC .', 'In contrast , an increase in the mean extraction efficiency , and a decrease in its heterogeneity , led to a higher auROC .', 'Thus , much like in the case of ChIP-seq fidelity and motif inference , extraction efficiency has a greater impact on the auROC of ChIP-seq read count ratios than PCR efficiency .', 'Fig 3E shows that as expected , the auROC increases with the ratio of the mean occupancy of true and false positive binding sites .', 'We illustrate in Methods that this mean occupancy ratio can be considered as a proxy for the specificity of the ChIP antibody used , where a higher mean occupancy ratio translates to a higher antibody specificity .', 'A mean occupancy ratio of ∼10 gives an auROC of around 0 . 9 , which translates to a true positive rate of 0 . 8 at a false positive rate of 0 . 07 .', 'For a fixed mean occupancy ratio , an increase in the sequencing depth of both ChIP and input samples leads to an increase in the auROC ( Fig 3F ) .', 'We further evaluated the impact of variation in chromatin accessibility across true and false positive binding sites .', 'Our model of TF-DNA binding in the presence of chromatin assumes that TFs occupy a genomic location only when chromatin is accessible , i . e . , the probability of occupancy depends on the energy of the binding site present and the probability that chromatin is accessible at that location ( Eq ( 8 ) ) .', 'We found that chromatin accessibility had a low impact on the auROC when the mean accessibility was at least 0 . 2 ( S6 Fig ) .', 'We note that Eq ( 8 ) implicitly assumes that the probability of occupancy at a genomic location is zero in cells where chromatin is inaccessible .', 'This does not hold for TFs such as pioneer TFs [32] , which are capable of binding inaccessible chromatin and opening it up for further binding by other TFs .', 'For these TFs , our model under-estimates the occupancy probability at a given genomic location , with the result that the auROC that we obtain in our simulation represents a lower bound .', 'The target TF ( say , A ) of a ChIP-seq experiment is said to be indirectly bound to DNA if it does not bind DNA directly but instead binds a second TF ( say , B ) that in turn directly binds DNA .', 'Thus , an indirectly bound locus lacks a binding site for the target TF , and the number of reads mapped to this location do not depend on the binding energy of A . This can be seen in the top panel of Fig 4A , where the read count ratios of locations directly bound by A depend on the binding energy of A ( R2 = −0 . 73 , p < 10−115 ) whereas the read counts from indirectly bound regions do not ( R2 = 0 . 004 , p = 0 . 933 ) .', 'In contrast to indirect binding , genomic locations that are cooperatively bound by the target TF and the second TF contain a binding site for each TF .', 'The magnitude of the cooperative effect is represented by the interaction energy ( Δ , in units of kB T ) .', 'Lower ( negative ) values of Δ represent a larger cooperative effect , with Δ = 0 representing independent binding ( see Eq ( 4 ) ) .', 'In the bottom panel of Fig 4A , read count ratios of some cooperatively bound locations ( Δ = −6kB T ) are higher than independently bound regions .', 'Unlike in the case of indirect binding , however , read count ratios from cooperatively bound regions depend on the binding energy of the target TF ( R2 = −0 . 66 , p < 10−38 ) , in addition to depending on the binding energy of the second TF and the interaction energy between both TFs .', 'In Fig 4B , we separately evaluated the impact of indirect and cooperative binding on the inferred PWM of the target TF .', 'A detailed description of PWM inference in the presence of indirect and cooperative binding is given in Methods .', 'The baseline PWM of the TF Tye7 is first computed with direct binding alone .', 'We then computed the mean K-L distance between this baseline PWM and the PWM inferred in the presence of either cooperative binding or indirect binding with a second TF .', 'We found that the inclusion of a given fraction of indirectly bound sequences distorted the inferred PWM more than the inclusion of the same fraction of cooperatively bound sequences .', 'The divergence of the inferred PWM from the baseline depended both on the binding energy range of the indirectly or cooperatively binding second TF , as well as the strength of the cooperative effect .', 'Similarly , we found that a given fraction of indirectly bound locations had a greater impact on ChIP-seq fidelity than the same fraction of cooperatively bound locations ( Fig 4C ) .', 'This was true even though the cooperativity is high ( Δ = −6kB T ) .', 'Indirect and cooperative binding particularly affected the fidelity of read count ratios from the strongest binding sites ( the top 25th percentile of read count ratios ) .', 'Next , we assessed the accuracy with which the binding energy at a locus can be estimated using ChIP-seq read count ratios , and the impact of replicates on this accuracy .', 'If we denote the true binding energy of the i-th locus as ϵ ( i ) and the estimate of the binding energy as ϵ ¯ ( i ) , we expect ϵ ¯ ( i ) to move closer to ϵ ( i ) , on average , as more replicates of ChIP-seq are performed .', 'The first step to quantifying the accuracy of ϵ ¯ ( i ) is to compute its posterior probability after having observed read count ratios R ( 1 ) , R ( 2 ) , … , R ( n ) from n biological replicates , which is the conditional probability P ( ϵ ¯ ( i ) = ϵ | R ( 1 ) , R ( 2 ) , … , R ( n ) ) .', 'We define the absolute uncertainty or uncertainty of the binding energy estimate to be the 95% Bayesian credible interval of this posterior distribution [33] , which we take to be the difference between its 97 . 5th and 2 . 5th quantiles .', 'The posterior probability , and hence , the uncertainty of ϵ ¯ ( i ) after n + 1 replicates are performed , can be computed from P ( ϵ ¯ ( i ) = ϵ | R ( 1 ) , R ( 2 ) , … , R ( n ) ) using Bayes’ rule , as shown in Methods .', 'The expected decrease in the uncertainty of ϵ ¯ ( i ) as more replicates are performed can be seen in Fig 5A , where the area under posterior probability of ϵ ¯ ( i ) becomes more concentrated around ϵ ( i ) with each additional replicate .', 'We repeated this procedure to calculate the reduction in uncertainty at locations with binding energies between 2kB T and 6kB T in Fig 5B , and found that the uncertainty is high at sites with the lowest binding energies even after five replicates are performed .', 'At the lowest binding energies , we also checked if the uncertainty after five replicates sufficed to distinguish between two binding sites that differ by one base pair from each other .', 'From an analysis of binding energy matrices of 368 TFs in the BEEML database [13] , we estimated that 75% of single mutations change the binding energy of a site by at least 0 . 24kB T , while 50% of single mutations change the energy by at least 0 . 71kB T . Thus , even three replicates of ChIP-seq would allow us to distinguish between two binding sites that are a single mutation apart with a probability between 50% and 75% at only relatively low energy binding sites .', 'In addition to simulating read counts from ChIP and control samples as described so far , ChIPulate can output single-end and paired-end reads for each genomic location if its chromosome coordinates are specified .', 'Fig 6A illustrates the working of ChIPulate when it is run in FASTQ generation mode .', 'To generate sequence reads , a set of n genomic intervals ( b1 , e1 ) , ( b2 , e2 ) , … , ( bn , en ) , along with a set of “summits” s1 , s2 , … , sn is input to ChIPulate .', 'These summits can be thought of as the physical location of TF-DNA binding in each interval .', 'The total and unique read counts ( ri and ui , respectively ) at each interval are computed as described in Methods based on the binding energy ( or energies ) associated with each interval .', 'Fragments that are PCR duplicates are tracked and the reads arising from duplicate fragments are named appropriately to reflect their status as a PCR duplicate .', 'In the ChIP sample , the start positions of fragments in the i-th interval are sampled from a Gaussian distribution with mean si − d/2 and standard deviation j , where d is the fragment length and j is referred to as the fragment jitter .', 'This captures the notion that a genomic position close to a TF-DNA binding event is more likely to give rise to a fragment in the ChIP sample than a genomic position that is far away .', 'In the control sample , each position in an interval is assumed to be equally likely to give rise to a fragment .', 'This corresponds to the Poisson background model assumed by peak callers such as MACS2 [24] .', 'Thus , in the control sample , ChIPulate samples start positions in the i-th genomic interval from a Uniform ( bi , ei ) distribution .', 'The l base pairs at the ends of each fragment in both ChIP and control samples represents the read length .', 'Though Fig 6A illustrates the simulation of paired-end reads , ChIPulate can simulate both single-end reads and paired-end libraries , where the sequence of each read is assumed to be error-free .', 'Fig 6B is an IGV ( Integrated Genomics Viewer ) [35] snapshot of paired-end reads simulated from a genomic interval containing a GCN4 binding site in the S . cerevisiae genome .', 'In this example , we chose a set of 549 genomic regions that were found to be bound by GCN4 in an earlier publication [34] , and computed the binding energy of GCN4 in each region using its binding energy matrix from the BEEML database .', 'The starting position of the lowest energy binding site was set to be the summit of each region .', 'We then generated the start positions of fragments ( 200 bp in length ) in the ChIP sample by setting the fragment jitter to 50 bp .', 'The read length for the paired-end reads were set at 50 bp at a sequencing depth of 100 .', 'The sequence of each read was assigned based on the sequence of the S288C R64-2-1 reference of the S . cerevisiae genome [36] from the Saccharomyces Genome Database [37] .', 'The reads generated were aligned to this genome reference using the BWA aligner [38] .'], ['ChIPulate substantially and qualitatively extends previous statistical models of read counts obtained in a ChIP-seq experiment .', 'Early work by Zhang et .', 'al assumed a distribution of binding affinities across the genome and computed occupancy using a similar biophysical model to ours [24] , but did not simulate extraction and amplification .', 'Bao et .', 'al explicitly accounted for the effect of extraction efficiency on read counts [26] , but did not model the underlying biophysical occupancy of genomic loci based on their binding energies .', 'More recent work by Ruan et .', 'al , accounted for the effect of chemical potential on TF-DNA occupancy and its impact on motif inference [14] .', 'However , the authors evaluated the impact of experimental noise by adding Gaussian noise to the energies of binding sites rather than modeling the downstream processes of the ChIP-seq protocol .', 'Our pipeline combines a biophysical model of TF-DNA binding with a detailed simulation of fragment extraction , amplification and sequencing .', 'This allows us to analyze the impact of heterogeneity in fragment extraction and PCR amplification on the motif inferred from bound genomic loci , and the fidelity of read count ratios in discriminating relatively close binding energies .', 'We also evaluated how these factors , along with chromatin accessibility and sequencing depth , affect the false positive peak detection rate .', 'The role of biological factors such as indirect binding and cooperative binding , which are frequent occurrences in ChIP-seq datasets , were also evaluated for their impact on motif inference and fidelity .', 'Finally , we measured the accuracy of using read count ratios to infer the energy of a binding site and calculated the improvement in this accuracy after several biological replicates are performed .', 'In our simulations , fragment extraction encapsulates multiple steps in the ChIP-seq protocol , namely , cross-linking TF-DNA bound complexes , antibody-mediated pull down of fragments bound by the TF , and the removal of cross-links .', 'In addition to the mean extraction efficiency , the nature of the extraction heterogeneity dictates the magnitude of its impact on motif inference and fidelity , with power law distributions having a greater impact than normally distributed heterogeneity .', 'The effect of extraction heterogeneity on motif inference also depended on the structural class of the TF , where we found that a low mean extraction efficiency distorted motifs with a lower information content much more than those with high information content .
['Introduction', 'Results', 'Discussion', 'Methods']
['ChIP-seq ( Chromatin Immunoprecipitation followed by sequencing ) is a high-throughput technique to identify genomic regions that are bound in vivo by a particular protein , e . g . , a transcription factor ( TF ) .', 'Biological factors , such as chromatin state , indirect and cooperative binding , as well as experimental factors , such as antibody quality , cross-linking , and PCR biases , are known to affect the outcome of ChIP-seq experiments .', 'However , the relative impact of these factors on inferences made from ChIP-seq data is not entirely clear .', 'Here , via a detailed ChIP-seq simulation pipeline , ChIPulate , we assess the impact of various biological and experimental sources of variation on several outcomes of a ChIP-seq experiment , viz . , the recoverability of the TF binding motif , accuracy of TF-DNA binding detection , the sensitivity of inferred TF-DNA binding strength , and number of replicates needed to confidently infer binding strength .', 'We find that the TF motif can be recovered despite poor and non-uniform extraction and PCR amplification efficiencies .', 'The recovery of the motif is , however , affected to a larger extent by the fraction of sites that are either cooperatively or indirectly bound .', 'Importantly , our simulations reveal that the number of ChIP-seq replicates needed to accurately measure in vivo occupancy at high-affinity sites is larger than the recommended community standards .', 'Our results establish statistical limits on the accuracy of inferences of protein-DNA binding from ChIP-seq and suggest that increasing the mean extraction efficiency , rather than amplification efficiency , would better improve sensitivity .', 'The source code and instructions for running ChIPulate can be found at https://github . com/vishakad/chipulate .']
['DNA-binding proteins perform many key roles in biology , such as transcriptional regulation of gene expression and chromatin modification .', 'ChIP-seq ( Chromatin immunoprecipitation followed by high-throughput sequencing ) is a widely used experimental technique to identify DNA-binding sites of specific proteins of interest , within cells , genome-wide .', 'DNA fragments from genomic regions that are bound by a protein of interest , often a transcription factor ( TF ) , are selectively extracted using specific antibodies , amplified using PCR , and sequenced .', 'The sequences are mapped to the reference genome .', 'Regions where many sequences map , called “peaks” , are used to infer the location of TF-bound loci ( peaks ) , in vivo occupancy at those loci , and the sequence pattern ( motif ) to which the TF shows a binding affinity .', 'But measurements of TF occupancy and motif inference are vulnerable to several biological and experimental sources of variation that are poorly understood and difficult to assess directly .', 'Here , we simulate key steps of the ChIP-seq protocol with the aim of estimating the relative effects of various sources of variations on motif inference and binding affinity estimations .', 'Besides providing specific insights and recommendations , we provide a general framework to simulate sequence reads in a ChIP-seq experiment , which should considerably aid in the development of software aimed at analyzing ChIP-seq data .']
['simulation and modeling', 'probability distribution', 'mathematics', 'genome analysis', 'sequence motif analysis', 'molecular biology techniques', 'epigenetics', 'chromatin', 'research and analysis methods', 'sequence analysis', 'genomic libraries', 'genomics', 'artificial gene amplification and extension', 'bioinformatics', 'chromosome biology', 'gene expression', 'biological databases', 'molecular biology', 'probability theory', 'biochemistry', 'biochemical simulations', 'sequence databases', 'cell biology', 'normal distribution', 'database and informatics methods', 'genetics', 'biology and life sciences', 'physical sciences', 'computational biology', 'polymerase chain reaction']