def dup(args):
    """
    %prog dup frgscffile

    Use the frgscf posmap file as an indication of the coverage of the library.
    Large insert libraries are frequently victims of high levels of redundancy.
    """
    p = OptionParser(dup.__doc__)
    opts, args = p.parse_args(args)

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

    frgscffile, = args

    fp = open(frgscffile)
    data = [FrgScfLine(row) for row in fp]
    # we need to separate forward and reverse reads, because the position
    # collisions are handled differently
    forward_data = [x for x in data if x.orientation == '+']
    reverse_data = [x for x in data if x.orientation == '-']

    counts = defaultdict(int)
    key = lambda x: (x.scaffoldID, x.begin)
    forward_data.sort(key=key)
    for k, data in groupby(forward_data, key=key):
        data = list(data)
        count = len(data)
        counts[count] += 1

    key = lambda x: (x.scaffoldID, x.end)
    reverse_data.sort(key=key)
    for k, data in groupby(forward_data, key=key):
        data = list(data)
        count = len(data)
        counts[count] += 1

    prefix = frgscffile.split(".")[0]
    print("Duplication level in `{0}`".format(prefix), file=sys.stderr)
    print("=" * 40, file=sys.stderr)
    for c, v in sorted(counts.items()):
        if c > 10:
            break
        label = "unique" if c == 1 else "{0} copies".format(c)
        print("{0}: {1}".format(label, v), file=sys.stderr)
def dup ( args ) : p = OptionParser ( dup . __doc__ ) opts , args = p . parse_args ( args ) if len ( args ) != 1 : sys . exit ( p . print_help ( ) ) frgscffile , = args fp = open ( frgscffile ) data = [ FrgScfLine ( row ) for row in fp ] forward_data = [ x for x in data if x . orientation == '+' ] reverse_data = [ x for x in data if x . orientation == '-' ] counts = defaultdict ( int ) key = lambda x : ( x . scaffoldID , x . begin ) forward_data . sort ( key = key ) for k , data in groupby ( forward_data , key = key ) : data = list ( data ) count = len ( data ) counts [ count ] += 1 key = lambda x : ( x . scaffoldID , x . end ) reverse_data . sort ( key = key ) for k , data in groupby ( forward_data , key = key ) : data = list ( data ) count = len ( data ) counts [ count ] += 1 prefix = frgscffile . split ( "." ) [ 0 ] print ( "Duplication level in `{0}`" . format ( prefix ) , file = sys . stderr ) print ( "=" * 40 , file = sys . stderr ) for c , v in sorted ( counts . items ( ) ) : if c > 10 : break label = "unique" if c == 1 else "{0} copies" . format ( c ) print ( "{0}: {1}" . format ( label , v ) , file = sys . stderr )
%prog dup frgscffile

    Use the frgscf posmap file as an indication of the coverage of the library.
    Large insert libraries are frequently victims of high levels of redundancy.