Question: How to find a variant in a particular position in a bam file only using samtools
0
gravatar for tw617
15 months ago by
tw61740
tw61740 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 • 1.0k views
ADD COMMENTlink modified 15 months ago by swbarnes27.9k • written 15 months ago by tw61740
2
gravatar for swbarnes2
15 months ago by
swbarnes27.9k
United States
swbarnes27.9k wrote:

Look up how to use samtools mpileup

ADD COMMENTlink written 15 months ago by swbarnes27.9k

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 15 months ago • written 15 months ago by tw61740
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 15 months ago • written 15 months ago by h.mon30k
0
gravatar for h.mon
15 months ago by
h.mon30k
Brazil
h.mon30k 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 15 months ago by h.mon30k

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 15 months ago • written 15 months ago by tw61740
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: 963 users visited in the last hour