Mpileup failing

[kwyattm]

I have a tumor bam file that I'm wanting to use to identify variants. I have run samtools depth and the mean depth is around 7 (kinda low), but when I start the variant identification process with samtools mpileup <filename>, the program runs for about a minute with no output, then just stops. I see no errors.

I've looked through SeqAnswers and I can't seem to find anything about this sort of error, but without an error message, it's hard to be sure it hasn't been covered yet.

There are gaps in the coverage for this bam file, but there are also gaps in other tumor samples from this experiment that work fine.

ANY help that can be offered will be greatly appreciated!!

Thanks in advance, community!

Wyatt

[deyler]

Could you copy paste the command from your terminal?

I will recuse myself if you're running the latest and greatest version of samtools, I haven't upgraded yet.

[kwyattm]

Thank you for your prompt response, deyler!

I'm using 0.1.18 (it's what I've used for everything in this project - without issue in the other samples). The most relevant command is:/cm/shared/apps/samtools mpileup -f /cm/shared/apps/genomes/hg19_genome.fa -q 1 -B $NORM $TUM > $MP,

where $NORM and $TUM are bamfiles from the Normal and Tumor, but it also fails with a simple /cm/shared/apps/samtools mpileup -f /cm/shared/apps/genomes/hg19_genome.fa -q 1 -B $TUMOR
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000

[deyler]

It sounds like that particular bam file has a problem. While I don't know what the problem is here are some commands I would run.

#get the header of the bam file
samtools view -H BAM
#get the number of reads in the bam file
samtools view -c BAM
#get the number of mapped reads in the bam file
samtools view -c -F 4 BAM

I am hoping that samtools will fail with an informative error, or at least show you something useful like there being no mapped reads in the file or something.

[kwyattm]

Thanks for your help, deyler, but I tried all three of those and there were no errors at all, and I have similar numbers of mapped reads in that tumor as I do from the other two tumors from that patient - both of which worked.

[deyler]

Bummer. I am sorry I don't know the answer but I'm willing to assist with troubleshooting if you want the help.

The next thing I'd try would be to determine the point at which mpileup is stopping.

samtools mpileup -f [FASTA] -B -q 1 [BAM] > [FILE]
tail -n1 [FILE]

You should be able to see where in its walk through the genome it stops. There's two possible outcomes:
1) no output at all
2) gets part of the way through the genome and then quits

Perhaps you're saying that 1) is already happening, and I haven't understood that.

With either outcome, you have to wonder what's different about this bam file from the other ones. With outcome 1) I would look at the header; with outcome 2) I would look specifically at the region where it quits.

Have you tried doing the alignment again? Does this sample reproducibly produce a bam file that mpileup dies on?

[kwyattm]

You're an incredible help, deyler!! Thanks so much!

Yes, you're correct that 1) is already happening. I can't tell anything about the header (everything looks okay). The output of samtools view -H <BAM> is below for two of the tumors: The first output (PT3) is from the bam file that is not giving an mpileup output, while the second output (PT2) is from the bam file that is working just fine. They seem identical!

Thanks 1^6!! You're taking a lot of time to help me with this: hopefully we can help someone else who comes across the same problem.

[kmcmaho5@kvhead checkingbams]$ samtools view -H PT3
SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
@SQ SN:chr8 LN:146364022
@SQ SN:chr9 LN:141213431
@SQ SN:chr10 LN:135534747
@SQ SN:chr11 LN:135006516
@SQ SN:chr12 LN:133851895
@SQ SN:chr13 LN:115169878
@SQ SN:chr14 LN:107349540
@SQ SN:chr15 LN:102531392
@SQ SN:chr16 LN:90354753
@SQ SN:chr17 LN:81195210
@SQ SN:chr18 LN:78077248
@SQ SN:chr19 LN:59128983
@SQ SN:chr20 LN:63025520
@SQ SN:chr21 LN:48129895
@SQ SN:chr22 LN:51304566
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566
[kmcmaho5@kvhead checkingbams]$ samtools view -H PT2
@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
@SQ SN:chr8 LN:146364022
@SQ SN:chr9 LN:141213431
@SQ SN:chr10 LN:135534747
@SQ SN:chr11 LN:135006516
@SQ SN:chr12 LN:133851895
@SQ SN:chr13 LN:115169878
@SQ SN:chr14 LN:107349540
@SQ SN:chr15 LN:102531392
@SQ SN:chr16 LN:90354753
@SQ SN:chr17 LN:81195210
@SQ SN:chr18 LN:78077248
@SQ SN:chr19 LN:59128983
@SQ SN:chr20 LN:63025520
@SQ SN:chr21 LN:48129895
@SQ SN:chr22 LN:51304566
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566

[deyler]

I would realign that sample and see if the new bamfile also does it. The only other thing I can think of to try is to run mpileup with a specified region in the middle of the genome, e.g. -r chr12:50000-60000. Maybe that would let mpileup bypass whatever problem is in the bam file.

[kwyattm]

You've been great, deyler! Thanks for your help! I'll let you know what I ultimately end up doing!