mpileup: specified region truncated.

[lbthrice]

Hello gentlepeople,

I am using the samtools 'mpileup' program.
I wish to generate a pileup file for a specific region.
I have several separate .bam input files that I am generating the pileups for.
Unfortunately, mpileup truncates the region and I only get data for the first part of the region.
I have ran the same command on different .bam files and get different results.

When I command:

samtools mpileup -f 'hg19.fa' -r chr6:27114408-27115845 'experiment01.bam' > out_PILEUP2.pile
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000

I get a file with 746 lines.
I am expecting 1438 lines because that is the size of the interval I have specified with the '-r' flag
If I use mpileup with the same region but with a different .bam file I get a different size file.

samtools mpileup -f 'hg19.fa' -r chr6:27114408-27115845 'experiment02.bam' > out_PILEUP3.pile
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000

I get a file with 623 lines.
Again, I am expecting 1438 lines because that is the size of the interval I have specified with the '-r' flag.

The first line in both files is for chr6:27114408, as expected.
BUT, the files are truncated at different positions.
With experiment01.bam the final line of the pileup file is chr6:27115577.
With experiment02.bam the final line of the pileup file is chr6:27115344.

This must have something to do with the .bam file.
Can someone tell me which parameter needs to be adjusted?
I have consulted the manual but I could not identify anything that sounded applicable.

Thanks for you time,
Lionel (Lee) Brooks 3rd
Dartmouth Genetics Grad Student

[swbarnes2]

The simplest answer is that you just don't have any reads across part of your region.

Try getting the .sam file for the same region +- 500 bases, or try looking at the whole .bam in IGV. You can zoom in on your region.

[gringer]

Yeah, mpileup doesn't generate lines for regions it has no reads on. If you need to work out whole-region coverage, then you should be looking at the column that gives the base number. For covered regions, this will increase by 1 per line (assuming no inserts), but will jump multiple bases when there are no reads.

[lbthrice]

whoops, should have checked that...thanks!