def gapsplit(args):
    """
    %prog gapsplit gffile > split.gff

    Read in the gff (normally generated by GMAP) and print it out after splitting
    each feature into one parent and multiple child features based on alignment
    information encoded in CIGAR string.
    """
    p = OptionParser(gapsplit.__doc__)

    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(not p.print_help())

    gffile, = args

    gff = Gff(gffile)
    for g in gff:
        if re.match("EST_match", g.type):
            """
            hacky implementation:
            since the standard urlparse.parse_qsl() replaces all "+" symbols with spaces
            we will write a regex to check either for a "-" or a " " (space)
            """
            match = re.search(r'\S+ (\d+) \d+ ([\s{1}\-])', g.attributes["Target"][0])
            if match.group(2) == "-":
                strand = match.group(2)
            else:
                strand = "+"
                g.attributes["Target"][0] = " ".join(str(x) \
                        for x in [g.attributes["Target"][0].rstrip(), strand])

            if g.strand == "?":
                g.strand = strand
        else:
            match = re.match(r'\S+ (\d+) \d+', g.attributes["Target"][0])
        target_start = int(match.group(1))

        re_cigar = re.compile(r'(\D+)(\d+)');
        cigar = g.attributes["Gap"][0].split(" ")
        g.attributes["Gap"] = None

        parts = []
        if g.strand == "+":
            for event in cigar:
                match = re_cigar.match(event)
                op, count = match.group(1), int(match.group(2))
                if op in "IHS":
                    target_start += count
                elif op in "DN":
                    g.start += count
                elif op == "P":
                    continue
                else:
                    parts.append([g.start, g.start + count - 1, \
                            target_start, target_start + count - 1])
                    g.start += count
                    target_start += count
        else:
            for event in cigar:
                match = re_cigar.match(event)
                op, count = match.group(1), int(match.group(2))
                if op in "IHS":
                    target_start += count
                elif op in "DN":
                    g.end -= count
                elif op == "P":
                    continue
                else:
                    parts.append([g.end - count + 1, g.end, \
                            target_start, target_start + count - 1])
                    g.end -= count
                    target_start += count

        g.update_attributes()
        print(g)

        parent = g.attributes["Name"][0]
        g.type = "match_part"
        g.attributes.clear()

        for part in parts:
            g.start, g.end = part[0], part[1]
            g.score, g.strand, g.phase = ".", g.strand, "."

            if re.match("EST", g.type):
                target_list = [parent, part[2], part[3], g.strand]
            else:
                target_list = [parent, part[2], part[3]]
            target = " ".join(str(x) for x in target_list)

            g.attributes["Parent"] = [parent]
            g.attributes["Target"] = [target]

            g.update_attributes()
            print(g)
def gapsplit ( args ) : p = OptionParser ( gapsplit . __doc__ ) opts , args = p . parse_args ( args ) if len ( args ) != 1 : sys . exit ( not p . print_help ( ) ) gffile , = args gff = Gff ( gffile ) for g in gff : if re . match ( "EST_match" , g . type ) : match = re . search ( r'\S+ (\d+) \d+ ([\s{1}\-])' , g . attributes [ "Target" ] [ 0 ] ) if match . group ( 2 ) == "-" : strand = match . group ( 2 ) else : strand = "+" g . attributes [ "Target" ] [ 0 ] = " " . join ( str ( x ) for x in [ g . attributes [ "Target" ] [ 0 ] . rstrip ( ) , strand ] ) if g . strand == "?" : g . strand = strand else : match = re . match ( r'\S+ (\d+) \d+' , g . attributes [ "Target" ] [ 0 ] ) target_start = int ( match . group ( 1 ) ) re_cigar = re . compile ( r'(\D+)(\d+)' ) ; cigar = g . attributes [ "Gap" ] [ 0 ] . split ( " " ) g . attributes [ "Gap" ] = None parts = [ ] if g . strand == "+" : for event in cigar : match = re_cigar . match ( event ) op , count = match . group ( 1 ) , int ( match . group ( 2 ) ) if op in "IHS" : target_start += count elif op in "DN" : g . start += count elif op == "P" : continue else : parts . append ( [ g . start , g . start + count - 1 , target_start , target_start + count - 1 ] ) g . start += count target_start += count else : for event in cigar : match = re_cigar . match ( event ) op , count = match . group ( 1 ) , int ( match . group ( 2 ) ) if op in "IHS" : target_start += count elif op in "DN" : g . end -= count elif op == "P" : continue else : parts . append ( [ g . end - count + 1 , g . end , target_start , target_start + count - 1 ] ) g . end -= count target_start += count g . update_attributes ( ) print ( g ) parent = g . attributes [ "Name" ] [ 0 ] g . type = "match_part" g . attributes . clear ( ) for part in parts : g . start , g . end = part [ 0 ] , part [ 1 ] g . score , g . strand , g . phase = "." , g . strand , "." if re . match ( "EST" , g . type ) : target_list = [ parent , part [ 2 ] , part [ 3 ] , g . strand ] else : target_list = [ parent , part [ 2 ] , part [ 3 ] ] target = " " . join ( str ( x ) for x in target_list ) g . attributes [ "Parent" ] = [ parent ] g . attributes [ "Target" ] = [ target ] g . update_attributes ( ) print ( g )
%prog gapsplit gffile > split.gff

    Read in the gff (normally generated by GMAP) and print it out after splitting
    each feature into one parent and multiple child features based on alignment
    information encoded in CIGAR string.