Macs2 callpeak ValueError: max() arg is an empty sequence

[julie dubois]

Hi,

I have an intermitent problem with Macs2 callpeak when I use build model. It seems a problem with building model because I have no problem when I use --nomodel and a fixed extend size in the same sample
So I ran this command with MACS2 2.1.0.20151222 :
callpeak --name MACS2 -t treat.bam --format=BAM --gsize 2451960000 --bw=300 --ratio 1.0 --slocal 1000 --llocal 10000 --keep-dup 2 --bdg --pvalue 0.01
# ARGUMENTS LIST:
# name = MACS2
# format = BAM
# ChIP-seq file = ['treat.bam']
# control file = None
# effective genome size = 2.45e+09
# band width = 300
# model fold = [5, 50]
# pvalue cutoff = 1.00e-02
# qvalue will not be calculated and reported as -1 in the final output.
# Larger dataset will be scaled towards smaller dataset.
# Range for calculating regional lambda is: 10000 bps
# Broad region calling is off
 
INFO  @ Fri, 15 Sep 2017 11:11:27: #1 read tag files...
INFO  @ Fri, 15 Sep 2017 11:11:27: #1 read treatment tags...
INFO  @ Fri, 15 Sep 2017 11:11:29:  1000000
INFO  @ Fri, 15 Sep 2017 11:11:31:  2000000
INFO  @ Fri, 15 Sep 2017 11:11:33:  3000000
...
...
...
INFO  @ Fri, 15 Sep 2017 11:13:00:  48000000
INFO  @ Fri, 15 Sep 2017 11:13:02:  49000000
INFO  @ Fri, 15 Sep 2017 11:13:03:  50000000
INFO  @ Fri, 15 Sep 2017 11:13:04: #1 tag size is determined as 16 bps
INFO  @ Fri, 15 Sep 2017 11:13:04: #1 tag size = 16
INFO  @ Fri, 15 Sep 2017 11:13:04: #1  total tags in treatment: 50000000
INFO  @ Fri, 15 Sep 2017 11:13:04: #1 user defined the maximum tags...
INFO  @ Fri, 15 Sep 2017 11:13:04: #1 filter out redundant tags at the same location and the same strand by allowing at most 2 tag(s)
INFO  @ Fri, 15 Sep 2017 11:13:10: #1  tags after filtering in treatment: 50000000
INFO  @ Fri, 15 Sep 2017 11:13:10: #1  Redundant rate of treatment: 0.00
INFO  @ Fri, 15 Sep 2017 11:13:10: #1 finished!
INFO  @ Fri, 15 Sep 2017 11:13:10: #2 Build Peak Model...
INFO  @ Fri, 15 Sep 2017 11:13:18: #2 number of paired peaks: 9006
INFO  @ Fri, 15 Sep 2017 11:13:18: start model_add_line...
INFO  @ Fri, 15 Sep 2017 11:13:31: start X-correlation...

Traceback (most recent call last):
  File "/data/galaxy-dist/tool_dependency/macs2/2.1.0.20151222/iuc/package_macs2_2_1_0_20151222/e1370f7d5e2f/bin/macs2", line 4, in <module>
    __import__('pkg_resources').run_script('MACS2==2.1.0.20151222', 'macs2')
  File "/data/galaxy-dist/.venv/local/lib/python2.7/site-packages/pkg_resources/__init__.py", line 748, in run_script
    self.require(requires)[0].run_script(script_name, ns)
  File "/data/galaxy-dist/.venv/local/lib/python2.7/site-packages/pkg_resources/__init__.py", line 1524, in run_script
    exec(script_code, namespace, namespace)
  File "/data/galaxy-dist/tool_dependency/macs2/2.1.0.20151222/iuc/package_macs2_2_1_0_20151222/e1370f7d5e2f/lib/python/MACS2-2.1.0.20151222-py2.7-linux-x86_64.egg/EGG-INFO/scripts/macs2", line 614, in <module>
   
  File "/data/galaxy-dist/tool_dependency/macs2/2.1.0.20151222/iuc/package_macs2_2_1_0_20151222/e1370f7d5e2f/lib/python/MACS2-2.1.0.20151222-py2.7-linux-x86_64.egg/EGG-INFO/scripts/macs2", line 56, in main
   
  File "build/bdist.linux-x86_64/egg/MACS2/callpeak_cmd.py", line 175, in run
  File "MACS2/PeakModel.pyx", line 108, in MACS2.PeakModel.PeakModel.__init__ (MACS2/PeakModel.c:2385)
  File "MACS2/PeakModel.pyx", line 150, in MACS2.PeakModel.PeakModel.build (MACS2/PeakModel.c:2953)
  File "MACS2/PeakModel.pyx", line 230, in MACS2.PeakModel.PeakModel.__paired_peak_model (MACS2/PeakModel.c:4209)
