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?