BAMPE in MACS2

[Jesper Grud]

Hi,

I am having a problem with -f (--format) BAMPE in MACS2, which produces an error. Running the same file using -f BAM produces no errors. The bam file is query name sorted, so paired end reads follow each other in the file. The output from BAMPE is as following:

INFO  @ Wed, 02 Oct 2013 11:44:29: 

# Command line: callpeak -t /media/Data/Storage/FILE.bam --format BAMPE

# ARGUMENTS LIST:

# name = NA

# format = BAMPE

# ChIP-seq file = ['/media/Data/Storage/FILE.bam']

# control file = None

# effective genome size = 2.70e+09

# band width = 300

# model fold = [5, 50]

# qvalue cutoff = 5.00e-02

# Larger dataset will be scaled towards smaller dataset.

# Range for calculating regional lambda is: 10000 bps

# Broad region calling is off

 

INFO  @ Wed, 02 Oct 2013 11:44:29: #1 read fragment files... 

INFO  @ Wed, 02 Oct 2013 11:44:29: #1 read treatment fragments... 

Traceback (most recent call last):

  File "/usr/local/bin/macs2", line 514, in <module>

    main()

  File "/usr/local/bin/macs2", line 45, in main

    run( args )

  File "/usr/local/lib/python2.7/dist-packages/MACS2/callpeak.py", line 68, in run

    if options.PE_MODE: (treat, control) = load_frag_files_options (options)

  File "/usr/local/lib/python2.7/dist-packages/MACS2/callpeak.py", line 340, in load_frag_files_options

    treat = tp.build_petrack()

  File "cParser.pyx", line 1048, in MACS2.IO.cParser.BAMPEParser.build_petrack (MACS2/IO/cParser.c:13437)

  File "cParser.pyx", line 1052, in MACS2.IO.cParser.BAMPEParser.build_petrack (MACS2/IO/cParser.c:13388)

  File "cParser.pyx", line 1173, in MACS2.IO.cParser.BAMPEParser.__build_petrack_wo_pysam (MACS2/IO/cParser.c:14006)

  File "cPairedEndTrack.pyx", line 70, in MACS2.IO.cPairedEndTrack.PETrackI.add_loc (MACS2/IO/cPairedEndTrack.c:2480)

  File "cPairedEndTrack.pyx", line 85, in MACS2.IO.cPairedEndTrack.PETrackI.add_loc (MACS2/IO/cPairedEndTrack.c:2279)

IndexError: index 0 is out of bounds for axis 0 with size 0

My MACS2 is update. The version info is as following:

jespergrudskatmadsen@ADM103933:~$ macs2 --version

macs2 2.0.10.20130915 (tag:beta)

I am clueless as to why this error occurs. Hope you can help!

[AS]

Hi,

I am getting this exact same error with BAMPE using macs2 2.0.10.20130915. I did not get this error with the previous version I had been using, 2.0.10.20120913. I am running the exact same command I had used for the previous version: macs2 callpeak -t $SAMPLEFILE -c $CONTROLFILE -f BAMPE -n $NAME -g hs -B --trackline

INFO  @ Mon, 14 Oct 2013 12:45:05: #1 read fragment files... 

INFO  @ Mon, 14 Oct 2013 12:45:05: #1 read treatment fragments... 

Traceback (most recent call last):

  File "/apps/source/python-lib/python2.7/macs2/2.0.10.20130915//bin/macs2", line 514, in <module>

    main()

  File "/apps/source/python-lib/python2.7/macs2/2.0.10.20130915//bin/macs2", line 45, in main

    run( args )

  File "/apps/source/python-lib/python2.7/macs2/2.0.10.20130915/lib/python2.7/site-packages/MACS2/callpeak.py", line 68, in run

    if options.PE_MODE: (treat, control) = load_frag_files_options (options)

  File "/apps/source/python-lib/python2.7/macs2/2.0.10.20130915/lib/python2.7/site-packages/MACS2/callpeak.py", line 340, in load_frag_files_options

    treat = tp.build_petrack()

  File "cParser.pyx", line 1048, in MACS2.IO.cParser.BAMPEParser.build_petrack (MACS2/IO/cParser.c:13437)

  File "cParser.pyx", line 1052, in MACS2.IO.cParser.BAMPEParser.build_petrack (MACS2/IO/cParser.c:13388)

  File "cParser.pyx", line 1173, in MACS2.IO.cParser.BAMPEParser.__build_petrack_wo_pysam (MACS2/IO/cParser.c:14006)

  File "cPairedEndTrack.pyx", line 70, in MACS2.IO.cPairedEndTrack.PETrackI.add_loc (MACS2/IO/cPairedEndTrack.c:2480)

  File "cPairedEndTrack.pyx", line 85, in MACS2.IO.cPairedEndTrack.PETrackI.add_loc (MACS2/IO/cPairedEndTrack.c:2279)

