Peak callers for ChIP seq comparison

Question

[Wet&DryImmunology]

Hi, I am quite a novice for NGS analysis. I had a conversion with my colleagues about peak calling for ChIP seq data. I got confused about several comments she mentioned.

1. comment one "MACS is too old, and no one use it any more", is that true? It seems to me that it is still the most widely used, if not the best, peak caller now. I could still see it being used in the most recent high profile journals.

2. "Setting the parameter of peak calling is an art", OK, i made this up. She said it is essential to tailor parameters to fit each set of data. I understand that it is important to adjust the parameters based on whether the enrichment is broad (many histone modifications) or narrow (TFs). but she comments seems suggest that you could tailor you parameter as much as I want, as long as you apply the same parameter to the same dataset you are supposed to look at. Are there general rules or principles.

Another is that are there any peak caller optimized for low cell number ChIP seq, say 110^5 as input for H3K27ac ? I applied the MACS14 to two inputs using the same parameters, one has 10 times cell numbers ( 110^6 ) , the generated bed could be used for downstream analysis, but the low cell number inputs had some issue with downstream analysis.

Thank you for any suggestion and comments. (simple links to relavent literature would be appreciated)

[jotan]

I would agree with your colleague that peak calling is an "art". It's actually more like witchcraft.

MACS is quite old but as far as I can tell, none of the newer peak callers are much better. I use the peak callers to start, then I filter through a human eyeball attached to a brain, and I use the lab to verify. Select a range of peaks and ChIP-qPCR until my replicates start failing. That's my personal workflow.

The ENCODE guidelines are a good place to start.

Cell numbers requirements will vary between different types of samples. For example, you can get away with far fewer cells for TFs (point source) compared to histone mods (broad source).

[John]

I would avoid peak callers at every opportunity. Their very purpose in life doesn't make sense to me. So you want to take 2 dimensional data and make it 1 dimension? Ok. Why?

[dariober]

John (I assume your question is rhetorical), are you suggesting using the continuous signal throughout the genome as input for further analysis, like differential binding, motif analysis, etc? In principle it's not bad idea I think, but it's just impractical for most purposes.

Personally I tend to see the peak calling exercise as a first step to thrash 99.9% of the genome as "not interesting" and the rest (which is still a lot of stuff) as "worth a further look".

Peaks are still 2D as they have can be quantified as enrichment over flanking regions or input, isn't it?

[John]

Unfortunately, this time i'm not joking around. Their existence is, in my opinion (and i'm no one special, I don't even have a doctorate), a poor way of dealing with an otherwise hard problem. I feel that they exist only to reduce the complexity of the piled up signal to a list of "interesting regions", for no reason other than it makes working with the data more manageable. There is no scientific or mathematical utility in doing this, other than the existing tools we have at our disposal take regions of genome as input, not a BigWig/Bedgraph (with the assay dynamics encoded in them somehow).

When I first used peak callers 3 years ago, I thought they were silly. A quirk of a new line of research that hasn't settled in yet. I was quite verbose about it on Biostars and elsewhere too, but over time I realised that:

1. I wasn't going to change anyone's opinion any time soon. The alternative meant writing new analysis programs, and without an alternative in place my criticism was basically just me whining.
2. Maybe the pros in the time saved in doing comparisons quickly and simply would out-weight the cons of reducing the richness of the data.

I was wrong about point 2. Peak callers caused researchers to think and reason about ChIP-Seq as a collection of "peaks", which muddled many people's thinking. Are they "broad" peaks or "narrow" peaks? How much should two peaks have to overlap before they're considered overlapping? Should that cut-off be in bases or as a % of the peak width? Should I use peak confidence scores or sum of signal in the peak regions? And the list goes on and on and on..

But i'd be OK with it if people knew what they were doing and where happy with the process - but no, due to all the artificial complexities thrown up by that unnecessary abstraction, we now need an abstraction for our abstraction. Chromatin State callers. We don't even look at regions any more, we look at colours -- often more colours than we had dimensions of data to begin with.

[ivivek_ngs]

I would second that, peak calling might not be represented as envisioned but it will definitely help us pointing the 5% or 10% interesting regions in the genome where you have binding of important TFs/histones that might give an idea of the downstream biological effects. Obviously motif enrichment or differential binding is one way to pursue and avoid peak calling but the motivation lies in the user, if one wants to see for enriched regions around the TSS and that might be having motifs then stricter peaks can be used to model that, it is just a way of quantification to remove non important regions. At the end of the day it will be dependent on the biological question and the resources that can be used and exploited to either find enriched TFs or histone modifications across regions of interest in the genome that might drive the phenotype and then integrate it to expression data, in fact helps us to identify the impact of this regions on expression of your subsets of samples that have both ChIP-Seq and RNASeq. There are other ways of also finding promoters and enhancers from ChIP-Seq data rather than distance-based metrics without calling peaks but then it is something that will be sole interest of the study concerned.

