How to get a list of reads supporting a variant in R?

Question

[abascalfederico]

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!

[Sean Davis]

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

[abascalfederico]

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

Answer 1

[Devon Ryan]

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

[abascalfederico]

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?

[Devon Ryan]

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.

[abascalfederico]

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

[Devon Ryan]

Sometimes the hacky solution is the best solution :P