Procedure to generate control pileup for Paired-end data

[Aarthi Mohan]

Hi Tao and all,

I am looking for the steps to generate the control pileup for a paired-end library. 

From 'callpeak' command, I can see that macs2 does not use the fragment size detected from the control paired-end reads

mean fragment size is determined as 273 bp from treatment
INFO  @ Tue, 09 Aug 2016 23:53:43: #1 note: mean fragment size in control is 255 bp -- value ignored
INFO  @ Tue, 09 Aug 2016 23:53:43: #1 fragment size = 273

I have followed the https://github.com/taoliu/MACS/wiki/Advanced:-Call-peaks-using-MACS2-subcommands tutorial for generating pileups with single-end. 

It has been straight forward to generate the ChIP (treat) paired-end pileup, since most of my experiments have more fragments in control. However, I am not able to generate the control pileup correctly (i did try using the above 273bp and following the advanced tutorial, but the final pileup values don't match the ones' generated by callpeak run with -B option). 

It will be very very helpful, if the steps to generate a control pileup with paired-end is shared. 

Thanks,

Aarthi

[Moshe Olshansky]

Hello Aarthi,

I think that if you have paired end reads in MACS2 BEDPE format you can use bedtools genomecov command to generate the pileup (instead of macs2 pileup).

I will be very indebted if Tao comments on this and says whether this is true.

Moshe.

[Ian Donaldson]

This may be of help, see Aaron Quinlan's answer.

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

[Ian Donaldson]

https://www.biostars.org/p/18258/

[Aarthi Mohan]

Hi Ian and Moshe,

Thanks for your responses. I have followed the same Biostar post to generate the 'treatment pileup' and it worked fine. But the control piluep (lambda) generation follows different procedure, and there is not much documentation on 'How the control pileup will be generated for a Paired-end data'. 

So, I tried to figure it out by going through the code, and I guess the following is roughly the steps to generate the control pileup:

For control library, treat the PE as SE tags.

1. Extract the left end and right end of the fragment as two reads: left read (+ve strand) and right read -ve strand) [check function  'pileup_a_chromosome_c' in CallPeakUnit.pyx]. So, due to this we have twice the amount of fragments in control, and this is confirmed from this line in PeakDetect.pyx

if self.PE_MODE:
            d = self.treat.average_template_length
            control_total = self.control.total * 2 # in PE mode, entire fragment is counted as 1
                                                   # in treatment whereas both ends of fragment are counted in control/input.

2. Extend each tag generated from step1 towards the 3' direction, by d/2 (d is derived as average fragment length from treatment). [If you goto 'pileup_a_chromosome_c' function, the function by default performs directional read extension (calling function does not set the direction=False flag). This is different from how the control pileup is generated with a SE library]

3. Extension is also done for slocal(1kb regions), llocal (10kb region). 

4. lambda_bg (whole genome_depth), which is the baseline value to be filled for no coverage regions is computed as 

lambda_bg = (sum of fragment length from treatment)  /genome_size [this is the case if the number of fragments in treat is lesser than control; Since the control fragments are counted twice the original size, the control size is always the higher library is my case. Otherwise the formula switches: refer line 147-173 in PeakDetect.pyx]

5. For each window/bp, substitute the max value from pileup generated from step 2-4. 

6. If SPMR is given in callpeaks command, scale the pileup values to adjust for sequencing depth difference [value from this variable 'effective_depth_in_million' is used]

I have tried many times to generate each step with macs2 pileup seperately (like shown in the tutorial here for SE libraries: https://github.com/taoliu/MACS/wiki/Advanced:-Call-peaks-using-MACS2-subcommands ), but the final pileup values are not same as one would generate through 'macs2 callpeak' command. 

It will be really helpful if anyone can confirm or a the amend the above "steps for generating  control pileup for a Paired-end (PE) library". 

Thanks and regards,

Aarthi

On Wednesday, August 24, 2016 at 4:52:23 PM UTC+8, Ian wrote:

https://www.biostars.org/p/18258/

On 24 August 2016 at 09:51, Ian Donaldson <donald...@gmail.com> wrote:

This may be of help, see Aaron Quinlan's answer.

On 24 August 2016 at 09:29, Moshe Olshansky <moshe.o...@gmail.com> wrote:

Hello Aarthi,

I think that if you have paired end reads in MACS2 BEDPE format you can use bedtools genomecov command to generate the pileup (instead of macs2 pileup).

I will be very indebted if Tao comments on this and says whether this is true.

Moshe.

On Wednesday, August 10, 2016 at 2:04:17 AM UTC+10, Aarthi Mohan wrote:

Hi Tao and all,

I am looking for the steps to generate the control pileup for a paired-end library. 

From 'callpeak' command, I can see that macs2 does not use the fragment size detected from the control paired-end reads

mean fragment size is determined as 273 bp from treatment
INFO  @ Tue, 09 Aug 2016 23:53:43: #1 note: mean fragment size in control is 255 bp -- value ignored
INFO  @ Tue, 09 Aug 2016 23:53:43: #1 fragment size = 273

I have followed the https://github.com/taoliu/MACS/wiki/Advanced:-Call-peaks-using-MACS2-subcommands tutorial for generating pileups with single-end. 

It has been straight forward to generate the ChIP (treat) paired-end pileup, since most of my experiments have more fragments in control. However, I am not able to generate the control pileup correctly (i did try using the above 273bp and following the advanced tutorial, but the final pileup values don't match the ones' generated by callpeak run with -B option). 

It will be very very helpful, if the steps to generate a control pileup with paired-end is shared. 

Thanks,

Aarthi

--
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.