Hi, geek_y. Thank you so much. I rerun MACS2 without the --call-summits parameter and no subpeaks have the same peak boundaries. Do you mean it is because (1) I run the --call-summits parameter, (2) ATAC-Seq have many big peaks? However, I don't really understand the manual. Could you explain it more clearly for me? Thank you very very much.
Duplicated ATAC-Seq peaks called by MACS2
Hi,
I use MACS2 to do peak calling for ATAC-Seq data using the command below. However I found that some duplicated peaks called by MASC2, e.g. ChineseBulbulHead_peak_23a, ChineseBulbulHead_peak_23b, & ChineseBulbulHead_peak_23c (the detail below). I don’t meet the same problem before. Could you help me how to deal with this issue? Many thanks.
Best,
Gary
The MACS2 command I sued
macs2 callpeak -t ChineseBulbulHead.bam --format BAMPE --gsize 1000000000 --keep-dup 1 --outdir ChineseBulbulHead --name ChineseBulbulHead --verbose 3 --qvalue 0.00001 --cutoff-analysis --call-summits
The part result of ChineseBulbulHead_peaks.narrowPeak
chr1 876210 876801 ChineseBulbulHead_peak_18a 137 . 7.66032 16.78457 13.76007 143
chr1 876210 876801 ChineseBulbulHead_peak_18b 191 . 9.25623 22.30006 19.12043 517
chr1 934111 934340 ChineseBulbulHead_peak_19 117 . 7.02196 14.69713 11.74247 82
chr1 987537 987972 ChineseBulbulHead_peak_20a 400 . 12.07194 43.64410 40.03254 107
chr1 987537 987972 ChineseBulbulHead_peak_20b 295 . 9.97071 33.01239 29.59148 329
chr1 1118363 1119105 ChineseBulbulHead_peak_21 863 . 14.59966 90.58256 86.34797 330
chr1 1183894 1184262 ChineseBulbulHead_peak_22 175 . 6.85205 20.72878 17.59063 115
chr1 1509040 1510513 ChineseBulbulHead_peak_23a 204 . 5.30448 23.65111 20.43788 98
chr1 1509040 1510513 ChineseBulbulHead_peak_23b 529 . 8.87295 56.72340 52.91084 546
chr1 1509040 1510513 ChineseBulbulHead_peak_23c 291 . 6.36538 32.60884 29.19571 733
chr1 1562249 1562717 ChineseBulbulHead_peak_24 340 . 11.72668 37.55609 34.05026 149
chr1 1918313 1918729 ChineseBulbulHead_peak_25 351 . 11.95111 38.69747 35.17159 279
chr1 2086004 2086223 ChineseBulbulHead_peak_26 186 . 7.50532 21.82297 18.65607 99
chr1 2807646 2808024 ChineseBulbulHead_peak_27a 53 . 4.78770 8.04636 5.38921 65
chr1 2807646 2808024 ChineseBulbulHead_peak_27b 169 . 8.61786 20.04613 16.92561 287
• 4,094 views
•
link
updated
by
GouthamAtla
•
written
by
Gary
•
written
by
Gary
1 answer
Its not a duplicated peak. Its because of --call-summits and the nature of ATAC data.
From the manual:
--call-summits
MACS will now reanalyze the shape of signal profile (p or q-score depending on cutoff setting) to deconvolve
subpeaks within each peakcalled from general procedure. It's highly recommended to detect adjacent binding events. While used, the output subpeaks of a big peak region will have the same peak boundaries, and different scores and peak summit positions.
• 1 views
•
link
written
by
GouthamAtla
Hi, I would suggest to read a bit more about ATAC-Seq data, why it has sub-peaks, why ATAC has properties of both TFs and histone data, nucleosome positioning etc.
• 1 views
•
link
written
by
GouthamAtla
Hi geek_y, Many thanks for your suggestion. It is very helpful. Using the ChineseBulbulHead_peak_18a & 18b as an example (the figure below), I find that there are two sub-peaks. The first track is a bed file using MACS2 command below without the --call-summit parameter.
macs2 callpeak -t ChineseBulbulHead.bam --format BAMPE --gsize 1000000000 --keep-dup 1 --outdir MACS2_ChineseBulbulHead --name ChineseBulbulHead --verbose 3 --qvalue 0.00001 --cutoff-analysis
The second track is a narrowPeak file using MACS2 command below with the --call-summit parameter.
macs2 callpeak -t ChineseBulbulHead.bam --format BAMPE --gsize 1000000000 --keep-dup 1 --outdir SummitsChineseBulbulHead --name ChineseBulbulHead --verbose 3 --qvalue 0.00001 --cutoff-analysis --call-summits
The third track is a bigwig file using deepTools2 command below with the –extendReads=20. It means that I artificially add 200bp for each read.
bamCoverage --bam ChineseBulbulHead.bam --outFileName ChineseBulbulHead.bigwig --outFileFormat=bigwig --scaleFactor=1 --binSize=1 --extendReads=200 --numberOfProcessors=6 --normalizeUsingRPKM –verbose
The fourth and fifth tracks are bam files. The coverage of the bam file is more like two peaks. The bigwig file is more like a single peak.
May I say that it is because my low number of reads (~ 26 million 87bp paired-end reads) causes two sub-peaks. If I sequence more reads, the two sub-peaks will merge into single one peak, such as my bigwig file? Many thanks.
If you want to get one peak, don't use --call-summits.
• 1 views
•
link
written
by
GouthamAtla
Log in to answer this question.