IndexError: index 0 is out of bounds for axis 0 with size 0

Does anyone have any suggestions?

[Gabe Zentner]

This isn't a helpful reply (sorry!) but I just got the exact same error when attempting to run MACS2 using BAMPE.  I would appreciate a solution as well.  Thanks!

[Tao Liu]

Hi Gabe, AS and Jesper,

Could you check your PE bam validity by:

$samtools flagstat file.bam

, and post me the results? I tested some PE data locally generated by BWA and samtools, and haven't seen any error.

Could you also let me know how you get the PE bam files?

Best,

Tao Liu

--
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-announcem...@googlegroups.com.
To post to this group, send email to macs-ann...@googlegroups.com.
Visit this group at http://groups.google.com/group/macs-announcement.
For more options, visit https://groups.google.com/groups/opt_out.

[AS]

Hi Tao,

I generate my bam files using bwa. Should we be filtering for proper pairs etc before running MACS? 

Below is the samtools flagstat output for one of my files:

$ samtools flagstat 8_130821_AC26BDACXX_Illumina1.bam

64381396 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

62727462 + 0 mapped (97.43%:-nan%)

64381396 + 0 paired in sequencing

32190698 + 0 read1

32190698 + 0 read2

61364244 + 0 properly paired (95.31%:-nan%)

62249984 + 0 with itself and mate mapped

477478 + 0 singletons (0.74%:-nan%)

502172 + 0 with mate mapped to a different chr

278226 + 0 with mate mapped to a different chr (mapQ>=5)

[Gabe Zentner]

I generated SAM files with Bowtie, then converted them to BAM with Samtools.  Here's a flagstat result for one of my BAM files:

26072890 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

22565148 + 0 mapped (86.55%:-nan%)

26072890 + 0 paired in sequencing

13036445 + 0 read1

13036445 + 0 read2

22565148 + 0 properly paired (86.55%:-nan%)

22565148 + 0 with itself and mate mapped

0 + 0 singletons (0.00%:-nan%)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

[Tao Liu]

Hi AS and Gabe,

Thanks for your information! Now I've found the bug -- my fault. The file MACS2/IO/cPairedEndTrack.pyx missed one line at line #69:

68        self.rlengths = {}

69        self.buffer_size = buffer_size

If you downloaded source code from PyPI before and unpacked then installed by yourself, please replace two files MACS2/IO/cPairedEndTrack.pyx and MACS2/IO/cPairedEndTrack.c with the files I attached in this email.  I've updated PyPI release minutes ago (2.0.10.20130915.1). So you can also do a 'pip install -U macs2' to have your local version updated.

Sorry for the trouble!

Best,

Tao Liu

[AS]

Hi Tao,

BAMPE works now with the updated files, thanks!

Another issue though -- with my dataset (a transcriptional regulator in human NPCs compared to Input control), using version 2.0.10.20120913 I got ~9000 peaks; using version 2.0.10.20120915.1 with the same data, I got only ~3000 peaks. The ~3000 peaks mostly overlap the ~9000 peaks so I would assume that the smaller number are stronger peaks, but have you made changes to the implementation since 2.0.10.20120913 that could produce such a drastic difference in the number of peaks called?

I ran the same command both times: macs2 callpeak -t $SAMPLEFILE -c $CONTROLFILE -f BAMPE -n $NAME -g hs -B --trackline

Thanks

[Tao Liu]

I am sure there is no 2.0.10.20120913. Do you mean 0713 vs 0915.1?  There is no dramatic changes in between these two versions. Perhaps you use -f BAM instead of -f BAMPE before?  I just did some tests and saw no difference.

-Tao

# remove current version then install 2.0.10.20130713...

$macs2 callpeak -t treat.pe.0.005.bam -c input.pe.0.005.bam -f BAMPE -n v713

# remove 2.0.10.20130713 then install 2.0.10.20130915.1...

$macs2 callpeak -t treat.pe.0.005.bam -c input.pe.0.005.bam -f BAMPE -n v915.1

$grep version v713_peaks.xls

# This file is generated by MACS version 2.0.10.20130712 (tag:beta)

$grep version v915.1_peaks.xls

