From now on, I won't include any C codes, generated by Cython wrapper, in github. So you can't just pull github snapshot then run setup.py to install MACS2. Instead, you can install Cython then use setup_w_cython.py to install MACS2 or just stick to PyPI release: https://pypi.python.org/pypi/MACS2

from macs.

OK, thanks for the info. By the way, python setup.py install works from an
earlier version,

commit ca80653
Author: liqingtian [email protected]
Date: Mon Jan 27 00:44:29 2014 -0500

Update README.rst

On Thu, May 1, 2014 at 3:12 PM, Tao Liu (τν) [email protected]:

Closed #36 #36.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/36
.

from macs.

Right. I just removed C codes in recent commit.

Tao

On May 1, 2014, at 3:15 PM, Toby Dylan Hocking [email protected] wrote:

OK, thanks for the info. By the way, python setup.py install works from an
earlier version,

commit ca80653
Author: liqingtian [email protected]
Date: Mon Jan 27 00:44:29 2014 -0500

Update README.rst

On Thu, May 1, 2014 at 3:12 PM, Tao Liu (τν) [email protected]:

Closed #36 #36.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/36
.

—
Reply to this email directly or view it on GitHub.

from macs.

I am using MACS2-2.0.10.20131216.
I want to debug code by print statements.

For Instance I put a line:
print "------------------------------------------------------------"
in init() function of cPeakModel.pyx

Then I reistalled macs2 by running
python setup.py install

and I ran again:
macs2 callpeak -t ~/macs\ input/something.sam

but that line did not appear on console.

Can you please help?
Am I missing something.

from macs.

Could you post the entire error message? Also, could you provide the version of macs? Thanks!

On Jun 20, 2014, at 5:25 AM, Kuldeep Yadav [email protected] wrote:

While running following
./macs2 callpeak -t ~/macs\ input/something.sam

it shows error
'Namespace' object has no attribute 'shiftsize'

—
Reply to this email directly or view it on GitHub.

from macs.

I will be able to send that by monday, but I do remember version
It is latest 20140616
On Jun 20, 2014 10:12 PM, "Tao Liu (τν)" [email protected] wrote:

Could you post the entire error message? Also, could you provide the
version of macs? Thanks!

On Jun 20, 2014, at 5:25 AM, Kuldeep Yadav [email protected]
wrote:

While running following
./macs2 callpeak -t ~/macs\ input/something.sam

it shows error
'Namespace' object has no attribute 'shiftsize'

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#36 (comment).

from macs.

If you installed it through pip, then to call the "macs2" in current working folder as "./macs2" should not be used, since pip should have put executable file in system path. You may installed a new MACS2 but invoked an older "macs2" script. Note, shiftsize option has been replaced by extsize optioning the latest version.

Tao

On Jun 20, 2014, at 12:49 PM, Kuldeep Yadav [email protected] wrote:

I will be able to send that by monday, but I do remember version
It is latest 20140616
On Jun 20, 2014 10:12 PM, "Tao Liu (¦Ó¦Í)" [email protected] wrote:

Could you post the entire error message? Also, could you provide the
version of macs? Thanks!

On Jun 20, 2014, at 5:25 AM, Kuldeep Yadav [email protected]
wrote:

While running following
./macs2 callpeak -t ~/macs\ input/something.sam

it shows error
'Namespace' object has no attribute 'shiftsize'

¡ª
Reply to this email directly or view it on GitHub.

¡ª
Reply to this email directly or view it on GitHub
#36 (comment).

¡ª
Reply to this email directly or view it on GitHub.

from macs.

Thank You.

One more problem I am facing. How to debug .pyx code?

I want to print statements to check values of variables in .pyx files. So I
made the changes to files and reinstalled macs with following command:
python setup.py install

but when I run it, I don't see any prints in terminal output. That is, no
changes in console/terminal logs.
Am I missing something in installing steps?

P.S. I made changes in init() of cPeakDetect.pyx file. There is no way
that this function is skipped while running.

On Fri, Jun 20, 2014 at 10:24 PM, Tao Liu (τν) [email protected]
wrote:

If you installed it through pip, then to call the "macs2" in current
working folder as "./macs2" should not be used, since pip should have put
executable file in system path. You may installed a new MACS2 but invoked
an older "macs2" script. Note, shiftsize option has been replaced by
extsize optioning the latest version.

Tao

On Jun 20, 2014, at 12:49 PM, Kuldeep Yadav [email protected]
wrote:

I will be able to send that by monday, but I do remember version
It is latest 20140616
On Jun 20, 2014 10:12 PM, "Tao Liu (¦Ó¦Í)" [email protected]
wrote:

Could you post the entire error message? Also, could you provide the
version of macs? Thanks!

On Jun 20, 2014, at 5:25 AM, Kuldeep Yadav [email protected]

wrote:

While running following
./macs2 callpeak -t ~/macs\ input/something.sam

it shows error
'Namespace' object has no attribute 'shiftsize'

¡ª
Reply to this email directly or view it on GitHub.

¡ª
Reply to this email directly or view it on GitHub
#36 (comment).

¡ª
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

You need to install ‘cython’ in order to ‘translate’ pyx files to ‘c’ codes. Then you need ’setup_w_cython.py’ script to compile and install your new version. You can find ’setup_w_cython.py’ on MACS2 github site.

Tao

On Jun 20, 2014, at 1:13 PM, Kuldeep Yadav [email protected] wrote:

Thank You.

One more problem I am facing. How to debug .pyx code?

I want to print statements to check values of variables in .pyx files. So I
made the changes to files and reinstalled macs with following command:
python setup.py install

but when I run it, I don't see any prints in terminal output. That is, no
changes in console/terminal logs.
Am I missing something in installing steps?

P.S. I made changes in init() of cPeakDetect.pyx file. There is no way
that this function is skipped while running.

On Fri, Jun 20, 2014 at 10:24 PM, Tao Liu (τν) [email protected]
wrote:

If you installed it through pip, then to call the "macs2" in current
working folder as "./macs2" should not be used, since pip should have put
executable file in system path. You may installed a new MACS2 but invoked
an older "macs2" script. Note, shiftsize option has been replaced by
extsize optioning the latest version.

Tao

On Jun 20, 2014, at 12:49 PM, Kuldeep Yadav [email protected]
wrote:

I will be able to send that by monday, but I do remember version
It is latest 20140616
On Jun 20, 2014 10:12 PM, "Tao Liu (¦Ó¦Í)" [email protected]
wrote:

Could you post the entire error message? Also, could you provide the
version of macs? Thanks!

On Jun 20, 2014, at 5:25 AM, Kuldeep Yadav [email protected]

