Extracting a subset of bam files

[toreoe]

Hi everyone,

What I have:
A bam file with 10kb MP reads mapped back to an in house reference genome assembled in Allpaths-LG. Going to be used for visual inspection in IGV and Geneious to evaluate a possible scaffold break up suggested by linkage map.

What I want:
To extract a new bam file that contains the mapped reads for only one of the scaffolds in my reference genome.

Problem:
samtools view -b mybamfile.bam 'scaffold000046' > scf000046.bam

The resulting file lists all the original scaffolds in the header, like this:

@SQ SN:scaffold_0 LN:21965366
@SQ SN:scaffold_1 LN:18670197
@SQ SN:scaffold_2 LN:16657575
@SQ SN:scaffold_3 LN:15264954
@SQ SN:scaffold_4 LN:14680565
...

I may have misunderstood the output of the above command, but I expected to find only the target scaffold in the new bam-file header.

Any help would be much appreciated. Feel free to direct me if the solution for my problem has already been posted elsewhere. However, have not succeeded in finding one.

Best,
Tore

[jeales]

Index scf000046.bam ('samtools index scf000046.bam') and then see what 'samtools idxstats scf000046.bam' says
Are there reads reported for all the scaffolds or just scaffold000046?

[toreoe]

I extracted another, scaffold0008, from a different bam file. I then indexed scaffold00008.bam as suggested. samtools idxstats reported reads only for the target scaffold, none for the others:

scaffold00001 21966007 0 0
scaffold00002 18670670 0 0
scaffold00003 16657927 0 0
scaffold00004 15265287 0 0
scaffold00005 14680965 0 0
scaffold00006 14166276 0 0
scaffold00007 14112451 0 0
scaffold00008 13353334 2289885 55567
scaffold00009 13059713 0 0

scaffold00008 loaded properly in IGV, something I couldn't do with my previous extracted file. So it might be that there was something wrong with my original bam file.

Anyways, it seems like this is simply the way the samtools view extract command outputs the subset bam, i.e. contains all scaffolds from original bam in header. And it works fine.

Thank you very much for your help, jeales.

[jeales]

Thats good, if there were reads from other scaffolds an awful lot of my work would be wrong as well!

I don't have a good solution for trimming the header to only include scaffolds/chromosomes that have reads present in the bam, other than manual editing

But it's not realyl doing any harm as long as the counts are zero

[dpryan]

The header is left as-is for technical reasons. In BAM files (and internally by samtools and likely picard too) the chromosome/contig field of an alignment is just a numeric index into this list in the header. Trimming the header, then would result in needing to alter every read as it was output, which would be annoying (this is also why people need to be careful when they reheader BAM files, since changing the order of the contigs/chromosomes in the header will make all of the alignments wrong), though I guess not that hard to implement (it would affect performance, though).

[jeales]

Very informative thanks
So if all the alignments in a bam are mapped to the 10th chromosome/contig in the header their chromosome/contig will be 10 (or probably 9, if you start counting at 0), irrespective of chromosome name?

[dpryan]

Exactly. As you surmised, everything is 0-based internally, so the value is 9. For those curious, the value is unsigned, so a value of -1 used for unmapped reads.