Discrepancy in uniquely mapped reads between samtools and macs (1.4 or 2)

[cnov...@gmail.com]

Hi Tao,

I have noticed a discrepancy in the uniquely mapped reads identified by samtools or either version of macs. I was wondering if you had any ideas about what could be happening? I have posted all necessary outputs, and all was performed using the same bam file.

# samtools view -c -F 4 /Users/cterranova/SKCM2-TCGA-SKCM-M028-SR-P010-B-NC-CJT.sorted.bam

34653863

#samtools flagstat /Users/cterranova/SKCM2-TCGA-SKCM-M028-SR-P010-B-NC-CJT.sorted.bam

60275724 + 0 in total (QC-passed reads + QC-failed reads)

19300039 + 0 duplicates

34653863 + 0 mapped (57.49%:nan%)

#macs14 -t /Users/cterranova/SKCM2-TCGA-SKCM-M028-SR-P010-B-NC-CJT.sorted.bam -f BAM -g hs --nodmodel -n TEST1

INFO  @ Mon, 19 Sep 2016 10:14:30: #1 tag size is determined as 36 bps 

INFO  @ Mon, 19 Sep 2016 10:14:30: #1 tag size = 36 

INFO  @ Mon, 19 Sep 2016 10:14:30: #1  total tags in treatment: 15353824 

INFO  @ Mon, 19 Sep 2016 10:14:30: #1 user defined the maximum tags... 

INFO  @ Mon, 19 Sep 2016 10:14:30: #1 filter out redundant tags at the same location and the same strand by allowing at most 1 tag(s) 

INFO  @ Mon, 19 Sep 2016 10:14:33: #1  tags after filtering in treatment: 15353824 

INFO  @ Mon, 19 Sep 2016 10:14:33: #1  Redundant rate of treatment: 0.00 

INFO  @ Mon, 19 Sep 2016 10:14:33: #1 finished! 

INFO  @ Mon, 19 Sep 2016 10:14:33: #2 Build Peak Model... 

INFO  @ Mon, 19 Sep 2016 10:14:33: #2 Skipped... 

INFO  @ Mon, 19 Sep 2016 10:14:33: #2 Use 100 as shiftsize, 200 as fragment length 

INFO  @ Mon, 19 Sep 2016 10:14:33: #3 Call peaks... 

INFO  @ Mon, 19 Sep 2016 10:14:33: #3 shift treatment data 

INFO  @ Mon, 19 Sep 2016 10:14:36: #3 merge +/- strand of treatment data 

INFO  @ Mon, 19 Sep 2016 10:14:36: #3 call peak candidates 

INFO  @ Mon, 19 Sep 2016 10:16:39: #3 use self to calculate local lambda and  filter peak candidates... 

INFO  @ Mon, 19 Sep 2016 10:17:10: #3 Finally, 14077 peaks are called! 

INFO  @ Mon, 19 Sep 2016 10:17:10: #4 Write output xls file... TEST1_peaks.xls 

INFO  @ Mon, 19 Sep 2016 10:17:11: #4 Write peak bed file... TEST1_peaks.bed 

INFO  @ Mon, 19 Sep 2016 10:17:11: #4 Write summits bed file... TEST1_summits.bed 

INFO  @ Mon, 19 Sep 2016 10:17:11: #5 Done! Check the output files!

#macs2 callpeak -t /Users/cterranova/SKCM2-TCGA-SKCM-M028-SR-P010-B-NC-CJT.sorted.bam -f BAM -g hs --nomodel -n TEST2

INFO  @ Mon, 19 Sep 2016 10:23:02: #1 tag size is determined as 36 bps 

INFO  @ Mon, 19 Sep 2016 10:23:02: #1 tag size = 36 

INFO  @ Mon, 19 Sep 2016 10:23:02: #1  total tags in treatment: 15353824 

INFO  @ Mon, 19 Sep 2016 10:23:02: #1 user defined the maximum tags... 

INFO  @ Mon, 19 Sep 2016 10:23:02: #1 filter out redundant tags at the same location and the same strand by allowing at most 1 tag(s) 

INFO  @ Mon, 19 Sep 2016 10:23:02: #1  tags after filtering in treatment: 15353824 

INFO  @ Mon, 19 Sep 2016 10:23:02: #1  Redundant rate of treatment: 0.00 

INFO  @ Mon, 19 Sep 2016 10:23:02: #1 finished! 

INFO  @ Mon, 19 Sep 2016 10:23:02: #2 Build Peak Model... 

INFO  @ Mon, 19 Sep 2016 10:23:02: #2 Skipped... 

INFO  @ Mon, 19 Sep 2016 10:23:02: #2 Use 200 as fragment length 

INFO  @ Mon, 19 Sep 2016 10:23:02: #3 Call peaks... 

INFO  @ Mon, 19 Sep 2016 10:23:02: #3 Pre-compute pvalue-qvalue table... 

INFO  @ Mon, 19 Sep 2016 10:23:32: #3 Call peaks for each chromosome... 

INFO  @ Mon, 19 Sep 2016 10:23:51: #4 Write output xls file... TEST2_peaks.xls 

INFO  @ Mon, 19 Sep 2016 10:23:51: #4 Write peak in narrowPeak format file... TEST2_peaks.narrowPeak 

INFO  @ Mon, 19 Sep 2016 10:23:51: #4 Write summits bed file... TEST2_summits.bed 

INFO  @ Mon, 19 Sep 2016 10:23:51: Done!

Thanks very much for the help!!