# This file is generated by MACS version 2.0.10.20130915.1 (tag:beta)

$wc -l v713_peaks.xls v915.1_peaks.xls

  5144 v713_peaks.xls

  5144 v915.1_peaks.xls

--
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-announcem...@googlegroups.com.
To post to this group, send email to macs-ann...@googlegroups.com.
Visit this group at http://groups.google.com/group/macs-announcement.
For more options, visit https://groups.google.com/groups/opt_out.

Tao Liu

--
Assistant Professor
Department of Biochemistry
University at Buffalo
NY State Center of Excellence in Bioinformatics & Life Sciences

B2-163 COEBLS
(O) 716-829-2749
tl...@buffalo.edu
http://biomisc.org/

Mailing address:
University at Buffalo-COEBLS
701 Ellicott St, B2-163
Buffalo, NY 14203-1221

[Aarathi Sugathan]

That's strange...

Here is the version number for my older version of MACS2:

[az272@eris1n2 ~]$ module load python/2.7.3

[az272@eris1n2 ~]$ module load numpy-1.7

[az272@eris1n2 ~]$ module use /apps/modulefiles/test

[az272@eris1n2 ~]$ module load macs2/2.0.10

[az272@eris1n2 ~]$ macs2 --version

macs2 2.0.10.20120913 (tag:beta)

And the first few lines of the output .xls file:

# This file is generated by MACS version 2.0.10.20120913 (tag:beta)

# ARGUMENTS LIST:

# name = 60418.v.Input.MACS2

# format = BAMPE

# ChIP-seq file = ['/tmp/880543.tmpdir/8_130821_AC26BDACXX_Illumina3.bam']

# control file = ['/tmp/880543.tmpdir/8_130821_AC26BDACXX_Illumina1.bam']

# effective genome size = 2.70e+09

# band width = 300

# model fold = [5, 50]

# qvalue cutoff = 5.00e-02

[Gabe Zentner]

I'm still getting the same error with the updated files:

ES_FLAG

# format = BAMPE

# ChIP-seq file = ['/home/gzentner/../../cygdrive/z/Jason/Ccnd1/ES_FLAG_2.bam', '/home/gzentner/../../cygdrive/z/Jason/Ccnd1/ES_input.bam']

# control file = None

# effective genome size = 1.87e+09

# band width = 300

# model fold = [5, 50]

# qvalue cutoff = 5.00e-02

# Larger dataset will be scaled towards smaller dataset.

# Range for calculating regional lambda is: 10000 bps

# Broad region calling is off

INFO  @ Tue, 22 Oct 2013 08:54:41: #1 read fragment files...

INFO  @ Tue, 22 Oct 2013 08:54:41: #1 read treatment fragments...

Traceback (most recent call last):

  File "/usr/bin/macs2", line 514, in <module>

    main()

  File "/usr/bin/macs2", line 45, in main

    run( args )

  File "/usr/lib/python2.7/site-packages/MACS2/callpeak.py", line 68, in run

    if options.PE_MODE: (treat, control) = load_frag_files_options (options)

  File "/usr/lib/python2.7/site-packages/MACS2/callpeak.py", line 340, in load_frag_files_options

    treat = tp.build_petrack()

  File "cParser.pyx", line 1048, in MACS2.IO.cParser.BAMPEParser.build_petrack (MACS2/IO/cParser.c:13437)

  File "cParser.pyx", line 1052, in MACS2.IO.cParser.BAMPEParser.build_petrack (MACS2/IO/cParser.c:13388)

  File "cParser.pyx", line 1173, in MACS2.IO.cParser.BAMPEParser.__build_petrack_wo_pysam (MACS2/IO/cParser.c:14006)

  File "cPairedEndTrack.pyx", line 70, in MACS2.IO.cPairedEndTrack.PETrackI.add_loc (MACS2/IO/cPairedEndTrack.c:2480)

  File "cPairedEndTrack.pyx", line 85, in MACS2.IO.cPairedEndTrack.PETrackI.add_loc (MACS2/IO/cPairedEndTrack.c:2279)

IndexError: index (0) out of range (0<=index<0) in dimension 0

On Wednesday, October 2, 2013 2:53:39 AM UTC-7, Jesper Grud wrote:

[Tao Liu]

Hi Gabe,

You may need to remove old installation then reinstall. If you replaced the two files manually, you need to remove the directory named 'build/' under source code directory then install again.

Best,

Tao Liu

--
Assistant Professor
Department of Biochemistry
University at Buffalo
NY State Center of Excellence in Bioinformatics & Life Sciences

