That's great. Thanks Ryan! On 09/08/2014 04:40, Ryan wrote: > Hi Aliaksei, > > First of all, dividing by library size is not the way to go. It's wrong > in ChIP-seq for the same reason it's wrong in RNA-seq: compositional > bias. In order to properly normalize for sequencing depth, I would > recommend the strategy described here: > http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4066778/ > > You just need to look at the small section about calculating > normalization factors using 10kb bins. Once you calculate the > normalization factors, multiply them by the nominal library sizes to get > the normalized library sizes. Then divide the raw coverage by the > normalized library sizes to get the normalized coverage, which has been > adjusted for both sequencing depth and compositional bias. You should be > able to normalize input and ChIP samples together in this way as long as > your ChIP peaks are not omnipresent throughout the entire genome. > > From there, you could theoretically subtract from each sample the > corresponding input (or maybe the average of all inputs, since they > should all be equivalent) to get enrichment above input, or divide to > get fold enrichment over input. You are correct in your supposition that > arithmetic operations on Rle objects (i.e. the type of object returned > by coverage) work element-wise, just like arithmetic on ordinary R > numeric vectors. > > Of course you are correct that if the input coverage is zero at any > point, then dividing by it will result in an infinite ratio, which you > probably don't want. You could take another page out of the RNA-seq > playbook and add a small constant to each coverage object before > dividing to regularize the ratios. You could also smooth out the input > coverage by replacing each coverage value with a mean over a window > centered at that position. You can do this with the runmean function > from the caTools package (from CRAN, not Bioconductor), using options > endrule="mean" and align="center". Depending on the sparseness of your > coverage, sufficient smoothing might eliminate zero values. > > Anyway, that's about all I have to suggest. I hope you find something > useful in there. > > -Ryan > > On Fri Aug 8 10:19:39 2014, Aliaksei Holik wrote: >> Dear Bioconductors and especially *Ranges team, >> >> I would like your expert opinion on whether what I'm doing is a >> sensible thing to do. >> >> I wish to plot some ChIP read coverage data in Gviz. The way I do it >> is by reading a number of bam files with readGAlignments() and >> calculating the coverage with the likewise named function. I would >> also like to normalise the coverage for each bam file for library size >> and the Input sample. I do so by simply dividing the coverage objects >> by the respective library size and by the coverage object for the >> Input sample. It all works beautifully, and the graphs look ok. But >> I'm left with a nagging feeling, that it's too easy to be true. So I >> guess my question is: am I right to assume that by dividing one >> coverage object by another, the coverage at each position in one >> object would be divided by the coverage at the same very position in >> the other object. I'm almost expecting something to go wrong, if, for >> instance, I have 0 coverage at some position in my reference object. >> But perhaps I should have faith in Bioconductor. >> >> Many thanks, >> >> Aliaksei. >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: >> http://news.gmane.org/gmane.science.biology.informatics.conductor > >