wrote:

While running following
./macs2 callpeak -t ~/macs\ input/something.sam

it shows error
'Namespace' object has no attribute 'shiftsize'

¡ª
Reply to this email directly or view it on GitHub.

¡ª
Reply to this email directly or view it on GitHub
#36 (comment).

¡ª
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA
—
Reply to this email directly or view it on GitHub.

from macs.

I have already installed all prerequisite like cython, gcc etc. I will try setup_w_cython.py.
Thank you.

On 20-Jun-2014, at 10:45 pm, Tao Liu (τν) [email protected] wrote:

You need to install ‘cython’ in order to ‘translate’ pyx files to ‘c’ codes. Then you need ’setup_w_cython.py’ script to compile and install your new version. You can find ’setup_w_cython.py’ on MACS2 github site.

Tao

On Jun 20, 2014, at 1:13 PM, Kuldeep Yadav [email protected] wrote:

Thank You.

One more problem I am facing. How to debug .pyx code?

I want to print statements to check values of variables in .pyx files. So I
made the changes to files and reinstalled macs with following command:
python setup.py install

but when I run it, I don't see any prints in terminal output. That is, no
changes in console/terminal logs.
Am I missing something in installing steps?

P.S. I made changes in init() of cPeakDetect.pyx file. There is no way
that this function is skipped while running.

On Fri, Jun 20, 2014 at 10:24 PM, Tao Liu (τν) [email protected]
wrote:

If you installed it through pip, then to call the "macs2" in current
working folder as "./macs2" should not be used, since pip should have put
executable file in system path. You may installed a new MACS2 but invoked
an older "macs2" script. Note, shiftsize option has been replaced by
extsize optioning the latest version.

Tao

On Jun 20, 2014, at 12:49 PM, Kuldeep Yadav [email protected]
wrote:

I will be able to send that by monday, but I do remember version
It is latest 20140616
On Jun 20, 2014 10:12 PM, "Tao Liu (¦Ó¦Í)" [email protected]
wrote:

Could you post the entire error message? Also, could you provide the
version of macs? Thanks!

On Jun 20, 2014, at 5:25 AM, Kuldeep Yadav [email protected]

wrote:

While running following
./macs2 callpeak -t ~/macs\ input/something.sam

it shows error
'Namespace' object has no attribute 'shiftsize'

¡ª
Reply to this email directly or view it on GitHub.

¡ª
Reply to this email directly or view it on GitHub
#36 (comment).

¡ª
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA
—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub.

from macs.

Hello Tao, I want to translate MACS python code to java. So I am reading the python code and trying to understand algorithm. And I have questions on some points.

1. In functions __call_peak_wo_control() in file cPeakDetect.pyx, why did you took ‘d’ to be zero for paired end mode.

   if self.PE_MODE: d = 0

P.S. please send me reading material if you have anything related to the algorithm beside the MACS research paper. But that paper doesn’t explains the implementation subtilities and also t does not have pseudo code.

Thank you.

On 20-Jun-2014, at 10:54 pm, Kuldeep Yadav [email protected] wrote:

I have already installed all prerequisite like cython, gcc etc. I will try setup_w_cython.py.
Thank you.

On 20-Jun-2014, at 10:45 pm, Tao Liu (τν) [email protected] wrote:

You need to install ‘cython’ in order to ‘translate’ pyx files to ‘c’ codes. Then you need ’setup_w_cython.py’ script to compile and install your new version. You can find ’setup_w_cython.py’ on MACS2 github site.

Tao

On Jun 20, 2014, at 1:13 PM, Kuldeep Yadav [email protected] wrote:

Thank You.

One more problem I am facing. How to debug .pyx code?

I want to print statements to check values of variables in .pyx files. So I
made the changes to files and reinstalled macs with following command:
python setup.py install

but when I run it, I don't see any prints in terminal output. That is, no
changes in console/terminal logs.
Am I missing something in installing steps?

P.S. I made changes in init() of cPeakDetect.pyx file. There is no way
that this function is skipped while running.

On Fri, Jun 20, 2014 at 10:24 PM, Tao Liu (τν) [email protected]
wrote:

If you installed it through pip, then to call the "macs2" in current
working folder as "./macs2" should not be used, since pip should have put
executable file in system path. You may installed a new MACS2 but invoked
an older "macs2" script. Note, shiftsize option has been replaced by
extsize optioning the latest version.

Tao

On Jun 20, 2014, at 12:49 PM, Kuldeep Yadav [email protected]
wrote:

I will be able to send that by monday, but I do remember version
It is latest 20140616
On Jun 20, 2014 10:12 PM, "Tao Liu (¦Ó¦Í)" [email protected]
wrote:

Could you post the entire error message? Also, could you provide the
version of macs? Thanks!

On Jun 20, 2014, at 5:25 AM, Kuldeep Yadav [email protected]

wrote:

While running following
./macs2 callpeak -t ~/macs\ input/something.sam

it shows error
'Namespace' object has no attribute 'shiftsize'

¡ª
Reply to this email directly or view it on GitHub.

¡ª
Reply to this email directly or view it on GitHub
#36 (comment).

¡ª
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA
—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub.

from macs.

Waiting for your reply sir. One more doubt:

file: cFixWidthTrack.pyx

cpdef pileup_a_chromosome ( self, str chrom, list ds, list

scale_factor_s, float baseline_value = 0.0, bint directional = True, bint
halfextension = True ):
"""pileup a certain chromosome, return [p,v](end position and
value) list.

    ds             : tag will be extended to this value to 3' direction,
                     unless directional is False. Can contain multiple

extension
values. Final pileup will the maximum.
scale_factor_s : linearly scale the pileup value applied to each d
in ds. The list should have the same length as ds.
baseline_value : a value to be filled for missing values, and will
be the minimum pileup.
directional : if False, the strand or direction of tag will be
ignored, so that extenstion will be both sides with d/2.
halfextension : only make a fragment of d/2 size centered at
fragment center
"""
cdef:
long d
long five_shift, three_shift # adjustment to 5' end and 3' end
positions to make a fragment
dict chrlengths = self.get_rlengths ()
long rlength = chrlengths[chrom]
object ends
list five_shift_s = []
list three_shift_s = []
list tmp_pileup, prev_pileup

    assert len(ds) == len(scale_factor_s), "ds and scale_factor_s must

have the same length!"

    # adjust extension length according to 'directional' and