B2-163 COEBLS
(O) 716-829-2749
tl...@buffalo.edu

http://biomisc.org/

Mailing address:
University at Buffalo-COEBLS
701 Ellicott St, B2-163
Buffalo, NY 14203-1221

[Tao Liu]

My bad memory! I just checked my records: there is a thing related to BAMPE changed in this year -- the way to calculate local bias. Previously (your version) the entire fragment from input is recognized as a open chromatin bias, but later I modified it to consider ends of input fragment as open chromatin bias. I admit that I haven't done enough test on BAMPE since I don't have enough data to play with. I may over-estimate input bias by 2-folds in this way -- reason you saw fewer peaks now. I will review the codes again. 

By the way, can anyone recommend some published paired-end ChIP-seq data for evaluation? Thanks!

Best,

Tao Liu

--
Assistant Professor
Department of Biochemistry
University at Buffalo
NY State Center of Excellence in Bioinformatics & Life Sciences

B2-163 COEBLS
(O) 716-829-2749
tl...@buffalo.edu

http://biomisc.org/

Mailing address:
University at Buffalo-COEBLS
701 Ellicott St, B2-163
Buffalo, NY 14203-1221

[abk]

Tao,

You can use very deep 101 bp PE data for CTCF, SA1, POL2, PU1 and several histone marks from human LCLs from our recent paper http://www.sciencemag.org/content/early/2013/10/16/science.1242510 . The FASTQs are in GEO GSE50893.

If you want pre-mapped PE BAM files let me know. I can send you links to them.

Thanks,

-A

[Gabe Zentner]

That worked, but I still get this error later on:

/usr/lib/python2.7/site-packages/MACS2/callpeak.py:249: RuntimeWarning: overflow encountered in int_scalars

  peakdetect.call_peaks()

INFO  @ Tue, 15 Oct 2013 12:18:41: #3 pileup treatment data

INFO  @ Tue, 15 Oct 2013 12:18:58: #3 calculate local lambda from control data

INFO  @ Tue, 15 Oct 2013 12:22:03: #3 Build score track ...

INFO  @ Tue, 15 Oct 2013 12:22:59: #3 Calculate qvalues ...

Traceback (most recent call last):

  File "/usr/bin/macs2", line 362, in <module>

    main()

  File "/usr/bin/macs2", line 45, in main

    run( args )

  File "/usr/lib/python2.7/site-packages/MACS2/callpeak.py", line 249, in run

    peakdetect.call_peaks()

  File "cPeakDetect.pyx", line 111, in MACS2.cPeakDetect.PeakDetect.call_peaks (MACS2/cPeakDetect.c:1693)

  File "cPeakDetect.pyx", line 317, in MACS2.cPeakDetect.PeakDetect.__call_peaks_w_control (MACS2/cPeakDetect.c:4369)

  File "cScoreTrack.pyx", line 624, in MACS2.IO.cScoreTrack.scoreTrackII.change_score_method (MACS2/IO/cScoreTrack.c:10282)

  File "cScoreTrack.pyx", line 640, in MACS2.IO.cScoreTrack.scoreTrackII.change_score_method (MACS2/IO/cScoreTrack.c:10092)

  File "cScoreTrack.pyx", line 694, in MACS2.IO.cScoreTrack.scoreTrackII.compute_qvalue (MACS2/IO/cScoreTrack.c:10949)

  File "cScoreTrack.pyx", line 764, in MACS2.IO.cScoreTrack.scoreTrackII.make_pq_table (MACS2/IO/cScoreTrack.c:11321)

OverflowError: Python int too large to convert to C long

[Tao Liu]

Looks like running message from old version. Here is the correct link:

https://pypi.python.org/packages/source/M/MACS2/MACS2-2.0.10.20130915.1.tar.gz

Remember to remove "/usr/lib/python2.7/site-packages/MACS2/“.

Tao

[Gabe Zentner]

I am still getting that error.  I deleted the specified folder and reinstalled from the file you listed; is there something else I need to do?

[Tao Liu]

What error did you get? In current MACS2 codes, cScoreTrack.pyx is not used so I still suspect you haven’t re-install successfully.

Please share us the information on your MACS2 version:

# in your case

$ ls /usr/lib/python2.7/site-packages/MACS2*
$ macs2 —version

Also, please let us know the steps you took to 1. remove the folder; 2. install MACS2. Did you find any error messages while re-installing MACS2?

Tao

[Gabe Zentner]

Result of ls:

