How can I locate a specific nucleotide based on a bam file?

Question

[2689563392]

Hello, now I need to get the nucleotide (from the reads, but not the reference genome) at a specific position, for example, chr4:33567, from a corresponding .bam file, and I hope the output will be a specific base, such as "A" or "T". Somebody said that samtools can do that, but I don't know how to.
So, I wanna know if it's possible and the command line to do that. Thank you!

[cpad0112]

Try IGVtools count function with window size 1 and bases option. It would give counts per each base at that position.

Answer 1

[tothepoint]

This may work for you.

samtools view input.bam "chr4:33567" > output.bam

[ATpoint]

I am pretty suer this will resut in a SAM file without header so for bam add -b option or for sam add -h, but the idea is correct I think.

[2689563392]

Thank you, but the output is not what I want. Maybe I didn't clearly articulate my purpose, and the output I want is a specific base, such as "A" or "T", but not the reads that overlap the given position.

[ATpoint]

You can pipe the resulting bam from above command into: A: BAM/SAM to FASTA conversion

[ATpoint]

Or alternatively run the relevant intervals through a variant caller and then filter the VCF file for REF and ALT alleles at that specific spot of the genome. That has the advantage that you get information on how reliable observed changes are. Not all positions in a read are of equal quality, plus variant callers often have heuristics to determine variant quality based on a variety of metrics.