How to find a variant in a particular position in a bam file only using samtools

Question

[tw617]

I am trying to find a way to display the presence of a single mismatch at position chr7:17375390 on the reverse strand of an overlapping read in a bam file using only samtools. I sliced the section of the chromosome, and I can see that there is a C aligning to T in the reference at this position in IGV. But when I use samtools view and look at the cigar files, it only says 55M for this read (no information about the mismatch). I realize that M stands for matches and mismatches but how can I get around this. I think I might need to make a BED file and pass it to this bam file with -L and look directly at the sequences of this position but can't seem to get this to work.

$ samtools view -L bedfile.bed -h subset2.bam > out.sam

This seems to have correctly made the bed file, there is a blue line under my region of interest in IGV, the coverage says 34 and :

$     samtools view -c out.sam

34

So now I can look at just the sections of the reads that are in this region: How do I line these up and visualize single mismatch that I am seeing in IGV?

Answer 1

[swbarnes2]

Look up how to use samtools mpileup

[tw617]

THANK YOU!

 $ samtools mpileup -v -u out.bam | grep "17375390" | cut -f 8 > out.vcf

Now I just need to learn how to interpret a vcf file!

[h.mon]

In recent versions, mpileup has been moved from samtools to bcftools, and you don't need the grep, you can use -r:

 samtools mpileup -r chr7:17375390 out.sam

The resulting file is a pileup file, not a vcf file.

Answer 2

[h.mon]

Call variants for the region of interest only (-r chr:from-to), output vcf, and use the vcf to visualize the variant with IGV. Only problem is bcftools instead of samtools.

[tw617]

Thanks for the reply but I can already see it in IGV, I need to be able to show it is there using only the bam file. I think I might need to use a particular flag on my bed file modified bam but I am having trouble figuring out which one I need or if I need them at all.