/home/gzentner/../../cygdrive/c/cygwin/lib/python2.7/site-packages/MACS2:

__init__.py       bdgdiff.pyc      Constants.pyc    cStat.pyc      OptValidator.py   predictd.pyc

__init__.pyc      bdgpeakcall.py   cPeakDetect.dll  diffpeak.py    OptValidator.pyc  randsample.py

bdgbroadcall.py   bdgpeakcall.pyc  cPeakModel.dll   diffpeak.pyc   OutputWriter.py   randsample.pyc

bdgbroadcall.pyc  callpeak.py      cPileup.dll      filterdup.py   OutputWriter.pyc  refinepeak.py

bdgcmp.py         callpeak.pyc     cProb.dll        filterdup.pyc  pileup.py         refinepeak.pyc

bdgcmp.pyc        cArray.dll       cSignal.dll      hashtable.dll  pileup.pyc

bdgdiff.py        Constants.py     cStat.py         IO             predictd.py

Here is the version I got:

macs2 2.0.10.20130915.1 (tag:beta)

To uninstall, I removed the MACS2 folder from my site-packages directory.  I realize now that I did not remove the EGG-INFO file from the previous installations.  I tried reinstalling after removing both the directory and EGG-INFO file, and it seems to be working fine now (there is no cScoreTrack.pyx in the directory).  However, I just got the following error:

/usr/lib/python2.7/site-packages/MACS2/callpeak.py:254: RuntimeWarning: overflow encountered in int_scalars

  peakdetect.call_peaks()

INFO  @ Thu, 24 Oct 2013 11:11:49: #3 Pre-compute pvalue-qvalue table...

Traceback (most recent call last):

  File "/usr/bin/macs2", line 514, in <module>

    main()

  File "/usr/bin/macs2", line 45, in main

    run( args )

  File "/usr/lib/python2.7/site-packages/MACS2/callpeak.py", line 254, in run

    peakdetect.call_peaks()

  File "cPeakDetect.pyx", line 111, in MACS2.cPeakDetect.PeakDetect.call_peaks (MACS2/cPeakDetect.c:1706)

  File "cPeakDetect.pyx", line 254, in MACS2.cPeakDetect.PeakDetect.__call_peaks_w_control (MACS2/cPeakDetect.c:3301)

  File "cCallPeakUnit.pyx", line 648, in MACS2.IO.cCallPeakUnit.CallerFromAlignments.call_peaks (MACS2/IO/cCallPeakUnit.c:8819)

  File "cCallPeakUnit.pyx", line 668, in MACS2.IO.cCallPeakUnit.CallerFromAlignments.call_peaks (MACS2/IO/cCallPeakUnit.c:8272)

  File "cCallPeakUnit.pyx", line 621, in MACS2.IO.cCallPeakUnit.CallerFromAlignments.__cal_pvalue_qvalue_table (MACS2/IO/cCallPeakUnit.c:7850)

OverflowError: Python int too large to convert to C long

[Tao Liu]

Let me guess: it seems that you are using Windows. Perhaps your cygwin is just a 32-bit version so the python in this cygwin doesn’t support long integer (of your 64bit windows OS)? I am not familiar with windows though.  My suggestion for windows users is, either to have a linux OS (e.g., ubuntu) running within VirtualBox (https://www.virtualbox.org/wiki/Downloads)or to find a linux server then PuTTY to it. Cygwin of MinGW is not a perfect solution.

Best,

Tao

[Tao Liu]

Hi Gabe,

Another suggestion, although I am not sure if it would work, is to 1) install cython in your cygwin environment, 2) remove all *.c files in MACS2 source code folder;3) then use ’setup_w_cython.py’ instead of ’setup.py’ to let cython re-generate .c files from .pyx again.  

Currently, if one runs ’setup.py’, the installer will just compile all the existing .c files within the released package that were created by cython package in my Mac OSX machine. They worked fine in Linux and Mac but I don’t have windows machine so they may have trouble there.

Tao

[Gabe Zentner]

I got access to MACS2 on a local Linux server and everything is working fine.  Thanks for all your help!

[Jesper Grud]

Hi Tao

Have you had time to check this local bias yet? And if not would you recommend newest or older commit or MACS2?

Also there is a small regression, and at least some of it is in the this commit: 23b9cf47681af6539dc116d6d8f17771467ad65f

Using BAMPE without control data automatically specifies scaling of data towards the smallest dataset and callpeak encounters an error at lines 209/210 with c1 using the c1 variable which is not defined, as no control data is used. The error is patched by simply removing lines 209/210, but the real problem is probably the automatic and forced scaling of data in the absense of control data and I cannot figure out where this happens. 

