CUT&Tag data processing: peak caller choice and typical peak size profile

Question

[nihilior]

Greetings everyone,

I am currently analysing CUT&Tag data (a recent technique, variant of the ChIP-seq), capturing different histone variants. To do so, I tried all the different suggested tools for peaks calling: the more general MACS2 (broad peaks, suggested for histones) and the more method-specific SEACR, gopeaks, and Cut&RunTools (which is actually a pipeline evoking MACS2/SEACR and that adds some filters). The problem is that I obtain drammaticaly different peak profiles from each tool, regarding both peak size and peak number, with MACS2 being the most generous and SEACR/gopeaks the most stringent. How can I know which is the typical peak size profile for a certain histone and/or is there a general consensus on which the expected peak sizes should be and some common practices to filter the obtained peaks (considering that replicates handling was already done with IDR) because the downstream results are heavly influenced by such choice. Thank you in advance,

[ATpoint]

Which marks did you pull down? Broad peaks only makes sense for broad marks. ENCODE has a table that is helpful under "Target-specific Standards" at https://www.encodeproject.org/chip-seq/histone/ listing typical broad and narrow marks. Can you share a screenshot from the IGV to see whether data are good quality?

[nihilior]

Hi ATpoint, thank you for your answer! We have macroH2A data, which is not strictly specified in the reported table (while I see the H2A variants being specified as narrow, I am not sure it is the same, since I see a lot of different H3 variants being broad or narrow). What kind of IGV screenshot can I share?

[ATpoint]

I usually make a bedGraph with bedtools genomecov and then write this to a bigwig with bedGraphToBigWig from kentUtils to inspect data on the IGV.

Answer 1

[jared.andrews07]

From experience, SEACR kind of sucks and way overcalls peaks if you don't have absolutely sterling quality data. Peak sizes are all over the place for histone marks, but I'd definitely consider macroH2A a "broad" peak.

IDR and/or consensus peaks if you have replicates are pretty common approaches. We may be able to provide more feedback if you lay out the experimental design. As ATpoint said, check the data quality first.

[ATpoint]

I personally found that macs2 is sufficient for CUT&RUN (TAG probably the same) but the data I produced was not as "clean" as they always advertise in their methods papers. It essentially looked like a normal ChIP-seq with the usual noise. I think you only need specialized tools if data are pristine and there is essentially only signal and no noise. In that case macs2 might produce too low pvalues as its model relies on presence of some noise to get its stats going iirc.

[jared.andrews07]

We found much the same. Mostly we end up hand tweaking/filtering MACS2 results.

Similar to ChIP, it seems to just really depend on the antibody in our hands. Most histone marks are fairly simple, but a lot of TF antibodies are just not very good/specific.

[nihilior]

Hi Jared, thank you very much for your answer. Here:

I share a screenshot of the data we have from IGV (two samples, with two replicates for the first and three for the second). We already know that the data is not sterling quality, as they are data from a CUT&Tag experiment on iPSC cells, and we received the data in which replicates are not that consistent. Some suggestions about the MACS2 tweaking/filtering process? Also, Jared, may we have a reference for the "broad" nature of macroH2A just to know? Thanks a lot again.

[jared.andrews07]

I just looked at some published data to determine they were clearly "broad" domains. Figures 4 & 5 in this paper show it pretty clearly.

In terms of tweaks to MACS2, it's mostly just playing with the q-value threshold and sometimes post hoc filtering based on the scores that are output. This vignette has a nice explanation of how adjusting q-value threshold affects not only whether a peak is called but also the width of peaks.

[nihilior]

thank you very much jared!

[nihilior]

We called peaks with MACS2 with the default q-value cutoff (0.1). Than we filtered the obtained peaks with a p-value threshold of 0.05. Our median peak size was around 30k bp, with a max value of 150k. Do you think that this peak sizes values are reasonable? We had some debate on wheather our histone peak size range profiles were compatible with typical macroH2A histone binding. We also don't know what the MACS2 score value exactly are, since the docs says: The score is for 'display' purpose only. It's for UCSC browser and will make stronger peaks darker, and weaker peaks lighter. May this score be used as another filter for our peaks? Sorry for all these questions.

