Question: How to find a variant in a particular position in a bam file only using samtools
0
gravatar for tw617
4 weeks ago by
tw61730
tw61730 wrote:

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?

sam samtools bam • 139 views
ADD COMMENTlink modified 4 weeks ago by swbarnes25.3k • written 4 weeks ago by tw61730
2
gravatar for swbarnes2
4 weeks ago by
swbarnes25.3k
United States
swbarnes25.3k wrote:

Look up how to use samtools mpileup

ADD COMMENTlink written 4 weeks ago by swbarnes25.3k

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 modified 4 weeks ago • written 4 weeks ago by tw61730
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 modified 4 weeks ago • written 4 weeks ago by h.mon24k
0
gravatar for h.mon
4 weeks ago by
h.mon24k
Brazil
h.mon24k wrote:

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 written 4 weeks ago by h.mon24k

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 modified 4 weeks ago • written 4 weeks ago by tw61730
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1521 users visited in the last hour