Pybedtools - Intersect Only Certain Feature Types

Question

[stolarek.ir]

Hi. I am doing the intersections on vcf files using pybedtools. However I'm interested in interescting only features of a given type (insertion in A with the position in B if this position is also an insertion). Also it would be nice to report while intersecting if the feature type did match or not between files.

I have written some code, but it's clumsy I think, as it retrieves the correct SVs types to separate lists, saves them as files, and then performs intersections.

import pybedtools   
    d=[]
    for f in a:
        if 'INS' in f[7]:
            d.append(f)
    g=[]
    for f in b:
        if 'INS' in f[7]:
            g.append(f)
    for i in d:
        aout.write(str(i)+'\n')
    for i in g:
        bout.write(str(i)+'\n')
    c = pybedtools.BedTool('/home/istolarek/aout')
    e = pybedtools.BedTool('/home/istolarek/bout')

    ab = c.intersect(e, f=0.9,r=False)

Answer 1

[Chris Whelan]

You can use pybedtools filters to get the right intervals out (untested code):

def insertion_filter(feature):
    return feature[7] == 'INS'

ab = a.filter(insertion_filter).intersect(b.filter(insertion_filter), f=0.9, r=False)

The second part of your question (reporting intersections that don't match) is slightly trickier, but you should be able to do it by intesecting a and b with the "write both" option and then running a filter similar to the one above that compares the feature fields to make sure they match.

Pybedtools is great.