question about calling peak by macs2

Question

[ch8316f5eyu]

I found the result of macs2 is not perfect and I want to know why. Firstly, I briefly describe the ChIP I made: The sequencing files are paired-end. Because of poor sequencing depth, I merge two input and two IP sample into one input and one sample. After processing the raw fastq files, I used bowtie2 to map these reads. Then I used macs2 to call peak.

macs2 callpeak -t IP_mm10.bam -c Input.bam --outdir . -n IP -g 2652783500 -q 0.1 -f BAMPE --keep-dup all -B

Because I am interested repeat region, I retaind all duplicated. I found the effective genome size from deeptools for mm10. I think the process is OK, but the result is weird. The input file also has peaks, And the score is high.

The two files are bigwig files converted by bamcoverage. Is my code wrong?

[goodez]

I found the result of macs2 is not perfect and I want to know why.

I don't understand what the problem is. Your workflow sounds okay, but what exactly is your question? Also you can use -g mm instead of using -g 2652783500. MACS2 has a predefined size for common genomes like mm10 or hg38.

I think the process is OK, but the result is weird.

I don't think anyone can help much when we don't know what you mean by "weird". We also don't know what kind of ChIP this is, so I have no idea what to expect.

[ch8316f5eyu]

I mean the good peak is that the IP sample has enrichment and input hasn't enrichment. In the two examples I show, the input sample also has enrichment on the peak called by macs2. I consider it strange. I am sorry for ambiguity.

[ch8316f5eyu]

Can you see the image I upload? The images show the peaks where IP and input both are enriched. This is my first time to upload image.

[goodez]

Yes I do see the images although they are very small. I agree it does look like there is not much enrichment of IP over the Input. However I noticed the scales for your IGV tracks are different. You should highlight both tracks, then right click and select the group autoscale. This way you can actually compare them on the same scale.

You did use -q 0.1, which allows more insignificant peaks than the default value setting of -q 0.05. Even then, I don't think you should get peaks for the regions which you showed. Can you also switch to -g mm because your genome size seems too big.

[ch8316f5eyu]

I order the table by negative log10qvalue, and these examples are the second one and the third one. The first one is also strange. I would not worry about some exception. I do change the parameter as you said and the result remains the same with less peak number. I think the calling peak is more about the peak shape than the reads number. This is the fourth one:

So I really think there is something wrong with my code and I do not know.

Answer 1

[goodez]

I would not expect MACS2 to call peaks for the regions you show, based on the tracks in your pictures. This makes me think you might have some errors in the track files themselves.

The two files are bigwig files converted by bamcoverage.

So maybe the peaks are correct, which were found using the original bam files. If you post the code used to generate the bigwigs file, the problem may become clear.

[ch8316f5eyu]

It confused me either. This time I show the bam file directly:
mass of alcl3

mass of alcl3

Here is my bowtie2 code:

bowtie2 -5 1 --very-sensitive -p 20 -I 0 -X 1000 --no-mixed --no-discordant --un-conc-gz Input_R%_un.fq.gz --met-file Input.met -x mm10 -1 Input_R1.fq.gz -2 Input_R2.fq.gz | samtools view -b | samtools sort > Input.bam&&

samtools index Input.bam

[goodez]

So what I meant was your bam files are probably okay, which is why macs2 calls peaks. BUT when you made the bigwig files, which you are viewing on IGV, you may have mixed something up. I think that both of your bigwig files may be using the input bam file in some way.