Call peaks on modified pileup file, instead of calling on modified BAMPE file

[jeanmonet]

Hi,

The initial peak calling is performed with the following command:

macs3 callpeak -f BAMPE -t $bam -g hs -n $name -B -q 0.01 --outdir $outdir

After some filtering, I'm modifying the sample_treat_pileup.bdg pileup file (and the fragments file data), and would like to perform peak calling form this intermediary step (the filtered pileup file), instead of filtering the BAMPE file and and performing peak calling on that filtered BAMPE (which would re-build the pile-up, etc.).

I am aware of the bdgpeakcall command, however it does not perform the model building and peak calling as per the pipeline in callpeak, so peaks obtained are not the same: ie. performing peak calling on the initial pileup file results in different peaks.

Any other way to achieve this?

[taoliu]

@jeanmonet Check this tutorial on how to do this: https://macs3-project.github.io/MACS/docs/Advanced_Step-by-step_Peak_Calling.html Since you don't have control data, you can use a single value for the whole genome background. In your case, the background = number_of_pairs * average_insertion_size / genome_size. The number_of_pairs and average_insertion_size should be found from your previous callpeak run. Then manually make a bedGraph file where for each chromosome and every position the value is the same. For example, the chromosome 1 has 123,456bps and we plan to set the background as 0.789, and the chromosome 2 has 101,112 bps:

chr1 0 123456 0.0789
chr2 0 101112 0.0789

Then follow the instruction of computing qscores with bdgcmp and the step to call peaks using bdgpeakcall.

[jeanmonet]

I was able to achieve this by constructing a PETrackI object using the fragments data (start & stop).
There are very slight differences in the pile-up resulting from processing the BAM file versus the pile-up obtained from the fragments data directly (by constructing the PETrackI).
Fragments data either comes from CellRanger or from conversion of BAM to fragments with sinto. The source BAM was filtered first for min MAPQ score to ensure same reads are present.
Somehow it seems fragments extracted from BAM by MACS3 are mostly the same but not exactly the same as those extracted with sinto from the same file. I could not yet determine why this is the case.

BAM --> MACS3 --> PETrackI --> pile-up
vs
BAM --> sinto --> fragments --> PETrackI --> pile-up