Answer 1

[Ming Tommy Tang]

1. MACS14 is still good and are widely used.
2. yes, one should tune the parameters for peak calling. narrow vs wide

I put some of my notes here https://github.com/crazyhottommy/ChIP-seq-analysis

I think 10^5 cells are still a lot of cells although usually it is 5 million cells per ChIP. one should still get very good enrichment with 10^5 cells.

[dariober]

Hi- I agree with this answer. Just a couple of comments: MACS14 is pretty much replaced by MACS2, I would use that instead.

Peak calling is frustrating. Sometimes you see peaks that look identical, one is identified by macs (or else) with high confidence, the other is ignored. Some called peaks do not look like peaks at all.

I don't think peak callers are designed to work with high or low cell number in mind. What they "see" is the signal to noise, which is typically higher if you work more cells.

[Ming Tommy Tang]

MACS2 is not published yet? I think MACS14 is very stable and very good when calling point source TFs and histone marks (H3K4me3 and H3K27ac). I tried MACS2 on H3K27ac and the results are not as good as MACS14. Of course, it depends on tuning the parameters. Just my personal experiences on H3K27ac, I found MACS14 is good at at a Pvlaue cut-off of 10-e-8. MACS2 if use a pvalue of 0.01, too many peaks, if a qvalue of 0.01, too few peaks.

[ivivek_ngs]

even though MACS2 is not published it has been well used by a lot of publications, MACS14 did not have broad peak calling. It is there in the new MACS2, there is also a change in the default parameter set up in MACS2 from MACS14 and that is the reason it does not work with pvalue 0.01 , it infact works for q value, when you put p value irrespective of significance all peaks are given. So it is better to run MACS2 with q value and also if you run broad use both q value and broadcutoff , these things are not very clearly mentioned in the documentation. I have tested various scenarios between the different version. The latest MACS2.1 is quite stable. Also the number peaks is not a definitive goal for the enrichment. I would rather go what downstream inference my peak gives be it ontology enrichment or motifs extracted from the peaks. One should also be carefully try to see to what extent the strand cross correlation is effective and to what extent the experiment is good enough to go for peak calling. That can be done with SPP or CHASE and these infact defines also the number of peaks that we get post peak calling. MACS2 also needs the fragment size and most people put a fragment size of 146 for humans depending on ENCODE which is not true, it is important to determine the fragment size of the experiments and then use it for better sensitive peak calling.

[Ming Tommy Tang]

Thanks for sharing. what's the difference between q value and broadcutoff when calling broad peaks with MACS2? I do use phantomPeakqual to get the fragment size for MACS2.

[ivivek_ngs]

If you look at this link , there is a very nice issue that was brought up using macs2.1, the p value and q value corresponds to narrow peaks inside broad and once you have the broad peak there is another level of peak calling to nearby regions and then assigned a significance, so I would believe it is a combinatorial approach. This is the reason I always felt the documentation is not that rich and I started playing with different peak callers specially for histone modifications.

Answer 2

[ivivek_ngs]

For narrow peak calls already MACS2.1 is there which was updated last year, you can take a look here . You have to run them for narrow and broad peaks separately choosing the parameter of fragment size and q-value thresholds. I would widely say to use MACS2.1, it is coming out as better for histone modifications , but still if you are having doubt , now it is the best time to play with other methods like SICER (link) and PePr, I am kinda fond of PePr (link) as it empirically estimates the fragment length/shift size and sliding window width which usually is a parameter to be supplied for others and estimates dispersion over local genomic area. and uses it for downstream analysis of finding peak considering read counts across replicates and between groups with a negative binomial distribution. MACS2.1 for histone is usually used with -nomodel paramter. So you can play around and check. The blog mentioned above is amazing and he did a great job. You can play with your data using the various parameters. Everything relies at the end on the S/N ratio and the number of reads mapped on your genome of interest and based on these you make the peak calling. So proper QC of the inputs and treatment samples should be done at level of experimental quality and also post alignment prior to peak calling. This will help in mitigating the noise and use fair parameters for peak calling thus giving much high confidence peaks.