Does samtools view -b -h <input.bam> -L <input.bed> keep the original header?

Question

[rubic]

Simple question

Seems like the command:

samtools view -b -h <input.bam> -L <input.bed>

for subsetting the bam file to keep only those reads intersecting the seqnames and their regions in input.bed is retaining the original header in the output bam file, including those seqnames not present in input.bed.

Is there a samtools flag that can be used to tell samtools to also subset the header to keep only those seqnames in input.bed?

Answer 1

[d-cameron]

Is there a samtools flag that can be used to tell samtools to also subset the header to keep only those seqnames in input.bed?

No. You would need to pipe to samtools reheader to do that.

The BAM file format differs from SAM in that RNAME is stored not as a string but as the zero-based integer offset of the match @SQ header instead of the seqname string. This means that a) BAM files must have @SQ headers, and b) any subsetting of @SQ requires rewriting all BAM records with the new offsets. samtools view does not have this capability.