[Tao Liu]

Hi Jesper,

Sorry for the late reply! I just fixed the bug with using BAMPE in the absence of control. The newest release can be found at pypi:

https://pypi.python.org/pypi/MACS2/2.0.10.20131216

It’s just a stupid mistake from me… :  https://github.com/taoliu/MACS/commit/46cb25d1e6dc3b202a12818928d02f0f042fc04a

Best,

Tao Liu

[wendy]

But, when I run the updated macs2 (20131216), "Command line: callpeak -t test.bam -f BAMPE -n test -g mm -B", I still got errors:

INFO  @ Wed, 18 Dec 2013 01:22:58: #3 Call peaks...

Traceback (most recent call last):

  File "/home/liu/macs21216/bin/macs2", line 557, in <module>
    main()
  File "/home/liu/macs21216/bin/macs2", line 56, in main
    run( args )
  File "/home/liu/local/lib/python2.7/site-packages/MACS2/callpeak.py", line 210, in run
    c1 = c1 * 2
UnboundLocalError: local variable 'c1' referenced before assignment

I just don't know where the problem is? Can anyone do me a favor?
Thank you very much.
Jinhua

[Tao Liu]

Hi Wendy,

Apparently, your executable macs2 script is in "/home/liu/macs21216/bin/“, and libraries in "/home/liu/local/lib/python2.7/site-packages/MACS2/“ that indicates you messed up the installation process. Please do not install MACS2 to its source code folder. Remove everything in /home/liu/local/lib/python2.7/site-packages/MACS2* and all your macs2 executable scripts. Then follow this instruction: https://github.com/taoliu/MACS/wiki/Install-macs2

I highly recommend using PyPI to install python libraries. 

Best,

Tao Liu

--
Assistant Professor
Department of Biochemistry

Department of Biostatistics

University at Buffalo
NY State Center of Excellence in Bioinformatics & Life Sciences

B2-163 COEBLS
(O) 716-829-2749
tl...@buffalo.edu

http://biomisc.org/

Mailing address:
University at Buffalo-COEBLS
701 Ellicott St, B2-163
Buffalo, NY 14203-1221

[AS]

Hi Tao,

Thank you for continuing to update this package.

I wanted to clarify the change you made for BAMPE. In your message earlier in this thread in response my question about getting many fewer peaks with version 2.0.10.20120915.1 compared to version 2.0.10.20120913, you had said, "there is a thing related to BAMPE changed in this year -- the way to calculate local bias. Previously (your version) the entire fragment from input is recognized as a open chromatin bias, but later I modified it to consider ends of input fragment as open chromatin bias. I admit that I haven't done enough test on BAMPE since I don't have enough data to play with. I may over-estimate input bias by 2-folds in this way" --- This sounds like the difference is in how you treat the input control samples.

But in your latest message about version 2.0.10.20131216, you said the bug was in "using BAMPE in the absence of control."

Could you confirm whether 20131216 makes any changes in how BAMPE reads are treated when there is an input control?

Another question, how does MACS2 deal with singletons, pairs mapped to different chromosomes, etc, in BAMPE?

Thanks.

[Tao Liu]

Hi Aarathi,

On Jan 7, 2014, at 4:02 PM, AS <aarathi....@gmail.com> wrote:

I wanted to clarify the change you made for BAMPE. In your message earlier in this thread in response my question about getting many fewer peaks with version 2.0.10.20120915.1 compared to version 2.0.10.20120913, you had said, "there is a thing related to BAMPE changed in this year -- the way to calculate local bias. Previously (your version) the entire fragment from input is recognized as a open chromatin bias, but later I modified it to consider ends of input fragment as open chromatin bias. I admit that I haven't done enough test on BAMPE since I don't have enough data to play with. I may over-estimate input bias by 2-folds in this way" --- This sounds like the difference is in how you treat the input control samples.

About this, the difference is the normalization factor for input control. 

But in your latest message about version 2.0.10.20131216, you said the bug was in "using BAMPE in the absence of control."

It's another bug, sorry :) When control is absent, some versions of MACS2 would throw exception and halt.

Could you confirm whether 20131216 makes any changes in how BAMPE reads are treated when there is an input control?

There is no significant thing fixed in Dec 2013 version comparing to Sept 2013 version regarding to BAMPE, except for the control issue. But perhaps, the fixed smoothing function will change the summit location and scores a bit. 

