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!
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)