Finding Fragment Size chip-seq paired end reads

Question

[luckysardar171]

Hello currently I am working on Chip-seq analysis i have paired end data 35-130 PE reads I have aligned with bowtie2, after mapping generated the bam. Using MACS2 i wanted call the peaks(H3K4me3 histone mark) when I run run_spp.R for Modelling Fragment Size i got the below output result i am little bit confused which value should i use from the pdf generated and the tab generated file

command used

run_spp.R -c=output.bam -rf -out=params.out

plot https://ibb.co/85qBRWt

tab file output

"chipSampleMaster.tagAlign.gz 48643337 160,345 0.409606223030329,0.372288179379363 135 0.3858776 600 0.3589484 1.141129 1.881149 2"

macs2 callpeak -t chip.bam -n H3K4me3 --gsize hs -c input.bam --nomodel --shiftsize=?

Thank you

Answer 1

[ATpoint]

Use -f BAMPE and macs2 will infer the fragment size based on the TLEN field from the BAM file. I would also remove --nomodel --shiftsize, both are not options you would want in ChIP-seq. It is intended for non-ChIP assays such as ATAC-seq where you have to turn off the shifting model it uses by default. H3K4me3 is a narrow mark so just use the defaults and -f BAMPE for paired end data.

[luckysardar171]

Okay Thank you for your suggestions, can you please tell me usage of --SPMR in MACS2 because encode project uses in its pipeline

[ATpoint]

This option has nothing to do with peak calling, it only saves the pipeups so you can display it on a genome browser, but for this I usually use deeptools bamCoverage.

macs2 callpeak -t chip.bam -n H3K4me3 --gsize hs -c input.bam -f BAMPE

This is probably what you need, no magic here.

[luckysardar171]

Okay Thank you. for my experiment I have one wt and knockout data i wanted compare this two and show the H3K4me3 status how I have to proceed, do the differential analysis or what, any examples and papers will be better to understand

Thank you