'halfextension' setting.
for d in ds:
if directional:
# only extend to 3' side
if halfextension:
five_shift_s.append(d/-4) # five shift is used to move
cursor towards 5' direction to find the start of fragment
three_shift_s.append(d*3/4) # three shift is used to
move cursor towards 3' direction to find the end of fragment
else:
five_shift_s.append(0)
three_shift_s.append(d)
else:
# both sides
if halfextension:
five_shift_s.append(d/4)
three_shift_s.append(d/4)
else:
five_shift_s.append(d/2)
three_shift_s.append(d - d/2)

    prev_pileup = None

    for i in range(len(ds)):
        five_shift = five_shift_s[i]
        three_shift = three_shift_s[i]
        scale_factor = scale_factor_s[i]

        tmp_pileup = se_all_in_one_pileup ( self.__locations[chrom][0],

self.__locations[chrom][1], five_shift, three_shift, rlength, scale_factor,
baseline_value )

        if prev_pileup:
            prev_pileup = max_over_two_pv_array ( prev_pileup,

tmp_pileup )
else:
prev_pileup = tmp_pileup
return prev_pileup

When direction is false, we expand tags equally both five end side and
three end side equally. But when it is true then we give more priority to
three end.
Why did we need it?
What is reason for doing so?

Regards.

On Sun, Jun 22, 2014 at 3:12 PM, Kuldeep Yadav [email protected]
wrote:

Hello Tao, I want to translate MACS python code to java. So I am reading
the python code and trying to understand algorithm. And I have questions on
some points.

1. In functions __call_peak_wo_control() in file cPeakDetect.pyx, why did
   you took ‘d’ to be zero for paired end mode.

if self.PE_MODE: d = 0

P.S. please send me reading material if you have anything related to the
algorithm beside the MACS research paper. But that paper doesn’t explains
the implementation subtilities and also t does not have pseudo code.

Thank you.

On 20-Jun-2014, at 10:54 pm, Kuldeep Yadav [email protected] wrote:

I have already installed all prerequisite like cython, gcc etc. I will try
setup_w_cython.py.
Thank you.

On 20-Jun-2014, at 10:45 pm, Tao Liu (τν) [email protected]
wrote:

You need to install ‘cython’ in order to ‘translate’ pyx files to ‘c’
codes. Then you need ’setup_w_cython.py’ script to compile and install your
new version. You can find ’setup_w_cython.py’ on MACS2 github site.

Tao

On Jun 20, 2014, at 1:13 PM, Kuldeep Yadav [email protected]
wrote:

Thank You.

One more problem I am facing. How to debug .pyx code?

I want to print statements to check values of variables in .pyx files.
So I
made the changes to files and reinstalled macs with following command:
python setup.py install

but when I run it, I don't see any prints in terminal output. That is,
no
changes in console/terminal logs.
Am I missing something in installing steps?

P.S. I made changes in init() of cPeakDetect.pyx file. There is no
way
that this function is skipped while running.

On Fri, Jun 20, 2014 at 10:24 PM, Tao Liu (τν) [email protected]

wrote:

If you installed it through pip, then to call the "macs2" in current
working folder as "./macs2" should not be used, since pip should have
put
executable file in system path. You may installed a new MACS2 but
invoked
an older "macs2" script. Note, shiftsize option has been replaced by
extsize optioning the latest version.

Tao

On Jun 20, 2014, at 12:49 PM, Kuldeep Yadav <
[email protected]>
wrote:

I will be able to send that by monday, but I do remember version
It is latest 20140616
On Jun 20, 2014 10:12 PM, "Tao Liu (¦Ó¦Í)" [email protected]

wrote:

Could you post the entire error message? Also, could you provide
the
version of macs? Thanks!

On Jun 20, 2014, at 5:25 AM, Kuldeep Yadav <
[email protected]>

wrote:

While running following
./macs2 callpeak -t ~/macs\ input/something.sam

it shows error
'Namespace' object has no attribute 'shiftsize'

¡ª
Reply to this email directly or view it on GitHub.

¡ª
Reply to this email directly or view it on GitHub
#36 (comment).

¡ª
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

—
Reply to this email directly or view it on GitHub.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

While dealing with Paired-end data, 'd' is not used to extend sequencing reads, so it's set as 0. Direction for extension is different for ChIP and control. ChIP reads need to be extended towards 3' direction, while control reads only represent 'genomic bias' at the sonication cutting or digestion sites.

from macs.

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) [email protected]
wrote:

While dealing with Paired-end data, 'd' is not used to extend sequencing
reads, so it's set as 0. Direction for extension is different for ChIP and
control. ChIP reads need to be extended towards 3' direction, while control
reads only represent 'genomic bias' at the sonication cutting or digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

Hello Taoliu, I want to know whether master branch code for MACS is stable
or not? If not then which is latest stable version.

On Wed, Jul 9, 2014 at 10:57 AM, Kuldeep Yadav [email protected]
wrote:

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) [email protected]
wrote:

While dealing with Paired-end data, 'd' is not used to extend sequencing
reads, so it's set as 0. Direction for extension is different for ChIP and
control. ChIP reads need to be extended towards 3' direction, while control
reads only represent 'genomic bias' at the sonication cutting or digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

Hi Kuldeep,

Check: https://github.com/taoliu/MACS/tags

Use the tag for the recent release: 2.1.0.20140616

Best,
Tao

On Jul 10, 2014, at 8:47 PM, Kuldeep Yadav [email protected] wrote:

Hello Taoliu, I want to know whether master branch code for MACS is stable
or not? If not then which is latest stable version.

On Wed, Jul 9, 2014 at 10:57 AM, Kuldeep Yadav [email protected]
wrote:

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) [email protected]
wrote:

While dealing with Paired-end data, 'd' is not used to extend sequencing
reads, so it's set as 0. Direction for extension is different for ChIP and
control. ChIP reads need to be extended towards 3' direction, while control
reads only represent 'genomic bias' at the sonication cutting or digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

Thank you Tao,

There is a difference in code snippet for finding 'd' by correlation in
version 2.0.10.20131216 and 2.1.0.20140616.

In 2.0.10, the value in xcorr was taken as 'd' where ycorr was maximum. And
all local maximas were alternative 'd'.
What new thing you are doing in 2.1.0 ?

2.0.10 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2, num=len(ycorr))

    # smooth correlation values to get rid of local maximums from small

fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by default 11.

    # best cross-correlation point
    self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]

#+self.tag_expansion_size
# all local maximums could be alternative ds.
i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] & np.r_[ycorr[:-1] >
ycorr[1:], False]
self.alternative_d = map(int,xcorr[i_l_max])
# get rid of the last local maximum if it's at the right end of
curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

2.1.0 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2, num=len(ycorr))

    # smooth correlation values to get rid of local maximums from small

fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by default 11.

    # all local maximums could be alternative ds.
    i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] & np.r_[ycorr[:-1] >

ycorr[1:], False]
tmp_cor_alternative_d = ycorr[ i_l_max ]
tmp_alternative_d = xcorr[ i_l_max ]
cor_alternative_d = tmp_cor_alternative_d [ tmp_alternative_d > 0 ]
self.alternative_d = map( int, tmp_alternative_d[ tmp_alternative_d

0 ] )

    # best cross-correlation point
    self.d = xcorr[ np.where( ycorr== max( cor_alternative_d ) )[0][0] ]
    #self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]

#+self.tag_expansion_size

    # get rid of the last local maximum if it's at the right end of

curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

On Fri, Jul 11, 2014 at 7:12 PM, Tao Liu (τν) [email protected]
wrote:

Hi Kuldeep,

Check: https://github.com/taoliu/MACS/tags

Use the tag for the recent release: 2.1.0.20140616

Best,
Tao

On Jul 10, 2014, at 8:47 PM, Kuldeep Yadav [email protected]
wrote:

Hello Taoliu, I want to know whether master branch code for MACS is
stable
or not? If not then which is latest stable version.

On Wed, Jul 9, 2014 at 10:57 AM, Kuldeep Yadav [email protected]
wrote:

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) [email protected]

wrote:

While dealing with Paired-end data, 'd' is not used to extend
sequencing
reads, so it's set as 0. Direction for extension is different for
ChIP and
control. ChIP reads need to be extended towards 3' direction, while
control
reads only represent 'genomic bias' at the sonication cutting or
digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

Tao, One more doubt.

In cPeakModel.pyx we have a function __find_pair_center ( pluspeaks,
minuspeaks )

It find pairs of peaks (and center) such that distance between them is no
more than peaksize. Also a naive peak may be in more than one pair. For
example, I have a positive naive peak at 20 and two negative peak at 60 and
80. So I will have two centers that are 40 and 50 for pairs (20,60) and
(20,80).
Is it Correct?

code snippet:
cdef __find_pair_center (self, pluspeaks, minuspeaks):
ip = 0 # index for plus peaks
im = 0 # index for minus peaks
im_prev = 0 # index for minus peaks in previous plus
peak
pair_centers = array(BYTE4,[])
ip_max = len(pluspeaks)
im_max = len(minuspeaks)
flag_find_overlap = False
while ip<ip_max and im<im_max:
(pp,pn) = pluspeaks[ip] # for (peakposition, tagnumber in peak)
(mp,mn) = minuspeaks[im]
if pp-self.peaksize > mp: # move minus
im += 1
elif pp+self.peaksize < mp: # move plus
ip += 1
im = im_prev # search minus peaks from previous index
flag_find_overlap = False
else: # overlap!
if not flag_find_overlap:
flag_find_overlap = True
im_prev = im # only the first index is recorded
if float(pn)/mn < 2 and float(pn)/mn > 0.5: # number tags
in plus and minus peak region are comparable...
if pp < mp:
pair_centers.append((pp+mp)/2)
#self.debug ( "distance: %d, minus: %d, plus: %d" %
(mp-pp,mp,pp))
im += 1
return pair_centers

On Sat, Jul 12, 2014 at 9:10 AM, Kuldeep Yadav [email protected]
wrote:

Thank you Tao,

There is a difference in code snippet for finding 'd' by correlation in
version 2.0.10.20131216 and 2.1.0.20140616.

In 2.0.10, the value in xcorr was taken as 'd' where ycorr was maximum.
And all local maximas were alternative 'd'.
What new thing you are doing in 2.1.0 ?

2.0.10 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

    # smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by default
11.

    # best cross-correlation point
    self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]

#+self.tag_expansion_size
# all local maximums could be alternative ds.
i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] & np.r_[ycorr[:-1]

ycorr[1:], False]
self.alternative_d = map(int,xcorr[i_l_max])
# get rid of the last local maximum if it's at the right end of
curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

2.1.0 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

    # smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by default
11.

    # all local maximums could be alternative ds.
    i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] & np.r_[ycorr[:-1]

ycorr[1:], False]
tmp_cor_alternative_d = ycorr[ i_l_max ]
tmp_alternative_d = xcorr[ i_l_max ]
cor_alternative_d = tmp_cor_alternative_d [ tmp_alternative_d > 0
]
self.alternative_d = map( int, tmp_alternative_d[
tmp_alternative_d > 0 ] )

    # best cross-correlation point
    self.d = xcorr[ np.where( ycorr== max( cor_alternative_d ) )[0][0]

]
#self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]
#+self.tag_expansion_size

    # get rid of the last local maximum if it's at the right end of

curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

On Fri, Jul 11, 2014 at 7:12 PM, Tao Liu (τν) [email protected]
wrote:

Hi Kuldeep,

Check: https://github.com/taoliu/MACS/tags

Use the tag for the recent release: 2.1.0.20140616

Best,
Tao

On Jul 10, 2014, at 8:47 PM, Kuldeep Yadav [email protected]
wrote:

Hello Taoliu, I want to know whether master branch code for MACS is
stable
or not? If not then which is latest stable version.

On Wed, Jul 9, 2014 at 10:57 AM, Kuldeep Yadav [email protected]
wrote:

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) <
[email protected]>
wrote:

While dealing with Paired-end data, 'd' is not used to extend
sequencing
reads, so it's set as 0. Direction for extension is different for
ChIP and
control. ChIP reads need to be extended towards 3' direction, while
control
reads only represent 'genomic bias' at the sonication cutting or
digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

Hii Tao,

Have you used "alternative_d" anywhere? I traversed code and found that you
computed it but did not use anywhere.

On Sat, Jul 12, 2014 at 9:32 AM, Kuldeep Yadav [email protected]
wrote:

Tao, One more doubt.

In cPeakModel.pyx we have a function __find_pair_center ( pluspeaks,
minuspeaks )

It find pairs of peaks (and center) such that distance between them is no
more than peaksize. Also a naive peak may be in more than one pair. For
example, I have a positive naive peak at 20 and two negative peak at 60 and
80. So I will have two centers that are 40 and 50 for pairs (20,60) and
(20,80).
Is it Correct?

