KeyError in MACS2

Question

[Kasthuri]

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'

Answer 1

[Devon Ryan]

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]

[Kasthuri]

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]

[Devon Ryan]

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.

[Kasthuri]

Done! Added this comment in GitHub and closed the ticket.

Answer 2

[James Ashmore]

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.

[Kasthuri]

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.