Another question, how does MACS2 deal with singletons, pairs mapped to different chromosomes, etc, in BAMPE?

Only successfully mapped pairs are used in BAMPE. Singletons will be ignored.

Best,

Tao

[AS]

Thanks for the clarification about absence of control.

So if I understand correctly, in the presence of control, there should still be a discrepancy in the number of peaks called using version 20120913 compared to the later versions, because in the former, "the entire fragment from input is recognized as open chromatin bias", whereas the latter versions "consider ends of input fragment as open chromatin bias", resulting in "over-estimate input bias by 2-folds"?

In other words, in the Dec 2013 version, do you consider open chromatin bias to be the entire input fragment or the ends of input fragment? And what about the ChIP fragments, do you consider their ends or the entire fragment?

And finally, how does MACS2 deal with PE reads where the mates map to different chromosomes?

Thanks!

[梁芳]

Dear Tao,

We have some questions when using your software macs2 2.0.10.20120913 (tag:beta).

For the command "/software/biosoft/software/python2.7.2/bin/macs2 callpeak":

1)  How to set the -s value if the length of the reads are different after we remove the adaptor sequences?

2)  How to set the --bw value if the size of degraded fragments is a range from 100 to 300bp? 

Thanks for your time!

Sincerely,

Fang Liang

-----原始邮件-----
发件人:"Tao Liu" <vladim...@gmail.com>
发送时间:2014-01-08 22:03:28 (星期三)
收件人: "macs-ann...@googlegroups.com announcement" <macs-ann...@googlegroups.com>
抄送:
主题: Re: [macs-announscement] BAMPE in MACS2

--

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+unsub...@googlegroups.com.

[hash]

Hi Tao,

I am seeing the same error with the latest version of macs2. I am trying to call peaks in BAMPE mode with two treatments and two controls:

macs-2.1.0_20150420/bin/macs2 callpeak --verbose=2 -t <TREAT_BAM1> <TREAT_BAM2> --format=BAMPE --gsize=hs --outdir=OUTDIR --name=NAME -c <CONTROL_BAM1> <CONTROL_BAM2> --keep-dup=all --qvalue=0.05 --cutoff-analysis 

It is the latest version macs-2.1.0_20150420.

ERROR from sysout file:

# effective genome size = 2.70e+09

# band width = 300

# model fold = [5, 50]

# qvalue cutoff = 5.00e-02

# Larger dataset will be scaled towards smaller dataset.

# Range for calculating regional lambda is: 1000 bps and 10000 bps

# Broad region calling is off

 

INFO  @ Wed, 20 May 2015 13:46:06: #1 read fragment files... 

INFO  @ Wed, 20 May 2015 13:46:06: #1 read treatment fragments... 

INFO  @ Wed, 20 May 2015 13:46:20:  1000000 

INFO  @ Wed, 20 May 2015 13:46:34:  2000000 

INFO  @ Wed, 20 May 2015 13:46:49:  3000000 

....

INFO  @ Wed, 20 May 2015 14:12:52:  112000000 

INFO  @ Wed, 20 May 2015 14:13:05:  113000000 

INFO  @ Wed, 20 May 2015 14:13:19:  114000000 

Traceback (most recent call last):

  File "/farm/babs/redhat6/software/macs-2.1.0_20150420/bin/macs2", line 4, in <module>

    __import__('pkg_resources').run_script('MACS2==2.1.0.20150420', 'macs2')

  File "/farm/babs/redhat6/software/python/lib/python2.7/site-packages/pkg_resources/__init__.py", line 729, in run_script

    self.require(requires)[0].run_script(script_name, ns)

  File "/farm/babs/redhat6/software/python/lib/python2.7/site-packages/pkg_resources/__init__.py", line 1642, in run_script

    exec(code, namespace, namespace)

  File "/farm/babs/redhat6/software/macs-2.1.0_20150420/lib/python2.7/site-packages/MACS2-2.1.0.20150420-py2.7-linux-x86_64.egg/EGG-INFO/scripts/macs2", line 614, in <module>

    main()

  File "/farm/babs/redhat6/software/macs-2.1.0_20150420/lib/python2.7/site-packages/MACS2-2.1.0.20150420-py2.7-linux-x86_64.egg/EGG-INFO/scripts/macs2", line 56, in main

    run( args )

  File "/farm/babs/redhat6/software/macs-2.1.0_20150420/lib/python2.7/site-packages/MACS2-2.1.0.20150420-py2.7-linux-x86_64.egg/MACS2/callpeak_cmd.py", line 69, in run

    if options.PE_MODE: (treat, control) = load_frag_files_options (options)

  File "/farm/babs/redhat6/software/macs-2.1.0_20150420/lib/python2.7/site-packages/MACS2-2.1.0.20150420-py2.7-linux-x86_64.egg/MACS2/callpeak_cmd.py", line 361, in load_frag_files_options

    treat = tp.append_petrack( treat )

  File "MACS2/IO/Parser.pyx", line 1170, in MACS2.IO.Parser.BAMPEParser.append_petrack (MACS2/IO/Parser.c:16632)

  File "MACS2/IO/Parser.pyx", line 1208, in MACS2.IO.Parser.BAMPEParser.append_petrack (MACS2/IO/Parser.c:16426)

  File "MACS2/IO/PairedEndTrack.pyx", line 76, in MACS2.IO.PairedEndTrack.PETrackI.add_loc (MACS2/IO/PairedEndTrack.c:2293)

  File "MACS2/IO/PairedEndTrack.pyx", line 93, in MACS2.IO.PairedEndTrack.PETrackI.add_loc (MACS2/IO/PairedEndTrack.c:2131)

