# Problems about value log2(IP/Input) less than zero in ChIP-seq?

I have a question about the normalization for ChIP-seq. I used CPM to normalize my bam files of each IP and Input. Then I calculate the coverage of gene bodies for all genes on the genome. I have WT and different mutants. Take H3K4me1 as an example, I first calculate the matrix over gene bodies for H3K4me1 IP and Input. Then I used the matrix of H3K4me1 divided by the matrix of input. then log2(H3K4me1/Input) and plot this value over entire gene bodies. I found all the regions are lower than 0 and it seems that my IP does not work. In fact, this is not the case since I have checked the antibody's efficiency via ChIP-qPCR and it works well. Moreover, the distribution of H3K4me1 is consistent with the published data. But the problem is the profile of log2(H3K4me1/Input) is lower than 0. I don't know how to deal with such problems? I would appreciate it if someone can give me some suggestions on that. Thanks!

Seconding what Devon says:

1. Use a proper normalization technique such as TMM from edgeR. For applying this to bigwig files you could follow my tutorial at biostars.org.

My two cents beyond that:

1. I would not use enrichment over input as input is hard to normalize and can be sparse depending on sequencing depth. Therefore, you might get biased enrichments that do not well reflect reality, especially if counts in IP or input (of both) are low. I personally only use input for peak calling, and for plots like this then compare IPs among each other.

So the strategy would be to first normalize your bigwigs (I guess you use these for the plots?) using the strategy as linked above, and then use these robust CPMs to make the plots without any input samples, only using the IPs. You can use MA-plots to explore whether samples are properly normalized. I made a little code suggestion here in this biostars post on how to do this starting basically from a matrix of raw counts using edgeR. This is basically (from the TMM normalization standpoint) the same as in the first link, but then doing pairwise MA-plots to see whether TMM properly centers the majority of points along y-axis being somewhat zero. Normalization is key in any experiment, therefore I always recommend to use these kinds of plots for diagnostics.

• Dear ATpoint, sorry for replying to your message so late and thank you very much for your suggestions! Yes, you are right, I first want to normalize the raw bigwig files once I get the normalized bigwig files, I will try to calcualte the matrix of genebodies for genes interested in each sample. I will have a try to use TMM method you recommended to see whether it works on my data or not! I will let you know once I got the results. Thanks! – Lingling Yang Jun 17 at 11:43
• Hi, ATpoint, I read the link you sent to me. I have a question about the regions in the matrix . If I understand properly, for example, I have 5 strains and each strain have a peak region list, you means that pooling the peak coordinates of each sample as one region list, then use the feature count to calculate the coverage of each region for all samples, if I am wrong, please correct me. I am looking forward to your reply! Best. – Lingling Yang Jun 17 at 13:45

This is a common issue caused by the normalization method not being particularly robust (most likely it's using the mean signal for normalization, which will be a bit upward skewed when you have regions of high enrichment in the ChIP). The basic ways to fix this are:

• Use a more robust normalization (e.g., the median signal, or even just manually select un-enriched regions to calculate the normalization factors)
• Manually tweak the normalization factors so the background (i.e., log2 ratio in regions where you don't expect any enrichment) is 1. That will generally make the regions you're plotting more biologically realistic.

Note that there can be a negative enrichment in some areas, due to things like exclusion of histones, so don't be overly zealous if you manually tweak the normalization factors.

• Dear Devon, Thanks for your kind reply and valuable suggestions! I have tried to use the median value of the non-enriched region to calculate the scaled factor of each IP vs their input. I found it works for all modifications but failed for H3 IP. I guess this is due to the peaks for H3 are widely spread on the genome and when I use MACS2 to do the peak calling for H3, it won't work efficiently. Do you have any suggestions for H3 IP normalization? I am looking forward to your reply. Thanks! – Lingling Yang Jun 17 at 11:36
• You may have to manually tweak normalization values to get reasonable results. – Devon Ryan Jun 17 at 11:38
• What I did: first did peak annotation, then extracted the genes(name it as gene list 1) with enriched peaks. Then I calculated the matrix of the gene body of genes (not included in gene list 1 ), further calculate the median value for each sample (for both IP and input). In the peak annotation step, I tried Homer and Rgmatch, unfortunately, the results of these two software are quite different. – Lingling Yang Jun 17 at 12:16
• Therefore, I recalculated the median value given the coverage of the non-peak region (just escaped the peak annotation step). It seems the results of using Homer were much closer to the results directly given the scaled factor for the non-peaks region. I checked the default parameter of Homer, they defined TSS as from -1kb to +100bp and TES as -100bp to +1kb. Since S.cerevisiae. genome is not big like human or mouse, I am not sure how much effect using the default parameter. I am pretty interested in the 500bp up-and downstream of gene bodies. – Lingling Yang Jun 17 at 12:22
• The Homer source code is based on perl and I haven't used Perl before. Do you have any other suggestions? Thanks！ – Lingling Yang Jun 17 at 12:22