Question: Macs2 callpeak error: IndexError: list index out of range
0
gravatar for Gary
22 months ago by
Gary450
Taiwan/Taichung/China Medical University Hospital
Gary450 wrote:

Hi,

I run macs2 for peak calling and get an error message below. I know there are similar questions on line, but I couldn't find a resolution for me. The below is my detail information. Could you help me to deal with this problem? Many thanks.

I. I run Bowtie2 on Partek Flow to align 40 bp paired-end ATAC-Seq reads onto the chicken galGal4 reference genome using adjusted parameters below.

Presets very-sensitive-local (Default: false) Alignment score cutoff G,20.0,8.0 (Default: L,-0.6,-0.6) Disable mixed alignments true (Default: false) Disable discordant alignments true (Default: false)

II. I use samtools to remove reads on mitochondrial chromosome (the command line was copied from on-line, because I am not familiar with samtools)

samtools idxstats Sox2_G7_6.bam | cut -f 1 | grep -v chrM | xargs samtools view -b Sox2_G7_6.bam > Sox2_G7_6_woChrM.bam

III. I use macs2 for peak calling, and get an error message below. Totally, I have 18 samples, 15 of them are fine, but 3 of them cannot be finished. The macs2 produced the same error message for the 3 problematic samples (i.e. IndexError: list index out of range). I show one example that sample cannot be done, and one example that sample can be done.

Example: sample cannot be done

gary > macs2 callpeak -t Sox2_G7_6_woChrM.bam --format BAMPE --gsize 930000000 --keep-dup 1 --outdir Sox2_G7_6_woChrM --name Sox2_G7_6_woChrM --verbose 3 --qvalue 0.00001 --cutoff-analysis --call-summits

INFO @ Sun, 11 Dec 2016 07:52:36:

Command line: callpeak -t Sox2_G7_6_woChrM.bam --format BAMPE --gsize 930000000 --keep-dup 1 --outdir Sox2_G7_6_woChrM --name Sox2_G7_6_woChrM --verbose 3 --qvalue 0.00001 --cutoff-analysis --call-summits

ARGUMENTS LIST:

name = Sox2_G7_6_woChrM

format = BAMPE

ChIP-seq file = ['Sox2_G7_6_woChrM.bam']

control file = None

effective genome size = 9.30e+08

band width = 300

model fold = [5, 50]

qvalue cutoff = 1.00e-05

Larger dataset will be scaled towards smaller dataset.

Range for calculating regional lambda is: 10000 bps

Broad region calling is off

Searching for subpeak summits is on

INFO @ Sun, 11 Dec 2016 07:52:36: #1 read fragment files... INFO @ Sun, 11 Dec 2016 07:52:36: #1 read treatment fragments... INFO @ Sun, 11 Dec 2016 07:52:46: 1000000 INFO @ Sun, 11 Dec 2016 07:52:54: 2000000 INFO @ Sun, 11 Dec 2016 07:53:03: 3000000 INFO @ Sun, 11 Dec 2016 07:53:12: 4000000 INFO @ Sun, 11 Dec 2016 07:53:22: 5000000 INFO @ Sun, 11 Dec 2016 07:53:32: 6000000 INFO @ Sun, 11 Dec 2016 07:53:43: 7000000 INFO @ Sun, 11 Dec 2016 07:53:53: 8000000 INFO @ Sun, 11 Dec 2016 07:54:03: 9000000 INFO @ Sun, 11 Dec 2016 07:54:14: 10000000 INFO @ Sun, 11 Dec 2016 07:54:25: 11000000 INFO @ Sun, 11 Dec 2016 07:54:36: 12000000 INFO @ Sun, 11 Dec 2016 07:54:48: 13000000 INFO @ Sun, 11 Dec 2016 07:54:58: 14000000 INFO @ Sun, 11 Dec 2016 07:55:09: 15000000 INFO @ Sun, 11 Dec 2016 07:55:20: 16000000 INFO @ Sun, 11 Dec 2016 07:55:30: 17000000 INFO @ Sun, 11 Dec 2016 07:55:41: 18000000 Traceback (most recent call last): File "/Users/gary/anaconda/bin/macs2", line 614, in <module> main() File "/Users/gary/anaconda/bin/macs2", line 56, in main run( args ) File "/Users/gary/anaconda/lib/python2.7/site-packages/MACS2/callpeak_cmd.py", line 69, in run if options.PE_MODE: (treat, control) = load_frag_files_options (options) File "/Users/gary/anaconda/lib/python2.7/site-packages/MACS2/callpeak_cmd.py", line 355, in load_frag_files_options treat = tp.build_petrack() File "MACS2/IO/Parser.pyx", line 1041, in MACS2.IO.Parser.BAMPEParser.build_petrack (MACS2/IO/Parser.c:14907) File "MACS2/IO/Parser.pyx", line 1045, in MACS2.IO.Parser.BAMPEParser.build_petrack (MACS2/IO/Parser.c:14856) File "MACS2/IO/Parser.pyx", line 1165, in MACS2.IO.Parser.BAMPEParser.__build_petrack_wo_pysam (MACS2/IO/Parser.c:15475) IndexError: list index out of range

