Pybedtools - Intersect Only Certain Feature Types
1
1
Entering edit mode
11.4 years ago
stolarek.ir ▴ 700

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)
python vcf • 2.9k views
ADD COMMENT
2
Entering edit mode
11.4 years ago
Chris Whelan ▴ 570

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.

ADD COMMENT

Login before adding your answer.

Traffic: 2333 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6