code snippet:
cdef __find_pair_center (self, pluspeaks, minuspeaks):
ip = 0 # index for plus peaks
im = 0 # index for minus peaks
im_prev = 0 # index for minus peaks in previous plus
peak
pair_centers = array(BYTE4,[])
ip_max = len(pluspeaks)
im_max = len(minuspeaks)
flag_find_overlap = False
while ip<ip_max and im<im_max:
(pp,pn) = pluspeaks[ip] # for (peakposition, tagnumber in peak)
(mp,mn) = minuspeaks[im]
if pp-self.peaksize > mp: # move minus
im += 1
elif pp+self.peaksize < mp: # move plus
ip += 1
im = im_prev # search minus peaks from previous index
flag_find_overlap = False
else: # overlap!
if not flag_find_overlap:
flag_find_overlap = True
im_prev = im # only the first index is recorded
if float(pn)/mn < 2 and float(pn)/mn > 0.5: # number tags
in plus and minus peak region are comparable...
if pp < mp:
pair_centers.append((pp+mp)/2)
#self.debug ( "distance: %d, minus: %d, plus: %d"
% (mp-pp,mp,pp))
im += 1
return pair_centers

On Sat, Jul 12, 2014 at 9:10 AM, Kuldeep Yadav [email protected]
wrote:

Thank you Tao,

There is a difference in code snippet for finding 'd' by correlation in
version 2.0.10.20131216 and 2.1.0.20140616.

In 2.0.10, the value in xcorr was taken as 'd' where ycorr was maximum.
And all local maximas were alternative 'd'.
What new thing you are doing in 2.1.0 ?

2.0.10 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

    # smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by default
11.

    # best cross-correlation point
    self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]

#+self.tag_expansion_size
# all local maximums could be alternative ds.
i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] & np.r_[ycorr[:-1]

ycorr[1:], False]
self.alternative_d = map(int,xcorr[i_l_max])
# get rid of the last local maximum if it's at the right end of
curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

2.1.0 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

    # smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by default
11.

    # all local maximums could be alternative ds.
    i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] & np.r_[ycorr[:-1]

ycorr[1:], False]
tmp_cor_alternative_d = ycorr[ i_l_max ]
tmp_alternative_d = xcorr[ i_l_max ]
cor_alternative_d = tmp_cor_alternative_d [ tmp_alternative_d >
0 ]
self.alternative_d = map( int, tmp_alternative_d[
tmp_alternative_d > 0 ] )

    # best cross-correlation point
    self.d = xcorr[ np.where( ycorr== max( cor_alternative_d )

)[0][0] ]
#self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]
#+self.tag_expansion_size

    # get rid of the last local maximum if it's at the right end of

curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

On Fri, Jul 11, 2014 at 7:12 PM, Tao Liu (τν) [email protected]
wrote:

Hi Kuldeep,

Check: https://github.com/taoliu/MACS/tags

Use the tag for the recent release: 2.1.0.20140616

Best,
Tao

On Jul 10, 2014, at 8:47 PM, Kuldeep Yadav [email protected]
wrote:

Hello Taoliu, I want to know whether master branch code for MACS is
stable
or not? If not then which is latest stable version.

On Wed, Jul 9, 2014 at 10:57 AM, Kuldeep Yadav [email protected]

wrote:

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) <
[email protected]>
wrote:

While dealing with Paired-end data, 'd' is not used to extend
sequencing
reads, so it's set as 0. Direction for extension is different for
ChIP and
control. ChIP reads need to be extended towards 3' direction, while
control
reads only represent 'genomic bias' at the sonication cutting or
digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

Hii Tao,
I figured out above doubts myself. But now I have a new :) . That is, in
__model_add_line() why did you use psize_adjusted instead of peaksize to
select tags located near by a peak center. It should have been peaksize
keeping window size in consideration.

cdef __model_add_line (self, object pos1, np.ndarray pos2, np.ndarray

start, np.ndarray end): #, int plus_strand=1):
"""Project each pos in pos2 which is included in
[pos1-self.peaksize,pos1+self.peaksize] to the line.

    pos1: paired centers -- array.array
    pos2: tags of certain strand -- a numpy.array object
    line: numpy array object where we pileup tags

    """
    cdef int i1, i2, i2_prev, i1_max, i2_max, last_p2, psize_adjusted1,

psize_adjusted2, p1, p2, max_index, s, e

    i1 = 0                  # index for pos1
    i2 = 0                  # index for pos2
    i2_prev = 0             # index for pos2 in previous pos1
                            # [pos1-self.peaksize,pos1+self.peaksize]
                            # region
    i1_max = len(pos1)
    i2_max = pos2.shape[0]
    last_p2 = -1
    flag_find_overlap = False

    max_index = start.shape[0] - 1

    psize_adjusted1 = self.peaksize + self.tag_expansion_size / 2 #

half window

    while i1<i1_max and i2<i2_max:
        p1 = pos1[i1]
        #if plus_strand:
        #    p2 = pos2[i2]
        #else:
        #    p2 = pos2[i2] - self.tag_expansion_size

        p2 = pos2[i2] #- self.tag_expansion_size/2

        if p1-psize_adjusted1 > p2: # move pos2
            i2 += 1
        elif p1+psize_adjusted1 < p2: # move pos1
            i1 += 1
            i2 = i2_prev    # search minus peaks from previous index
            flag_find_overlap = False
        else:               # overlap!
            if not flag_find_overlap:
                flag_find_overlap = True
                i2_prev = i2 # only the first index is recorded
            # project
            #for i in

range(p2-p1+self.peaksize,p2-p1+self.peaksize+self.tag_expansion_size):
s = max(p2-self.tag_expansion_size/2-p1+psize_adjusted1, 0)
start[s] += 1
e = min(p2+self.tag_expansion_size/2-p1+psize_adjusted1,
max_index)
end[e] -= 1
#line[s:e] += 1
#for i in range(s,e):
# #if i>=0 and i<length_l:
# line[i]+=1
i2+=1
return

On Mon, Jul 14, 2014 at 8:29 PM, Kuldeep Yadav [email protected]
wrote:

Hii Tao,

Have you used "alternative_d" anywhere? I traversed code and found that
you computed it but did not use anywhere.

On Sat, Jul 12, 2014 at 9:32 AM, Kuldeep Yadav [email protected]
wrote:

Tao, One more doubt.

In cPeakModel.pyx we have a function __find_pair_center ( pluspeaks,
minuspeaks )

It find pairs of peaks (and center) such that distance between them is no
more than peaksize. Also a naive peak may be in more than one pair. For
example, I have a positive naive peak at 20 and two negative peak at 60 and
80. So I will have two centers that are 40 and 50 for pairs (20,60) and
(20,80).
Is it Correct?