Example: sample can be done

gary > macs2 callpeak -t Sox2_G7_4_woChrM.bam --format BAMPE --gsize 930000000 --keep-dup 1 --outdir Sox2_G7_4_woChrM --name Sox2_G7_4_woChrM --verbose 3 --qvalue 0.00001 --cutoff-analysis --call-summits

INFO @ Sun, 11 Dec 2016 06:09:36:

Command line: callpeak -t Sox2_G7_4_woChrM.bam --format BAMPE --gsize 930000000 --keep-dup 1 --outdir Sox2_G7_4_woChrM --name Sox2_G7_4_woChrM --verbose 3 --qvalue 0.00001 --cutoff-analysis --call-summits

ARGUMENTS LIST:

name = Sox2_G7_4_woChrM

format = BAMPE

ChIP-seq file = ['Sox2_G7_4_woChrM.bam']

control file = None

effective genome size = 9.30e+08

band width = 300

model fold = [5, 50]

qvalue cutoff = 1.00e-05

Larger dataset will be scaled towards smaller dataset.

Range for calculating regional lambda is: 10000 bps

Broad region calling is off

Searching for subpeak summits is on

INFO @ Sun, 11 Dec 2016 06:09:36: #1 read fragment files... INFO @ Sun, 11 Dec 2016 06:09:36: #1 read treatment fragments... INFO @ Sun, 11 Dec 2016 06:09:52: 1000000 INFO @ Sun, 11 Dec 2016 06:10:10: 2000000 INFO @ Sun, 11 Dec 2016 06:10:28: 3000000 INFO @ Sun, 11 Dec 2016 06:10:45: 4000000 INFO @ Sun, 11 Dec 2016 06:11:02: 5000000 INFO @ Sun, 11 Dec 2016 06:11:19: 6000000 INFO @ Sun, 11 Dec 2016 06:11:36: 7000000 INFO @ Sun, 11 Dec 2016 06:11:53: 8000000 INFO @ Sun, 11 Dec 2016 06:12:10: 9000000 INFO @ Sun, 11 Dec 2016 06:12:27: 10000000 INFO @ Sun, 11 Dec 2016 06:12:44: 11000000 INFO @ Sun, 11 Dec 2016 06:13:01: 12000000 INFO @ Sun, 11 Dec 2016 06:13:18: 13000000 INFO @ Sun, 11 Dec 2016 06:13:35: 14000000 INFO @ Sun, 11 Dec 2016 06:13:52: 15000000 INFO @ Sun, 11 Dec 2016 06:14:09: 16000000 INFO @ Sun, 11 Dec 2016 06:14:26: 17000000 INFO @ Sun, 11 Dec 2016 06:14:43: 18000000 INFO @ Sun, 11 Dec 2016 06:15:26: #1 mean fragment size is determined as 181 bp from treatment INFO @ Sun, 11 Dec 2016 06:15:26: #1 fragment size = 181 INFO @ Sun, 11 Dec 2016 06:15:26: #1 total fragments in treatment: 18703295 INFO @ Sun, 11 Dec 2016 06:15:26: #1 user defined the maximum fragments... INFO @ Sun, 11 Dec 2016 06:15:26: #1 filter out redundant fragments by allowing at most 1 identical fragment(s) INFO @ Sun, 11 Dec 2016 06:16:30: #1 fragments after filtering in treatment: 17084975 INFO @ Sun, 11 Dec 2016 06:16:30: #1 Redundant rate of treatment: 0.09 INFO @ Sun, 11 Dec 2016 06:16:30: #1 finished! INFO @ Sun, 11 Dec 2016 06:16:30: #2 Build Peak Model... INFO @ Sun, 11 Dec 2016 06:16:30: #2 Skipped... INFO @ Sun, 11 Dec 2016 06:16:30: #2 Use 181 as fragment length INFO @ Sun, 11 Dec 2016 06:16:30: #3 Call peaks... INFO @ Sun, 11 Dec 2016 06:16:30: #3 Going to call summits inside each peak ... INFO @ Sun, 11 Dec 2016 06:16:30: #3 Pre-compute pvalue-qvalue table... INFO @ Sun, 11 Dec 2016 06:16:30: #3 Cutoff will be automatically decided! DEBUG @ Sun, 11 Dec 2016 06:16:30: Start to calculate pvalue stat... DEBUG @ Sun, 11 Dec 2016 06:18:32: access pq hash for 83653408 times INFO @ Sun, 11 Dec 2016 06:18:32: #3 Analysis of cutoff vs num of peaks or total length has been saved in Sox2_G7_4_woChrM/Sox2_G7_4_woChrM_cutoff_analysis.txt INFO @ Sun, 11 Dec 2016 06:18:32: #3 Suggest a cutoff... INFO @ Sun, 11 Dec 2016 06:18:32: #3 -10log10pvalue cutoff 1.00 will call approximately 1 bps regions as significant regions INFO @ Sun, 11 Dec 2016 06:18:32: #3 Call peaks for each chromosome... INFO @ Sun, 11 Dec 2016 06:19:26: #4 Write output xls file... Sox2_G7_4_woChrM/Sox2_G7_4_woChrM_peaks.xls INFO @ Sun, 11 Dec 2016 06:19:26: #4 Write peak in narrowPeak format file... Sox2_G7_4_woChrM/Sox2_G7_4_woChrM_peaks.narrowPeak INFO @ Sun, 11 Dec 2016 06:19:26: #4 Write summits bed file... Sox2_G7_4_woChrM/Sox2_G7_4_woChrM_summits.bed INFO @ Sun, 11 Dec 2016 06:19:26: Done!

atac-seq samtools macs2 bowtie2 • 1.3k views
ADD COMMENTlink modified 22 months ago • written 22 months ago by Gary450

Hello--I want this to be more informative than critical. I find your post fairly hard to read. Biostars now uses markdown in question and comment fields (similar to stackexchange and github). Meta: Markdown in Biostars (now enabled)

ADD REPLYlink written 22 months ago by ariel.balter130
0
gravatar for Gary
22 months ago by
Gary450
Taiwan/Taichung/China Medical University Hospital
Gary450 wrote:

I found that if I sort bam files first (below command), then macs2 can run well. I don’t know why, but it could be related to what Goutham Atla said (https://www.biostars.org/p/127465/). Maybe no reads on chrMt is an issue for macs2.

samtools sort Sox2_G7_6_woChrM.bam Sox2_G7_6_woChrM.sorted

ADD COMMENTlink modified 22 months ago • written 22 months ago by Gary450
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: 2155 users visited in the last hour