Biostars beta testing.
Question: How to find a variant in a particular position in a bam file only using samtools
0
Entering edit mode

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?

ADD COMMENTlink 10 months ago tw617 • 30 • updated 10 months ago swbarnes2 5.7k
2
Entering edit mode

Look up how to use samtools mpileup

ADD COMMENTlink 10 months ago swbarnes2 5.7k
Entering edit mode
0

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!

ADD REPLYlink 10 months ago
tw617
• 30
Entering edit mode
1

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.

ADD REPLYlink 10 months ago
h.mon
25k
0
Entering edit mode

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.

ADD COMMENTlink 10 months ago h.mon 25k
Entering edit mode
0

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.

ADD REPLYlink 10 months ago
tw617
• 30

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.0