code snippet:
cdef __find_pair_center (self, pluspeaks, minuspeaks):
ip = 0 # index for plus peaks
im = 0 # index for minus peaks
im_prev = 0 # index for minus peaks in previous plus
peak
pair_centers = array(BYTE4,[])
ip_max = len(pluspeaks)
im_max = len(minuspeaks)
flag_find_overlap = False
while ip<ip_max and im<im_max:
(pp,pn) = pluspeaks[ip] # for (peakposition, tagnumber in
peak)
(mp,mn) = minuspeaks[im]
if pp-self.peaksize > mp: # move minus
im += 1
elif pp+self.peaksize < mp: # move plus
ip += 1
im = im_prev # search minus peaks from previous index
flag_find_overlap = False
else: # overlap!
if not flag_find_overlap:
flag_find_overlap = True
im_prev = im # only the first index is recorded
if float(pn)/mn < 2 and float(pn)/mn > 0.5: # number tags
in plus and minus peak region are comparable...
if pp < mp:
pair_centers.append((pp+mp)/2)
#self.debug ( "distance: %d, minus: %d, plus: %d"
% (mp-pp,mp,pp))
im += 1
return pair_centers

On Sat, Jul 12, 2014 at 9:10 AM, Kuldeep Yadav [email protected]
wrote:

Thank you Tao,

There is a difference in code snippet for finding 'd' by correlation in
version 2.0.10.20131216 and 2.1.0.20140616.

In 2.0.10, the value in xcorr was taken as 'd' where ycorr was maximum.
And all local maximas were alternative 'd'.
What new thing you are doing in 2.1.0 ?

2.0.10 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

    # smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by default
11.

    # best cross-correlation point
    self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]

#+self.tag_expansion_size
# all local maximums could be alternative ds.
i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] &
np.r_[ycorr[:-1] > ycorr[1:], False]
self.alternative_d = map(int,xcorr[i_l_max])
# get rid of the last local maximum if it's at the right end of
curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

2.1.0 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

    # smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by default
11.

    # all local maximums could be alternative ds.
    i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] &

np.r_[ycorr[:-1] > ycorr[1:], False]
tmp_cor_alternative_d = ycorr[ i_l_max ]
tmp_alternative_d = xcorr[ i_l_max ]
cor_alternative_d = tmp_cor_alternative_d [ tmp_alternative_d >
0 ]
self.alternative_d = map( int, tmp_alternative_d[
tmp_alternative_d > 0 ] )

    # best cross-correlation point
    self.d = xcorr[ np.where( ycorr== max( cor_alternative_d )

)[0][0] ]
#self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]
#+self.tag_expansion_size

    # get rid of the last local maximum if it's at the right end of

curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

On Fri, Jul 11, 2014 at 7:12 PM, Tao Liu (τν) [email protected]
wrote:

Hi Kuldeep,

Check: https://github.com/taoliu/MACS/tags

Use the tag for the recent release: 2.1.0.20140616

Best,
Tao

On Jul 10, 2014, at 8:47 PM, Kuldeep Yadav [email protected]
wrote:

Hello Taoliu, I want to know whether master branch code for MACS is
stable
or not? If not then which is latest stable version.

On Wed, Jul 9, 2014 at 10:57 AM, Kuldeep Yadav [email protected]

wrote:

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) <
[email protected]>
wrote:

While dealing with Paired-end data, 'd' is not used to extend
sequencing
reads, so it's set as 0. Direction for extension is different for
ChIP and
control. ChIP reads need to be extended towards 3' direction,
while control
reads only represent 'genomic bias' at the sonication cutting or
digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

Hello Tao,
Please help.

I am in final phase of work and that is computing all parameters like
lambda, scale factors etc.
Why did You need to multiply count of total control tags by 2.

https://github.com/taoliu/MACS/blob/master/MACS2/cPeakDetect.pyx#L149

Can't I do
ratio_treat2control = (treat.length/treat.average_template_length) /
(control.length/control.average_template_length)
This is actually ratio of counts of treatment tags to control tags.

On Tue, Jul 22, 2014 at 4:55 PM, Kuldeep Yadav [email protected]
wrote:

Hii Tao,
I figured out above doubts myself. But now I have a new :) . That is, in
__model_add_line() why did you use psize_adjusted instead of peaksize to
select tags located near by a peak center. It should have been peaksize
keeping window size in consideration.

cdef __model_add_line (self, object pos1, np.ndarray pos2, np.ndarray

start, np.ndarray end): #, int plus_strand=1):
"""Project each pos in pos2 which is included in
[pos1-self.peaksize,pos1+self.peaksize] to the line.

    pos1: paired centers -- array.array
    pos2: tags of certain strand -- a numpy.array object
    line: numpy array object where we pileup tags

    """
    cdef int i1, i2, i2_prev, i1_max, i2_max, last_p2,

psize_adjusted1, psize_adjusted2, p1, p2, max_index, s, e

    i1 = 0                  # index for pos1
    i2 = 0                  # index for pos2
    i2_prev = 0             # index for pos2 in previous pos1
                            # [pos1-self.peaksize,pos1+self.peaksize]
                            # region
    i1_max = len(pos1)
    i2_max = pos2.shape[0]
    last_p2 = -1
    flag_find_overlap = False

    max_index = start.shape[0] - 1

    psize_adjusted1 = self.peaksize + self.tag_expansion_size / 2 #

half window

    while i1<i1_max and i2<i2_max:
        p1 = pos1[i1]
        #if plus_strand:
        #    p2 = pos2[i2]
        #else:
        #    p2 = pos2[i2] - self.tag_expansion_size

        p2 = pos2[i2] #- self.tag_expansion_size/2

        if p1-psize_adjusted1 > p2: # move pos2
            i2 += 1
        elif p1+psize_adjusted1 < p2: # move pos1
            i1 += 1
            i2 = i2_prev    # search minus peaks from previous index
            flag_find_overlap = False
        else:               # overlap!
            if not flag_find_overlap:
                flag_find_overlap = True
                i2_prev = i2 # only the first index is recorded
            # project
            #for i in

range(p2-p1+self.peaksize,p2-p1+self.peaksize+self.tag_expansion_size):
s = max(p2-self.tag_expansion_size/2-p1+psize_adjusted1, 0)
start[s] += 1
e = min(p2+self.tag_expansion_size/2-p1+psize_adjusted1,
max_index)
end[e] -= 1
#line[s:e] += 1
#for i in range(s,e):
# #if i>=0 and i<length_l:
# line[i]+=1
i2+=1
return

On Mon, Jul 14, 2014 at 8:29 PM, Kuldeep Yadav [email protected]
wrote:

Hii Tao,

Have you used "alternative_d" anywhere? I traversed code and found that
you computed it but did not use anywhere.

On Sat, Jul 12, 2014 at 9:32 AM, Kuldeep Yadav [email protected]
wrote:

Tao, One more doubt.