IndexError: index 9350002 is out of bounds for axis 0 with size 9350002

It seems to fail after reading in the first bam file. If I run macs2 on individual treatment vs control bam files in BAMPE mode it works.

The bam files are sorted by position and I have already removed duplicate reads with Picard MarkDuplicates. I have then filtered for proper pairs and removed any singleton reads. 

Samtools flagstat for one of the bam files:

188374238 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 secondary

0 + 0 supplementary

0 + 0 duplicates

188374238 + 0 mapped (100.00%:-nan%)

188374238 + 0 paired in sequencing

94187119 + 0 read1

94187119 + 0 read2

188289616 + 0 properly paired (99.96%:-nan%)

188374238 + 0 with itself and mate mapped

0 + 0 singletons (0.00%:-nan%)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

Any help would be much appreciated.

Many Thanks,

Harshil

On Wednesday, 2 October 2013 10:53:39 UTC+1, Jesper Grud wrote:

Hi,

I am having a problem with -f (--format) BAMPE in MACS2, which produces an error. Running the same file using -f BAM produces no errors. The bam file is query name sorted, so paired end reads follow each other in the file. The output from BAMPE is as following:

INFO  @ Wed, 02 Oct 2013 11:44:29: 

# Command line: callpeak -t /media/Data/Storage/FILE.bam --format BAMPE

# ARGUMENTS LIST:

# name = NA

# format = BAMPE

# ChIP-seq file = ['/media/Data/Storage/FILE.bam']

# control file = None

# effective genome size = 2.70e+09

# band width = 300

# model fold = [5, 50]

# qvalue cutoff = 5.00e-02

# Larger dataset will be scaled towards smaller dataset.

# Range for calculating regional lambda is: 10000 bps

# Broad region calling is off

 

INFO  @ Wed, 02 Oct 2013 11:44:29: #1 read fragment files... 

INFO  @ Wed, 02 Oct 2013 11:44:29: #1 read treatment fragments... 

Traceback (most recent call last):

  File "/usr/local/bin/macs2", line 514, in <module>

    main()

  File "/usr/local/bin/macs2", line 45, in main

    run( args )

  File "/usr/local/lib/python2.7/dist-packages/MACS2/callpeak.py", line 68, in run

    if options.PE_MODE: (treat, control) = load_frag_files_options (options)

  File "/usr/local/lib/python2.7/dist-packages/MACS2/callpeak.py", line 340, in load_frag_files_options

    treat = tp.build_petrack()

  File "cParser.pyx", line 1048, in MACS2.IO.cParser.BAMPEParser.build_petrack (MACS2/IO/cParser.c:13437)

  File "cParser.pyx", line 1052, in MACS2.IO.cParser.BAMPEParser.build_petrack (MACS2/IO/cParser.c:13388)

  File "cParser.pyx", line 1173, in MACS2.IO.cParser.BAMPEParser.__build_petrack_wo_pysam (MACS2/IO/cParser.c:14006)

  File "cPairedEndTrack.pyx", line 70, in MACS2.IO.cPairedEndTrack.PETrackI.add_loc (MACS2/IO/cPairedEndTrack.c:2480)

  File "cPairedEndTrack.pyx", line 85, in MACS2.IO.cPairedEndTrack.PETrackI.add_loc (MACS2/IO/cPairedEndTrack.c:2279)

IndexError: index 0 is out of bounds for axis 0 with size 0