[jared.andrews07]

No worries on the questions - it's good that you're actually thinking about this. I am not surprised by those peak sizes, but I consider peak boundaries pretty fickle beasts at the best of times.

You're right that you shouldn't use the scores for anything important, I forgot how they were derived. I'd look in IGV at some of your peaks and see if you believe them. If the boundaries seem off or you have lots of spurious calls, adjust the q-value cutoff and take a look again.

[nihilior]

Okay jared, thank you so much again for your time and help :)

[ATpoint]

You have to zoom in to individual loci, like individual genes, a 50kb window, something like that. Looking at it at that resolution you show allows no judgement of quality.

[nihilior]

[gglim]

Hi nihilior,

Did you figure it out? Recently I get bw files from CUT&TAG data (generated using deepTools bamcoverage) look like yours, and I don't know if it's okay to proceed.

The H3K4me3 and RNA polymerase II have beautiful signal tracks, but FLAG and control group look strange (auto-scaled). I've checked the alignment info, which shows that they all aligned well (>95%) but possess different sequencing depth. Could it be the reason and why?

[ATpoint]

The FLAG looks like noise, at least in this genomic region, the pol and histone mark is nice. Could you call any significant peaks from the FLAG?

[gglim]

Hi, ATpoint

Thanks for the reply! I got ~15000 peaks using default MACS2 parameter (with negative control as control group), but they show quite high q-values, mostly > 0.1. Does this indicate low sample quality ?

[ATpoint]

The default macs2 q-value cutoff is 0.05 -- so if you have many with higher q then indeed that fits the browser track that this particular IP did not work well.

[gglim]

Thanks a lot...It's strange that we obtained significant result from ChIP-qPCR, but ChIP-seq and CUT&Tag show nothing. This problem has been reported by others and remains unsolved. Do you have any idea?

[Alexandra]

Are there any ways to filter out the SEACR peaks in Cut&Run experiment? In our case they are too wide and, compared to bigwig, sometimes seem to consist of several peaks. The IgG data contained a lot of noise.

[ATpoint]

Have you tried MACS2 as suggested here?

The IgG data contained a lot of noise

Of course, it's noise by definition.

[Alexandra]

Yes, I tried both, but macs2 gave too many peaks (some of them are definitely fake). So we thought about filtering seacr somehow, because by eye they are sometimes more reasonable. But there is no p-value or any equivalent. About IgG... as I understand, the background should be lower than in Chip-seq, but for our data it is not.

[ATpoint]

I cannot help without some material to illustrate plus code examples. As for CUTnRUN vs ChIP-seq, in our hands, and what I hear from others for example jared.andrews07, CnR does not guarantee to have less noise. It can just look as crappy as regular ChIP-seq. Just because the original paper sells it as the new gold standard does not mean it is noise free.

[jared.andrews07]

It also has disturbingly different profiles than ChIP-seq for the same TFs at times, which we aren't real sure how to interpret. Even profiles for some histone marks, e.g. H3K27me3 look quite different.

macs2 gave too many peaks (some of them are definitely fake)

You can filter by p-value to clean some of them up (or adjust the peak calling parameters to mitigate the erroneous calls). Peak calling can be a finicky process, and there's not really a one size fits all solution that's going to work seamlessly.

[ATpoint]

What I always do with macs is to a) require a minimum peak width, for TFs for example 150bp. b) Use a control if you have it (in your case IgG). Be stringent in the filtering towards significance. Standard padj < 0.05 will generate too many peaks. c) If you have replicates use IDR or similar setups that scan for consistent peaks. In my hands one can easily get thousands of peaks more or less depending on how you set cutoffs and choice of (IDR/consistency tool). In the end, I often just look at the peaks for like 20 loci randomly, and then ask myself if I could confidently stand up for these results in public. If not, keep filtering. IP experiments have considerable noise, so it is really critical to look at the data and apply filters as needed.

[jared.andrews07]

Yup, the eyeball test still reigns supreme.