This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How can I locate a specific nucleotide based on a bam file?

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!

alignment genome sequence
• 1,352 views • link
written by 2689563392 • 2689563392

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

• 0 views • link
written by cpad0112 cpad0112

1 answer

This may work for you.

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

• 0 views • link
written by tothepoint tothepoint

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.

• 0 views • link
written by ATpoint ATpoint

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.

• 0 views • link
written by 2689563392 • 2689563392

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.

• 0 views • link
written by ATpoint ATpoint

Log in to answer this question.

More posts like this