Get the positive and negative strands of reads in bam file for particular nucleotide position?

Question

[murali]

Get the positive and negative strands of reads in bam file for particular nucleotide position?

For example, If I have the bam file contains Chromosome:3 nucleotide_position:13226978, I want to find the no of reads mapped to particular nucleotide_position. It can be done by the pysamstats package using reads_all function.

Now, I want to find out positive and negative strands or reads to that particular nucleotide_position as well as gene-symbol associated with it?

How can I do this and can able extract relevant information?

I want to do it for list of bam files, so I am writing python script.

Any suggestions related to extracting strands as well as gene-symbol are really helpful to me. :)

Answer 1

[Devon Ryan]

For getting a count of the orientation of alignments overlapping a position, you would use pysam with the pileup() function. You would then iterate over the PileupReads returned for a given position and increment a counter if an alignment is reverse complemented.

For the annotation aspect, pybedtools would be the logical choice if you're already writing stuff in python.