In cPeakModel.pyx we have a function __find_pair_center ( pluspeaks,
minuspeaks )

It find pairs of peaks (and center) such that distance between them is
no more than peaksize. Also a naive peak may be in more than one pair. For
example, I have a positive naive peak at 20 and two negative peak at 60 and
80. So I will have two centers that are 40 and 50 for pairs (20,60) and
(20,80).
Is it Correct?

code snippet:
cdef __find_pair_center (self, pluspeaks, minuspeaks):
ip = 0 # index for plus peaks
im = 0 # index for minus peaks
im_prev = 0 # index for minus peaks in previous plus
peak
pair_centers = array(BYTE4,[])
ip_max = len(pluspeaks)
im_max = len(minuspeaks)
flag_find_overlap = False
while ip<ip_max and im<im_max:
(pp,pn) = pluspeaks[ip] # for (peakposition, tagnumber in
peak)
(mp,mn) = minuspeaks[im]
if pp-self.peaksize > mp: # move minus
im += 1
elif pp+self.peaksize < mp: # move plus
ip += 1
im = im_prev # search minus peaks from previous index
flag_find_overlap = False
else: # overlap!
if not flag_find_overlap:
flag_find_overlap = True
im_prev = im # only the first index is recorded
if float(pn)/mn < 2 and float(pn)/mn > 0.5: # number
tags in plus and minus peak region are comparable...
if pp < mp:
pair_centers.append((pp+mp)/2)
#self.debug ( "distance: %d, minus: %d, plus:
%d" % (mp-pp,mp,pp))
im += 1
return pair_centers

On Sat, Jul 12, 2014 at 9:10 AM, Kuldeep Yadav [email protected]
wrote:

Thank you Tao,

There is a difference in code snippet for finding 'd' by correlation in
version 2.0.10.20131216 and 2.1.0.20140616.

In 2.0.10, the value in xcorr was taken as 'd' where ycorr was maximum.
And all local maximas were alternative 'd'.
What new thing you are doing in 2.1.0 ?

2.0.10 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

    # smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by
default 11.

    # best cross-correlation point
    self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]

#+self.tag_expansion_size
# all local maximums could be alternative ds.
i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] &
np.r_[ycorr[:-1] > ycorr[1:], False]
self.alternative_d = map(int,xcorr[i_l_max])
# get rid of the last local maximum if it's at the right end of
curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

2.1.0 code
# cross-correlation
ycorr =
np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]
xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

    # smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by
default 11.

    # all local maximums could be alternative ds.
    i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] &

np.r_[ycorr[:-1] > ycorr[1:], False]
tmp_cor_alternative_d = ycorr[ i_l_max ]
tmp_alternative_d = xcorr[ i_l_max ]
cor_alternative_d = tmp_cor_alternative_d [ tmp_alternative_d

0 ]
self.alternative_d = map( int, tmp_alternative_d[
tmp_alternative_d > 0 ] )

    # best cross-correlation point
    self.d = xcorr[ np.where( ycorr== max( cor_alternative_d )

)[0][0] ]
#self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]
#+self.tag_expansion_size

    # get rid of the last local maximum if it's at the right end of

curve.

    assert len(self.alternative_d) > 0, "No proper d can be found!

Tweak --mfold?"

    self.ycorr = ycorr
    self.xcorr = xcorr

On Fri, Jul 11, 2014 at 7:12 PM, Tao Liu (τν) <[email protected]

wrote:

Hi Kuldeep,

Check: https://github.com/taoliu/MACS/tags

Use the tag for the recent release: 2.1.0.20140616

Best,
Tao

On Jul 10, 2014, at 8:47 PM, Kuldeep Yadav [email protected]
wrote:

Hello Taoliu, I want to know whether master branch code for MACS is
stable
or not? If not then which is latest stable version.

On Wed, Jul 9, 2014 at 10:57 AM, Kuldeep Yadav <
[email protected]>
wrote:

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) <
[email protected]>
wrote:

While dealing with Paired-end data, 'd' is not used to extend
sequencing
reads, so it's set as 0. Direction for extension is different for
ChIP and
control. ChIP reads need to be extended towards 3' direction,
while control
reads only represent 'genomic bias' at the sonication cutting or
digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA

from macs.

That code is for pair-end data. I believe I put comments there. In
treatment, two ends represent a fragment/template, whereas in control, I
assume each control tag represent a 'fragment', because we want to capture
the bias from control at 'cutting sites'.

Tao

On Thu, Aug 21, 2014 at 8:33 AM, Kuldeep Yadav [email protected]
wrote:

Hello Tao,
Please help.

I am in final phase of work and that is computing all parameters like
lambda, scale factors etc.
Why did You need to multiply count of total control tags by 2.

https://github.com/taoliu/MACS/blob/master/MACS2/cPeakDetect.pyx#L149

Can't I do
ratio_treat2control = (treat.length/treat.average_template_length) /
(control.length/control.average_template_length)
This is actually ratio of counts of treatment tags to control tags.

On Tue, Jul 22, 2014 at 4:55 PM, Kuldeep Yadav [email protected]
wrote:

Hii Tao,
I figured out above doubts myself. But now I have a new :) . That is, in
__model_add_line() why did you use psize_adjusted instead of peaksize to
select tags located near by a peak center. It should have been peaksize
keeping window size in consideration.

cdef __model_add_line (self, object pos1, np.ndarray pos2, np.ndarray
start, np.ndarray end): #, int plus_strand=1):
"""Project each pos in pos2 which is included in
[pos1-self.peaksize,pos1+self.peaksize] to the line.

pos1: paired centers -- array.array
pos2: tags of certain strand -- a numpy.array object
line: numpy array object where we pileup tags

"""
cdef int i1, i2, i2_prev, i1_max, i2_max, last_p2,
psize_adjusted1, psize_adjusted2, p1, p2, max_index, s, e

i1 = 0 # index for pos1
i2 = 0 # index for pos2
i2_prev = 0 # index for pos2 in previous pos1

[pos1-self.peaksize,pos1+self.peaksize]

region

i1_max = len(pos1)
i2_max = pos2.shape[0]
last_p2 = -1
flag_find_overlap = False

max_index = start.shape[0] - 1

psize_adjusted1 = self.peaksize + self.tag_expansion_size / 2 #
half window

while i1<i1_max and i2<i2_max:
p1 = pos1[i1]
#if plus_strand:

p2 = pos2[i2]

#else:

p2 = pos2[i2] - self.tag_expansion_size

p2 = pos2[i2] #- self.tag_expansion_size/2

