This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to get a list of reads supporting a variant in R?

For a given position, in which there is a variant called by some algorithm, I want to get the pileup and analyse the reads supporting that variant. Any idea on how to do this in R? I know how to do this in Python using pysam but I need an R implementation.

Thanks!

r bam
• 3,379 views • link
updated by Devon Ryan Devon Ryan • written by abascalfederico abascalfederico

What is the pysam functionality you want to reproduce, specifically? A few lines of code would help.

• 0 views • link
written by Sean Davis Sean Davis

For a given position, I want to read the pileup and when there is a match to the variant, I want to get the mapping coordinates of the read matching the variant. I do not have the exact code yet, but it would imply something like:

for pileupcolumn in samfile.pileup(chr,int(pos)-1,int(pos)):
    base = pileupcolumn.pos
    if(base == (int(pos)-1)) :
        print('SITE_OVERLAP:%s:%s:%s' % (str(chr),pos,position) )
        coverage = pileupcolumn.n
        ref_position = pileupcolumn.reference_pos
        ref_base = reffile.fetch(chr,ref_position,ref_position+1)
        for pileupread in pileupcolumn.pileups:
                qry_len = pileupread.alignment.query_length
                read_base = pileupread.alignment.query_sequence[pileupread.query_position]
                if read_base==variant_base: 
                    dist_to_3p = qry_len - pileupread.query_position - 1
                    dist_to_5p = pileupread.query_position + 1
                    ...

By doing that I can access information about each read and relate that to the presence/absence of a variant match

• 0 views • link
written by abascalfederico abascalfederico

1 answer

You're looking for the Rsamtools package, particularly the pileup() command. This and Rhtslib combine to form the R equivalent to pysam in python.

• 0 views • link
written by Devon Ryan Devon Ryan

Thanks Devon, I already had a look at that function, but it just provide a table of counts. I need to identify which reads support a given variant. This is relatively easy with pysam but I find no way to do this in Rsamtools. Maybe I missed something from the documentation?

• 0 views • link
written by abascalfederico abascalfederico

It looks like Rsamtools isn't exposing this, which is unfortunate. You can get around this by fetching all reads overlapping a position, but of course that puts a fair bit of the programmatic burden on you.

• 0 views • link
written by Devon Ryan Devon Ryan

Thanks, Devon. I don't want the programmatic burden ;-) I will instead make a call to python from R - not what I wanted though

• 0 views • link
written by abascalfederico abascalfederico

Sometimes the hacky solution is the best solution :P

• 0 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this