Q: using macs2 to call peaks on eCLIP-seq data get an error

[qicaibiology]

Hi:

I am trying to analyze eCLIP-seq data and my PI advised me to call peaks using MACS2. I read the manual on how macs2 can be used for peak calling when doing chip-seq. It looks easy to follow: I just need to provide the macs command bam files. So, I did the following:

$ macs2 callpeak -t GRCh37SRR5111379Aligned.sortedByCoord.out.bam -c GRCh37SRR5111939Aligned.sortedByCoord.out.bam -f BAM -g 1.3e+8 -n HNRNPKRep1 --outdir macs2/ >log1 2>error1

But I got error like this:

INFO  @ Sat, 21 Nov 2020 14:51:08: 
# Command line: callpeak -t GRCh37SRR5111379Aligned.sortedByCoord.out.bam -c GRCh37SRR5111939Aligned.sortedByCoord.out.bam -f BAM -g 1.3e+8 -n HNRNPKRep1 --outdir macs2/
# ARGUMENTS LIST:
# name = HNRNPKRep1
# format = BAM
# ChIP-seq file = ['GRCh37SRR5111379Aligned.sortedByCoord.out.bam']
# control file = ['GRCh37SRR5111939Aligned.sortedByCoord.out.bam']
# effective genome size = 1.30e+08
# band width = 300
# model fold = [5, 50]
# qvalue cutoff = 5.00e-02
# The maximum gap between significant sites is assigned as the read length/tag size.
# The minimum length of peaks is assigned as the predicted fragment length "d".
# Larger dataset will be scaled towards smaller dataset.
# Range for calculating regional lambda is: 1000 bps and 10000 bps
# Broad region calling is off
# Paired-End mode is off
 
INFO  @ Sat, 21 Nov 2020 14:51:08: #1 read tag files... 
INFO  @ Sat, 21 Nov 2020 14:51:08: #1 read treatment tags... 
INFO  @ Sat, 21 Nov 2020 14:51:11:  1000000 
INFO  @ Sat, 21 Nov 2020 14:51:14:  2000000 
INFO  @ Sat, 21 Nov 2020 14:51:17:  3000000 
INFO  @ Sat, 21 Nov 2020 14:51:20:  4000000 
INFO  @ Sat, 21 Nov 2020 14:51:23:  5000000 
INFO  @ Sat, 21 Nov 2020 14:51:24: #1.2 read input tags... 
INFO  @ Sat, 21 Nov 2020 14:51:27:  1000000 
INFO  @ Sat, 21 Nov 2020 14:51:29:  2000000 
INFO  @ Sat, 21 Nov 2020 14:51:32:  3000000 
INFO  @ Sat, 21 Nov 2020 14:51:34:  4000000 
INFO  @ Sat, 21 Nov 2020 14:51:37:  5000000 
INFO  @ Sat, 21 Nov 2020 14:51:40:  6000000 
INFO  @ Sat, 21 Nov 2020 14:51:43:  7000000 
INFO  @ Sat, 21 Nov 2020 14:51:46:  8000000 
INFO  @ Sat, 21 Nov 2020 14:51:49:  9000000 
INFO  @ Sat, 21 Nov 2020 14:51:52:  10000000 
INFO  @ Sat, 21 Nov 2020 14:51:55:  11000000 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1 tag size is determined as 42 bps 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1 tag size = 42.0 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1  total tags in treatment: 1942896 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1 user defined the maximum tags... 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1 filter out redundant tags at the same location and the same strand by allowing at most 1 tag(s) 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1  tags after filtering in treatment: 1617769 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1  Redundant rate of treatment: 0.17 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1  total tags in control: 3284905 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1 user defined the maximum tags... 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1 filter out redundant tags at the same location and the same strand by allowing at most 1 tag(s) 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1  tags after filtering in control: 1846916 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1  Redundant rate of control: 0.44 
INFO  @ Sat, 21 Nov 2020 14:51:58: #1 finished! 
INFO  @ Sat, 21 Nov 2020 14:51:58: #2 Build Peak Model... 
INFO  @ Sat, 21 Nov 2020 14:51:58: #2 looking for paired plus/minus strand peaks... 
INFO  @ Sat, 21 Nov 2020 14:51:58: #2 number of paired peaks: 2 
WARNING @ Sat, 21 Nov 2020 14:51:58: Too few paired peaks (2) so I can not build the model! Broader your MFOLD range parameter may erase this error. If it still can't build the model, we suggest to use --nomodel and --extsize 147 or other fixed number instead. 
WARNING @ Sat, 21 Nov 2020 14:51:58: Process for pairing-model is terminated!

Does any one can help me correct the command to make it work?

Thanks,

[taoliu]

Hi @qicaibiology you can skip the model building by --nomodel --extsize 200. I don't think the tag shifting model for ChIP-seq can be assumed true for RNA data.

[qicaibiology]

Hi @qicaibiology you can skip the model building by --nomodel --extsize 200. I don't think the tag shifting model for ChIP-seq can be assumed true for RNA data.

Thanks, I got the bed file using your advice.