if p1-psize_adjusted1 > p2: # move pos2
i2 += 1
elif p1+psize_adjusted1 < p2: # move pos1
i1 += 1
i2 = i2_prev # search minus peaks from previous index
flag_find_overlap = False
else: # overlap!
if not flag_find_overlap:
flag_find_overlap = True
i2_prev = i2 # only the first index is recorded

project

#for i in
range(p2-p1+self.peaksize,p2-p1+self.peaksize+self.tag_expansion_size):
s = max(p2-self.tag_expansion_size/2-p1+psize_adjusted1, 0)
start[s] += 1
e = min(p2+self.tag_expansion_size/2-p1+psize_adjusted1,
max_index)
end[e] -= 1
#line[s:e] += 1
#for i in range(s,e):

#if i>=0 and i<length_l:

line[i]+=1

i2+=1
return

On Mon, Jul 14, 2014 at 8:29 PM, Kuldeep Yadav [email protected]
wrote:

Hii Tao,

Have you used "alternative_d" anywhere? I traversed code and found that
you computed it but did not use anywhere.

On Sat, Jul 12, 2014 at 9:32 AM, Kuldeep Yadav [email protected]
wrote:

Tao, One more doubt.

In cPeakModel.pyx we have a function __find_pair_center ( pluspeaks,
minuspeaks )

It find pairs of peaks (and center) such that distance between them is
no more than peaksize. Also a naive peak may be in more than one pair.
For
example, I have a positive naive peak at 20 and two negative peak at
60 and
80. So I will have two centers that are 40 and 50 for pairs (20,60)
and
(20,80).
Is it Correct?

code snippet:
cdef __find_pair_center (self, pluspeaks, minuspeaks):
ip = 0 # index for plus peaks
im = 0 # index for minus peaks
im_prev = 0 # index for minus peaks in previous plus
peak
pair_centers = array(BYTE4,[])
ip_max = len(pluspeaks)
im_max = len(minuspeaks)
flag_find_overlap = False
while ip<ip_max and im<im_max:
(pp,pn) = pluspeaks[ip] # for (peakposition, tagnumber in
peak)
(mp,mn) = minuspeaks[im]
if pp-self.peaksize > mp: # move minus
im += 1
elif pp+self.peaksize < mp: # move plus
ip += 1
im = im_prev # search minus peaks from previous index
flag_find_overlap = False
else: # overlap!
if not flag_find_overlap:
flag_find_overlap = True
im_prev = im # only the first index is recorded
if float(pn)/mn < 2 and float(pn)/mn > 0.5: # number
tags in plus and minus peak region are comparable...
if pp < mp:
pair_centers.append((pp+mp)/2)
#self.debug ( "distance: %d, minus: %d, plus:
%d" % (mp-pp,mp,pp))
im += 1
return pair_centers

On Sat, Jul 12, 2014 at 9:10 AM, Kuldeep Yadav [email protected]

wrote:

Thank you Tao,

There is a difference in code snippet for finding 'd' by correlation
in
version 2.0.10.20131216 and 2.1.0.20140616.

In 2.0.10, the value in xcorr was taken as 'd' where ycorr was
maximum.
And all local maximas were alternative 'd'.
What new thing you are doing in 2.1.0 ?

2.0.10 code

cross-correlation

ycorr =

np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]

xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by
default 11.

best cross-correlation point

self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]
#+self.tag_expansion_size

all local maximums could be alternative ds.

i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] &
np.r_[ycorr[:-1] > ycorr[1:], False]
self.alternative_d = map(int,xcorr[i_l_max])

get rid of the last local maximum if it's at the right end of

curve.

assert len(self.alternative_d) > 0, "No proper d can be found!
Tweak --mfold?"

self.ycorr = ycorr
self.xcorr = xcorr

2.1.0 code

cross-correlation

ycorr =

np.correlate(minus_data,plus_data,mode="full")[window_size-self.peaksize:window_size+self.peaksize]

xcorr = np.linspace(len(ycorr)//2*-1, len(ycorr)//2,
num=len(ycorr))

smooth correlation values to get rid of local maximums from

small fluctuations.
ycorr = smooth(ycorr, window="flat") # window size is by
default 11.

all local maximums could be alternative ds.

i_l_max = np.r_[False, ycorr[1:] > ycorr[:-1]] &
np.r_[ycorr[:-1] > ycorr[1:], False]
tmp_cor_alternative_d = ycorr[ i_l_max ]
tmp_alternative_d = xcorr[ i_l_max ]
cor_alternative_d = tmp_cor_alternative_d [ tmp_alternative_d

0 ]
self.alternative_d = map( int, tmp_alternative_d[
tmp_alternative_d > 0 ] )

best cross-correlation point

self.d = xcorr[ np.where( ycorr== max( cor_alternative_d )
)[0][0] ]
#self.d = xcorr[np.where(ycorr==max(ycorr))[0][0]]
#+self.tag_expansion_size

get rid of the last local maximum if it's at the right end of

curve.

assert len(self.alternative_d) > 0, "No proper d can be found!
Tweak --mfold?"

self.ycorr = ycorr
self.xcorr = xcorr

On Fri, Jul 11, 2014 at 7:12 PM, Tao Liu (τν) <
[email protected]

wrote:

Hi Kuldeep,

Check: https://github.com/taoliu/MACS/tags

Use the tag for the recent release: 2.1.0.20140616

Best,
Tao

On Jul 10, 2014, at 8:47 PM, Kuldeep Yadav [email protected]

wrote:

Hello Taoliu, I want to know whether master branch code for MACS
is
stable
or not? If not then which is latest stable version.

On Wed, Jul 9, 2014 at 10:57 AM, Kuldeep Yadav <
[email protected]>
wrote:

Thank you. That certainly helped.

On Wed, Jul 9, 2014 at 1:09 AM, Tao Liu (τν) <
[email protected]>
wrote:

While dealing with Paired-end data, 'd' is not used to extend
sequencing
reads, so it's set as 0. Direction for extension is different
for
ChIP and
control. ChIP reads need to be extended towards 3' direction,
while control
reads only represent 'genomic bias' at the sonication cutting
or
digestion
sites.

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
*Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<

http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>*

—
Reply to this email directly or view it on GitHub
#36 (comment).

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

Regards.

KULDEEP
+919535580366
Proud to be a donor of the AAP. Join the movement
today. http://donate.aamaadmiparty.org
<
http://messenger.aamaadmiparty.org/l/MxCuqEumiVEr9ClU9oqxZA/9KS50hXPn6DRjLRWMOJzkw/M4r8D763LgDtF7f16PFYqQeA>

—

Reply to this email directly or view it on GitHub
#36 (comment).

from macs.