How to keep the header when I'm extracting a specific region from a bam file

Question

[vasilislenis]

Hello,

I have used bowtie2 to map some newly sequenced reads to the sheep reference genome.

After that with samtools I converted the sam file to bam file. I'm interested for the reads that are aligned in the mitochondria, so I extracted the particular region from the bam fie. So far everything seems ok, except the fact that the new bam file has "invalid BAM binary header". How can I keep the header to the new bam file?

Thank you very much in advance

Vasilis

[GouthamAtla]

Can you show the commands used? How did you extract the region of interest?

[vasilislenis]

I used:

samtools view -h B2_ACAGTG_L001_PE_bwt.sorted.bam chrM > B2_ACAGTG_L001_PE_bwt.sorted.chrM.sam

and after that I converted again into bam.

samtools view -bS B2_ACAGTG_L001_PE_bwt.sorted.chrM.sam > B2_ACAGTG_L001_PE_bwt.sorted.chrM.bam

[GouthamAtla]

Just do

samtools view -b B2_ACAGTG_L001_PE_bwt.sorted.bam chrM > B2_ACAGTG_L001_PE_bwt.sorted.chrM.bam

[vasilislenis]

Yes, it works!

Thank you very much

Answer 1

[Pierre Lindenbaum]

Try to generate a BAM file, not a SAM file. for SAM, by default samtools-view ignores the header.

samtools view  -b -o new.bam old.bam MT

[vasilislenis]

Do you mean at the mapping stage to produce directly a bam file, and after that whatever I want to extract it will keep the header?

[Jorge Amigo]

Using option -b for binary output always implies writing the header at the output, so Pierre's answer is all you need.

[vasilislenis]

Thank you very much for the explanation, it works!