Does BBMap not include headers in BAM files?

Question

[Sinji]

Been working with some ChIPseq data, and am trying to run it through ChIPQC from Bioconductor. I continually get this error:

[bam_header_read] EOF marker is absent. The input is probably truncated.
  [bam_header_read] invalid BAM binary header (this is not a BAM file).
  Error in value[[3L]](cond) : 
    failed to open BamFile: SAM/BAM header missing or empty
    file: '/usr/local/anaconda2/lib/R/bin/exec/R'
  Calls: ChIPQCsample ... tryCatch -> tryCatchList -> tryCatchOne -> <Anonymous>
  Execution halted

I then samtools view | head my bam file and sure enough I don't get any header information:

NS500579:38:HC3MLBGXX:3:22605:6021:20111    16  chr1    9980    2   49M *   0   0   GCCTACGCGTCCGCGTCCGCGTAACCCTAACCCTAACCCTAACCCTAAC   //////A//A///E/E/////6/A<E6EAEAEEEEEEEEEEEEEAAAAA   XT:A:R  NM:i:21 AM:i:2
NS500579:38:HC3MLBGXX:3:21402:17333:9302    16  chr1    9981    2   49M *   0   0   CATACGTGTGCTCTTCCGATTAACCCTAACCCTAACCCTAACCCTAACC   6//E/A//E///E/E///////EEEEEEEEEEAEEEEEEEEEEEAAAAA   XT:A:R  NM:i:20 AM:i:2
NS500579:38:HC3MLBGXX:1:11212:16057:4322    16  chr1    9981    2   49M *   0   0   GTTACTCATCAGATTACGAGTAACCCTAACCCTAACCCTAACCCTAACC   <//6//////////A/////E/////EAAE/AEEEEEEEEEEEEAA/AA   XT:A:R  NM:i:20 AM:i:2
NS500579:38:HC3MLBGXX:3:23606:5285:10927    16  chr1    10000   2   49M *   0   0   GTAACCGTTACCGTAACCGTAACCCTAACCCTAACCCTAACCCTAACCC   /////E/////E/////E/////E////EEA6EEEEEAEEEEEEAAAAA   XT:A:R  NM:i:5  AM:i:2
NS500579:38:HC3MLBGXX:3:13511:5485:4267 16  chr1    10003   2   49M *   0   0   ACGCTTCCGCTACCGCTCACCCTAACCCTAACCCTAACCCAAACCCTAA   /E/////A/AE//E//E//EEEEE/EAAEEEEEEAEEEEE/EEAAAAAA   XT:A:R  NM:i:8  AM:i:2
NS500579:38:HC3MLBGXX:1:23206:6162:17684    0   chr1    10016   3   49M *   0   0   CCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCAAACCCAAAC   AAAAAEEEEEEEEEEEEEEEEEEAEE<AE6EAE/E/E///E/A/</E//   XT:A:R  NM:i:2  AM:i:3
NS500579:38:HC3MLBGXX:1:22311:17702:3502    0   chr1    10023   3   49M *   0   0   CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC   AAAAAEEEEEEEEEEEEAEEEEE/E//AE/E//AE/<<///E////EE/   XT:A:R  NM:i:0  AM:i:3
NS500579:38:HC3MLBGXX:4:21408:9090:4962 16  chr1    10024   3   49M *   0   0   CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC   ////EEE////E/E///EEEE/E/E/E/EEEEAEAEEAEEAAEE//AAA   XT:A:R  NM:i:0  AM:i:3
NS500579:38:HC3MLBGXX:3:11512:25155:2180    0   chr1    10026   2   49M *   0   0   AACCCAAACCCTAACCCTAACCCTAACCCAACCCCTAACCCAAACCATA   AAAA//EEEAAEEAAAEAE/E/AEE//EE/////<A///E</A/EA/A/   XT:A:R  NM:i:5  AM:i:2
NS500579:38:HC3MLBGXX:4:11610:17453:12297   16  chr1    10050   3   49M *   0   0   AACCCTAACCCTAACCCAAACCCTAACCCTAACCCTAACCCTAACCCTA   //////A//////A/<//E//AE//6EAEAE//EE/EEEAAEEEAA/AA   XT:A:R  NM:i:1  AM:i:3

I figure it's something about the way i'm mapping and then converting the files, or maybe something specific with sambama, i'm not entirely sure, I can go ahead and test it out but I figure perhaps someone form Biostars can answer this quickly enough that I don't have to bother.

Here's my pipeline:

 bbmap.sh in=${read1} outm=${id}.mapped.bam outu=${id}.unmapped.bam keepnames=t trd sam=1.3 maxindel=1 ambig=random statsfile=${id}.alignmentReport.txt minid=${params.minid} usemodulo

 sambamba sort --tmpdir $baseDir -t ${params.threads} -o ${id}.sorted.mapped.bam ${id}.mapped.bam

 sambamba index -t ${params.threads} ${id}.sorted.mapped.bam

Thanks in advance!

EDIT: As Devon pointed out, I DO have headers.

[GenoMax]

I doubt BBMap is the culprit here :) What version of samtools are you using?

[Sinji]

Am not using samtools, instead using sambamba 0.6.5.

EDIT: Since BBMap does use samtools to convert to bam automatically, samtools version is 1.3.1.

Answer 1

[Sinji]

The culprit was my R script for ChIPQC doing some funny stuff with arguments. BBMap, Sambamba, and Samtools were not the problem.

Answer 2

[Devon Ryan]

Your file has a header, otherwise you wouldn't see chromosome names with samtools view, where you forgot -h.

BTW, the error message suggests that you have a truncated file, since it didn't end with the magic number.

[Sinji]

Whoops. You're right. Any idea why ChIPQC isn't recognizing the header then?

EDIT: NVM, just saw your edit. Doesn't make any sense to me, the file maps just fine, and it works if I use Bowtie2 as my mapper. Odd.

[Devon Ryan]

My only guess is that sambamba produced a broken BAM file. Sort the BAM file with samtools and see what happens.

[Sinji]

After sorting and indexing with samtools, it appears I still have the same problem.