ValueError: max() arg is an empty sequence


Any idea ?

Thanks
Julie

[Moshe Olshansky]

I would try the latest version of MACS2 (can be downloaded from Github).

--
You received this message because you are subscribed to the Google Groups "MACS announcement" group.
To unsubscribe from this group and stop receiving emails from it, send an email to macs-announcement+unsubscribe@googlegroups.com.
To post to this group, send email to macs-announcement@googlegroups.com.
Visit this group at https://groups.google.com/group/macs-announcement.
For more options, visit https://groups.google.com/d/optout.

[julie dubois]

I've been tried the latest version of MACS2 : macs2 2.1.1.20160309, with same parameters, and I have exactly the same issue ! with little changes in number of lines.
I've tried too with other mfold values, other bandwith but I've obtained the same behaviour.

...
INFO  @ Mon, 18 Sep 2017 11:19:57: start model_add_line...
INFO  @ Mon, 18 Sep 2017 11:19:58: start X-correlation...

Traceback (most recent call last):

  File "/usr/local/bin/macs2", line 617, in <module>
    main()
  File "/usr/local/bin/macs2", line 101, in main
    run( args )
  File "/usr/local/lib/python2.7/dist-packages/MACS2/predictd_cmd.py", line 65, in run
    opt = options
  File "MACS2/PeakModel.pyx", line 108, in MACS2.PeakModel.PeakModel.__init__ (MACS2/PeakModel.c:2438)
  File "MACS2/PeakModel.pyx", line 150, in MACS2.PeakModel.PeakModel.build (MACS2/PeakModel.c:3010)
  File "MACS2/PeakModel.pyx", line 230, in MACS2.PeakModel.PeakModel.__paired_peak_model (MACS2/PeakModel.c:4266)

ValueError: max() arg is an empty sequence

Thanks in advance
Julie

[Moshe Olshansky]

Hi Julie,

I have never seen such an error.

I think that 16 bp long reads are unusual but I do not think that this caused the problem. I hope that Tao will say something.

One possibility is to estimate mean fragment length using some other tool and then run macs2 with nomodel using that fragment length.

What operating system are you using? If you wish you can send me a link so that I can download your file and see whether I also get this error.

Best regards,

Moshe.

On Fri, Sep 15, 2017 at 7:42 PM, julie dubois <dubj...@gmail.com> wrote:

--

[julie dubois]

Hi Moshe,
I don't think that the length of reads is the problem too. I have the same error with 20bp reads from ENCODE and certain datasets passed the building of model and others didn't pass it !

I use Ubuntu 14.04 and 16.04 in two different computers (respectively 32 and 64 G of RAM)
Here it is the link to download my dataset with reads of 16bp, if you want to test.
https://drive.google.com/file/d/0Bzb4P7gCnNAlR2tMWmJ1ZnF0STA/view?usp=sharing

Thanks
Julie

[Moshe Olshansky]

Hi Julie,

I downloaded your file and get same error both on my CentOS 7 laptop and our Ubuntu 16.04 server. So it really looks like MACS2 bug and I hope that Tao comments on this.

I converted your bam file to bed file (bedtools bamtobed) and got same error. However when I only took reads aligned to chr11 and ran MACS2 with same parameters I got no error and fragment size was 296. Just out of curiosity I then ran MACS2 again but now making genome size 135 millions (roughly the highest position in your chr11). Now the predicted fragment size was 15 bp, which is definitely wrong. Then I ran MACS2 with the complete bed file (50 M reads) but making the library size 2.45e10. I got no error and the predicted fragment size was 294 bp.

So probably you should run MACS@ first with extremely large genome size (and without the -bdg option) to get the fragment size and then use it with nomodel and correct genome size to get the results (of course using very large genome size will make the background very low and VERY many regions will become peaks).

It really looks like MACS2 bug and I hope that it will be corrected.

Best regards,

Moshe.

[julie dubois]

Dear  Moshe,

Thanks a lot for the time you took for this issue.
I hope Tao can explain this bug too.

Thanks again.

Julie