Q: I got an Error with PeakDetect.call_peaks while running callpeak with pair-end MEDIP seq data

[chantisakee]

my data is : Illumina paired-end MEDIP seq

I used -f BAMPE -g hs -p 1.00e-02
and I got this error. please give me any suggestions

---

effective genome size = 2.70e+09
band width = 300
model fold = [5, 50]
pvalue cutoff = 1.00e-02
qvalue will not be calculated and reported as -1 in the final output.
The maximum gap between significant sites is assigned as the read length/tag size.
The minimum length of peaks is assigned as the predicted fragment length "d".
Larger dataset will be scaled towards smaller dataset.
Range for calculating regional lambda is: 1000 bps and 10000 bps
Broad region calling is off
Paired-End mode is on

INFO @ Tue, 20 Apr 2021 15:21:07: #1 read fragment files...
INFO @ Tue, 20 Apr 2021 15:21:07: #1 read treatment fragments...
INFO @ Tue, 20 Apr 2021 15:21:09: 1000000
INFO @ Tue, 20 Apr 2021 15:21:12: 2000000
INFO @ Tue, 20 Apr 2021 15:21:14: 3000000
INFO @ Tue, 20 Apr 2021 15:21:17: 4000000
INFO @ Tue, 20 Apr 2021 15:21:19: 5000000
INFO @ Tue, 20 Apr 2021 15:21:21: 6000000
INFO @ Tue, 20 Apr 2021 15:21:24: 7000000
INFO @ Tue, 20 Apr 2021 15:21:26: 8000000
INFO @ Tue, 20 Apr 2021 15:21:28: 9000000
INFO @ Tue, 20 Apr 2021 15:21:30: 10000000
INFO @ Tue, 20 Apr 2021 15:21:33: 11000000
INFO @ Tue, 20 Apr 2021 15:21:47: 11517512 fragments have been read.
INFO @ Tue, 20 Apr 2021 15:21:55: #1.2 read input fragments...
INFO @ Tue, 20 Apr 2021 15:21:57: 1000000
INFO @ Tue, 20 Apr 2021 15:21:59: 2000000
INFO @ Tue, 20 Apr 2021 15:22:02: 3000000
INFO @ Tue, 20 Apr 2021 15:22:04: 4000000
INFO @ Tue, 20 Apr 2021 15:22:06: 5000000
INFO @ Tue, 20 Apr 2021 15:22:09: 6000000
INFO @ Tue, 20 Apr 2021 15:22:11: 7000000
INFO @ Tue, 20 Apr 2021 15:22:13: 8000000
INFO @ Tue, 20 Apr 2021 15:22:16: 9000000
INFO @ Tue, 20 Apr 2021 15:22:18: 10000000
INFO @ Tue, 20 Apr 2021 15:22:20: 11000000
INFO @ Tue, 20 Apr 2021 15:22:23: 12000000
INFO @ Tue, 20 Apr 2021 15:22:25: 13000000
INFO @ Tue, 20 Apr 2021 15:22:27: 14000000
INFO @ Tue, 20 Apr 2021 15:22:39: 14942679 fragments have been read.
INFO @ Tue, 20 Apr 2021 15:22:49: #1 mean fragment size is determined as 0.0 bp from treatment
INFO @ Tue, 20 Apr 2021 15:22:49: #1 note: mean fragment size in control is 0.0 bp -- value ignored
INFO @ Tue, 20 Apr 2021 15:22:49: #1 fragment size = 0.0
INFO @ Tue, 20 Apr 2021 15:22:49: #1 total fragments in treatment: 11517512
INFO @ Tue, 20 Apr 2021 15:22:49: #1 user defined the maximum fragments...
INFO @ Tue, 20 Apr 2021 15:22:49: #1 filter out redundant fragments by allowing at most 1 identical fragment(s)
INFO @ Tue, 20 Apr 2021 15:23:03: #1 fragments after filtering in treatment: 93
INFO @ Tue, 20 Apr 2021 15:23:03: #1 Redundant rate of treatment: 1.00
INFO @ Tue, 20 Apr 2021 15:23:03: #1 total fragments in control: 14942679
INFO @ Tue, 20 Apr 2021 15:23:03: #1 user defined the maximum fragments...
INFO @ Tue, 20 Apr 2021 15:23:03: #1 filter out redundant fragments by allowing at most 1 identical fragment(s)
INFO @ Tue, 20 Apr 2021 15:23:20: #1 fragments after filtering in control: 93
INFO @ Tue, 20 Apr 2021 15:23:20: #1 Redundant rate of control: 1.00
INFO @ Tue, 20 Apr 2021 15:23:20: #1 finished!
INFO @ Tue, 20 Apr 2021 15:23:20: #2 Build Peak Model...
INFO @ Tue, 20 Apr 2021 15:23:20: #2 Skipped...
INFO @ Tue, 20 Apr 2021 14:51:02: #3 Call peaks...
Traceback (most recent call last):
File "/ckeerati/.conda/envs/MACS2-CK/bin/macs2", line 653, in
main()
File "/ckeerati/.conda/envs/MACS2-CK/bin/macs2", line 51, in main
run( args )
File "/ckeerati/.conda/envs/MACS2-CK/lib/python3.6/site-packages/MACS2/callpeak_cmd.py", line 256, in run
peakdetect.call_peaks()
File "MACS2/PeakDetect.pyx", line 107, in MACS2.PeakDetect.PeakDetect.call_peaks
File "MACS2/PeakDetect.pyx", line 155, in MACS2.PeakDetect.PeakDetect.__call_peaks_w_control
ZeroDivisionError: float division

---

[chantisakee]

here is my samples' flagstat (hope it might help)

(Treatment sample)
19978542 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
11827028 + 0 mapped (59.20% : N/A)
7921795 + 0 paired in sequencing
3960070 + 0 read1
3961725 + 0 read2
0 + 0 properly paired (0.00% : N/A)
96766 + 0 with itself and mate mapped
212750 + 0 singletons (2.69% : N/A)
75534 + 0 with mate mapped to a different chr
46152 + 0 with mate mapped to a different chr (mapQ>=5)

(Control sample)
21281481 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
15450910 + 0 mapped (72.60% : N/A)
5734130 + 0 paired in sequencing
2855852 + 0 read1
2878278 + 0 read2
0 + 0 properly paired (0.00% : N/A)
210608 + 0 with itself and mate mapped
297623 + 0 singletons (5.19% : N/A)
185908 + 0 with mate mapped to a different chr
143984 + 0 with mate mapped to a different chr (mapQ>=5)

Thank you in advance :)

[taoliu]

@chantisakee As for reading paired-end BAM files, MACS will only keep those 'properly paried'. However, it seems that there are 0 properly paired reads....

[chantisakee]

Thanks a lot for your response.
According to your answer I've already change the alignment parameters and it works fine now!

Thank you :]

[RocesV]

Hello everyone! Actually facing the same problem. How did you change the alignment parameters in order to solve it?

Thanks in advance :)