I get the following error when using MACS2 on a histone mark in ChIP-seq. I saw online that this error is seen with few other people but I couldn't find any resolution. I have posted this in their Github forum as well. Any help will be much appreciated. Thanks.
This is the command I issue,
/Users/home/Tools/MACS-master/bin/macs2 callpeak --format BED --nomodel --broad --verbose 2 --qvalue 0.25 --g mm --name Input.vs.Exp --treatment treatment.bg --control control.bg --outdir .
and I get this error,
INFO @ Fri, 15 Jan 2016 21:12:49:
# Command line: callpeak --format BED --nomodel --broad --verbose 2 --qvalue 0.25 --g mm --name Input.vs.Exp --treatment treatment.bg --control control.bg --outdir .
# ARGUMENTS LIST:
# name = Input.vs.Exp
# format = BED
# ChIP-seq file = ['treatment.bg']
# control file = ['control.bg']
# effective genome size = 1.87e+09
# band width = 300
# model fold = [5, 50]
# qvalue cutoff for narrow/strong regions = 2.50e-01
# qvalue cutoff for broad/weak regions = 1.00e-01
# Larger dataset will be scaled towards smaller dataset.
# Range for calculating regional lambda is: 1000 bps and 10000 bps
# Broad region calling is on
INFO @ Fri, 15 Jan 2016 21:12:49: #1 read tag files...
INFO @ Fri, 15 Jan 2016 21:12:49: #1 read treatment tags...
INFO @ Fri, 15 Jan 2016 21:12:50: 1000000
INFO @ Fri, 15 Jan 2016 21:12:51: 2000000
INFO @ Fri, 15 Jan 2016 21:12:53: 3000000
INFO @ Fri, 15 Jan 2016 21:12:54: 4000000
INFO @ Fri, 15 Jan 2016 21:12:55: 5000000
INFO @ Fri, 15 Jan 2016 21:12:56: 6000000
INFO @ Fri, 15 Jan 2016 21:12:57: #1.2 read input tags...
INFO @ Fri, 15 Jan 2016 21:12:58: 1000000
INFO @ Fri, 15 Jan 2016 21:13:00: 2000000
INFO @ Fri, 15 Jan 2016 21:13:01: 3000000
INFO @ Fri, 15 Jan 2016 21:13:02: 4000000
INFO @ Fri, 15 Jan 2016 21:13:03: #1 tag size is determined as 41 bps
INFO @ Fri, 15 Jan 2016 21:13:03: #1 tag size = 41
INFO @ Fri, 15 Jan 2016 21:13:03: #1 total tags in treatment: 6600699
INFO @ Fri, 15 Jan 2016 21:13:03: #1 total tags in control: 4295400
INFO @ Fri, 15 Jan 2016 21:13:03: #1 finished!
INFO @ Fri, 15 Jan 2016 21:13:03: #2 Build Peak Model...
INFO @ Fri, 15 Jan 2016 21:13:03: #2 Skipped...
INFO @ Fri, 15 Jan 2016 21:13:03: #2 Use 200 as fragment length
INFO @ Fri, 15 Jan 2016 21:13:03: #3 Call peaks...
INFO @ Fri, 15 Jan 2016 21:13:03: #3 Call broad peaks with given level1 -log10qvalue cutoff and level2: 0.602060, 1.000000...
INFO @ Fri, 15 Jan 2016 21:13:03: #3 Pre-compute pvalue-qvalue table...
INFO @ Fri, 15 Jan 2016 21:13:22: #3 Call peaks for each chromosome...
Traceback (most recent call last):
File "/Users/kasthuri/Tools/MACS-master/bin/macs2", line 617, in
main()
File "/Users/kasthuri/Tools/MACS-master/bin/macs2", line 57, in main
run( args )
File "/Library/Python/2.7/site-packages/MACS2-2.1.0.20151222-py2.7-macosx-10.9-intel.egg/MACS2/callpeak_cmd.py", line 264, in run
peakdetect.call_peaks()
File "PeakDetect.pyx", line 105, in MACS2.PeakDetect.PeakDetect.call_peaks (MACS2/PeakDetect.c:1650)
File "PeakDetect.pyx", line 259, in MACS2.PeakDetect.PeakDetect.__call_peaks_w_control (MACS2/PeakDetect.c:3248)
File "CallPeakUnit.pyx", line 1428, in MACS2.IO.CallPeakUnit.CallerFromAlignments.call_broadpeaks (MACS2/IO/CallPeakUnit.c:18243)
File "CallPeakUnit.pyx", line 1501, in MACS2.IO.CallPeakUnit.CallerFromAlignments.call_broadpeaks (MACS2/IO/CallPeakUnit.c:17770)
File "PeakIO.pyx", line 158, in MACS2.IO.PeakIO.PeakIO.get_data_from_chrom (MACS2/IO/PeakIO.c:3309)
KeyError: 'chr12'
It looks like that will happen whenever there are no called peaks on a chromosome. The only way around that would be to remove that chromosome from the analysis (or fix the macs2 code, by making this function this:
def get_data_from_chrom (self, str chrom):
if not self.peaks.has_key:
return []
return self.peaks[chrom]
Awesome, I changed the macs2 code and it works after compiling. However, there is a slight change in the syntax. With the current MACS-master.zip checkout, it should be,
def get_data_from_chrom (self, str chrom):
if not self.peaks.has_key(chrom):
self.peaks[chrom]=[]
return self.peaks[chrom]
Thanks for finding my bug, I should have seen that :) If you haven't already done so, please add this to your issue on github so it can get fixed.
It looks like your treatment and control files are in bedgraph format. MACS2 can call peaks from this format but you have to use the bdgpeakcall command.
•
written
by
James Ashmore
Unfortunately, bdgpeakcall doesn't have the option to call with the control file. In fact, I had this error even with bam file with --format BAM option. Also, it ran perfectly fine with other bedgraph files that I use with the same exact command.
Log in to answer this question.
