Output BAM with SAM flags retained

197 views
Skip to first unread message

Lucy

unread,
Aug 28, 2018, 5:37:39 AM8/28/18
to deepTools
Hi,

I am using deepTools alignmentSieve to shift the reads in my ATAC-seq BAM files prior to peakcalling with MACS2. However, when I run MACS2 with --nomodel --shift -100 -extsize options, I get the following error:
"Exception: File is not of a valid BAM format! 0"

My BAM file looks as follows after deepTools:
NB501605:118:HM5MMBGX5:2:21211:25279:5516 163 chr1 10531 30 36M = 10570 74 * *
NB501605:130:HN53GBGX5:4:22611:19989:15930 163 chr1 10532 30 36M = 10556 58 * *
NB502048:39:HT3HMBGX3:3:23403:16439:14937 147 chr1 10551 30 34M = 10518 -67 * *
NB501605:103:HJTV3BGX5:2:22303:13575:13151 83 chr1 10556 30 34M = 10388 -202 * *
NB501605:103:HJTV3BGX5:2:23308:22908:2460 83 chr1 10556 30 34M = 10176 -414 * *
NB501605:103:HJTV3BGX5:4:23410:15093:20233 83 chr1 10556 30 32M = 10449 -139 * *
NB501605:118:HM5MMBGX5:2:23212:2944:8687 147 chr1 10556 30 32M = 10091 -497 * *
NB501605:130:HN53GBGX5:4:22611:19989:15930 83 chr1 10556 30 34M = 10532 -58 * *
NB501605:103:HJTV3BGX5:1:12311:19279:13650 83 chr1 10560 30 35M = 10376 -219 * *
NB501605:103:HJTV3BGX5:4:11404:11057:4980 83 chr1 10560 30 35M = 10518 -77 * *
NB502048:39:HT3HMBGX3:2:12211:4641:7529 83 chr1 10562 30 35M = 10375 -222 * *
NB501605:103:HJTV3BGX5:1:12303:8882:13722 83 chr1 10568 34 34M = 10513 -89 * *
NB501605:118:HM5MMBGX5:4:22609:20966:10318 147 chr1 10568 30 32M = 10450 -150 * *
NB501605:103:HJTV3BGX5:4:11404:4695:1867 83 chr1 10570 30 35M = 10445 -160 * *
NB501605:118:HM5MMBGX5:2:21211:25279:5516 83 chr1 10570 30 35M = 10531 -74 * *
NB501605:130:HN53GBGX5:4:11402:16163:19491 163 chr1 10582 30 35M = 11220 673 * *
NB501605:118:HM5MMBGX5:1:13101:20188:15959 99 chr1 10585 30 35M = 11226 675 * *
NB501605:118:HM5MMBGX5:3:13402:1878:15578 99 chr1 10585 30 35M = 11277 726 * *
NB501605:130:HN53GBGX5:3:12401:1570:15711 163 chr1 10585 30 35M = 10933 383 * *
NB501605:103:HJTV3BGX5:2:13107:10067:10868 163 chr1 10587 40 36M = 11169 617 * *
NB501605:103:HJTV3BGX5:3:13407:20587:16102 99 chr1 10587 30 36M = 11277 724 * *

The SAM flags are missing and this is why I am getting an error with MACS2 I think, because it can't detect which is the 1st and 2nd read in a pair. Is it possible to output a BAM file retaining the SAM flags after deepTools, or how would you suggest to overcome this problem?

Many thanks,
Lucy

Devon Ryan

unread,
Aug 28, 2018, 5:42:29 AM8/28/18
to lucyc...@gmail.com, deepTools
The SAM flags are present in the second column. Can you make the BAM
file (privately) available and provide the full macs2 command you're
using? We process a fair number of ATAC-seq samples in exactly this
manner and haven't run into this problem, so perhaps there's an odd
edge-case somewhere.

Devon
--
Devon Ryan, Ph.D.
Email: dpr...@dpryan.com
Data Manager/Bioinformatician
Max Planck Institute of Immunobiology and Epigenetics
Stübeweg 51
79108 Freiburg
Germany
> --
> You received this message because you are subscribed to the Google Groups "deepTools" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to deeptools+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Devon Ryan

unread,
Aug 28, 2018, 7:32:19 AM8/28/18
to lucyc...@gmail.com, deepTools
This is due to a bug in MACS2 where it uses incorrect fields in the
BAM file (specifically, it incorrectly looks at the read length to
determine the fragment length). Thankfully, this only happens in the
"sniffer" part, so if you specify `-f BAMPE` or `-f BAM` then you will
avoid the bug.

Devon

--
Devon Ryan, Ph.D.
Email: dpr...@dpryan.com
Data Manager/Bioinformatician
Max Planck Institute of Immunobiology and Epigenetics
Stübeweg 51
79108 Freiburg
Germany

On Tue, Aug 28, 2018 at 11:37 AM Lucy <lucyc...@gmail.com> wrote:
>

Lucy

unread,
Aug 28, 2018, 8:14:53 AM8/28/18
to deepTools
Ok, thank you.

For -BAMPE, --shift is set to 0 according to the MACS2 manual.
"Note, you can't set values other than 0 if format is BAMPE or BEDPE for paired-end data. Default is 0."

With -f BAM, is it possible to use --shift -100 i.e. values other than 0?

Devon Ryan

unread,
Aug 28, 2018, 8:17:34 AM8/28/18
to lucyc...@gmail.com, deepTools
That'd be a question for the MACS2 author, we only use BAMPE.

--
Devon Ryan, Ph.D.
Email: dpr...@dpryan.com
Data Manager/Bioinformatician
Max Planck Institute of Immunobiology and Epigenetics
Stübeweg 51
79108 Freiburg
Germany

Reply all
Reply to author
Forward
0 new messages