MACS2 Failing to produce a good amount peaks, causing failure at stage 5

[RJPNewell]

Hi!

I hope you are doing well, I just have run into an issue while using the mango pipeline.

I'm getting a failure at stage 5 of the mango pipeline with the following error:

[1] "counting reads per peak"
[1] "determining self-ligation distance"
[1] "self-ligation cutoff = 4498.4"
[1] "grouping PETs into interactions"
[1] "---Splitting PETs by chromosome"
[1] "---Intersecting PETs with peaks"
[1] "---Building putative interaction set"
[1] "modeling PETs based on peak depth and distance"
Error in splitdata[[i]] : subscript out of bounds
Calls: binmaker
Execution halted

This seems to be caused by the fact that MACS2 is only able to call a grand total of 30 peaks from the tagAlign file, which is alarming. The following is the output from MACS2:

INFO @ Sun, 04 Mar 2018 14:39:43: #1 read tag files...
INFO @ Sun, 04 Mar 2018 14:39:43: #1 read treatment tags...
INFO @ Sun, 04 Mar 2018 14:39:44: 1000000
INFO @ Sun, 04 Mar 2018 14:39:45: #1 tag size is determined as 20 bps
INFO @ Sun, 04 Mar 2018 14:39:45: #1 tag size = 20
INFO @ Sun, 04 Mar 2018 14:39:45: #1 total tags in treatment: 1724283
INFO @ Sun, 04 Mar 2018 14:39:45: #1 user defined the maximum tags...
INFO @ Sun, 04 Mar 2018 14:39:45: #1 filter out redundant tags at the same location and the same strand by allowing at most 1 tag(s)
INFO @ Sun, 04 Mar 2018 14:39:45: #1 tags after filtering in treatment: 1233394
INFO @ Sun, 04 Mar 2018 14:39:45: #1 Redundant rate of treatment: 0.28
INFO @ Sun, 04 Mar 2018 14:39:45: #1 finished!
INFO @ Sun, 04 Mar 2018 14:39:45: #2 Build Peak Model...
INFO @ Sun, 04 Mar 2018 14:39:45: #2 looking for paired plus/minus strand peaks...
INFO @ Sun, 04 Mar 2018 14:39:48: #2 number of paired peaks: 55103
INFO @ Sun, 04 Mar 2018 14:39:48: start model_add_line...
INFO @ Sun, 04 Mar 2018 14:39:48: start X-correlation...
INFO @ Sun, 04 Mar 2018 14:39:48: end of X-cor
INFO @ Sun, 04 Mar 2018 14:39:48: #2 finished!
INFO @ Sun, 04 Mar 2018 14:39:48: #2 predicted fragment length is 298 bps
INFO @ Sun, 04 Mar 2018 14:39:48: #2 alternative fragment length(s) may be 298 bps
INFO @ Sun, 04 Mar 2018 14:39:48: #2.2 Generate R script for model : L1/mm10_H3K4me1_cerebellum_L001_AC_model.r
INFO @ Sun, 04 Mar 2018 14:39:48: #3 Call peaks...
INFO @ Sun, 04 Mar 2018 14:39:48: #3 Pre-compute pvalue-qvalue table...
INFO @ Sun, 04 Mar 2018 14:39:50: #3 Call peaks for each chromosome...
INFO @ Sun, 04 Mar 2018 14:39:51: #4 Write output xls file... L1/mm10_H3K4me1_cerebellum_L001_AC_peaks.xls
INFO @ Sun, 04 Mar 2018 14:39:51: #4 Write peak in narrowPeak format file... L1/mm10_H3K4me1_cerebellum_L001_AC_peaks.narrowPeak
INFO @ Sun, 04 Mar 2018 14:39:51: #4 Write summits bed file... L1/mm10_H3K4me1_cerebellum_L001_AC_summits.bed
INFO @ Sun, 04 Mar 2018 14:39:51: Done!

It's that narrowPeak file that only contains 30 or so peaks. Now, I thought this would be due to a low amount of uniquely mapped reads during the alignment step but that doesn't seem to be the case. I'm currently getting ~40 million uniques mapped reads. However, during the construction of the tagAlign file it appears that a lot of duplicate BEDPE entries are removed. This results in only 1.8 million tags for MACS2 to use.

Is this an issue that you have come across before, and would you potentially have any workarounds for an issue like this?

Thanks for your time!
Cheers,

Rhys