Question: KeyError in MACS2
0
gravatar for Kasthuri
3.6 years ago by
Kasthuri260
United States
Kasthuri260 wrote:

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'

 

chip-seq macs2 • 1.6k views
ADD COMMENTlink modified 3.6 years ago by Devon Ryan91k • written 3.6 years ago by Kasthuri260
2
gravatar for Devon Ryan
3.6 years ago by
Devon Ryan91k
Freiburg, Germany
Devon Ryan91k wrote:

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]
ADD COMMENTlink written 3.6 years ago by Devon Ryan91k
1

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]
ADD REPLYlink modified 3.6 years ago • written 3.6 years ago by Kasthuri260
1

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.

ADD REPLYlink written 3.6 years ago by Devon Ryan91k
1

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

ADD REPLYlink written 3.6 years ago by Kasthuri260
0
gravatar for James Ashmore
3.6 years ago by
James Ashmore2.6k
UK/Edinburgh/MRC Centre for Regenerative Medicine
James Ashmore2.6k wrote:

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.

ADD COMMENTlink written 3.6 years ago by James Ashmore2.6k

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.

ADD REPLYlink written 3.6 years ago by Kasthuri260
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1135 users visited in the last hour