@@ -0,0 +1,2461 @@
+#!/usr/bin/env python%0A%0A%0Aimport os%0Aimport sys%0Aimport cPickle as pickle%0Aimport logging%0A%0Alogger = logging.getLogger(__name__)%0A%0Afrom collections import OrderedDict as od%0A%0Adef resolve_duplicates(insertions):%0A    ''' resolve instances where breakpoints occur %3E 1x in the insertion list '''%0A    ''' this can happen if intervals overlap, e.g. in  genome chunking '''%0A%0A    insdict = od() # --%3E index in insertions%0A%0A    for n, ins in enumerate(insertions):%0A        be1 = ins%5B'INFO'%5D%5B'chrom'%5D + ':' + str(ins%5B'INFO'%5D%5B'be1_breakpos'%5D)%0A        be2 = ins%5B'INFO'%5D%5B'chrom'%5D + ':' + str(ins%5B'INFO'%5D%5B'be2_breakpos'%5D)%0A        %0A        if be1 not in insdict:%0A            insdict%5Bbe1%5D = n %0A            insdict%5Bbe2%5D = n%0A        %0A        else:%0A            if prefer_insertion(ins, insertions%5Binsdict%5Bbe1%5D%5D):%0A                insdict%5Bbe1%5D = n%0A                insdict%5Bbe2%5D = n%0A%0A    return %5Binsertions%5Bn%5D for n in list(set(insdict.values()))%5D%0A%0A%0Adef prefer_insertion(ins1, ins2):%0A    ''' return true if ins1 has more evidence than ins2, false otherwise '''%0A    # prefer two-end support over one end%0A    if ins1%5B'INFO'%5D%5B'be1_breakpos'%5D != ins1%5B'INFO'%5D%5B'be2_breakpos'%5D and ins2%5B'INFO'%5D%5B'be1_breakpos'%5D == ins2%5B'INFO'%5D%5B'be2_breakpos'%5D:%0A        return True%0A%0A    # prefer higher split read count%0A    if ins1%5B'INFO'%5D%5B'be1_sr_count'%5D + ins1%5B'INFO'%5D%5B'be2_sr_count'%5D %3E ins2%5B'INFO'%5D%5B'be1_sr_count'%5D + ins2%5B'INFO'%5D%5B'be2_sr_count'%5D:%0A        return True%0A%0A    # prefer higher discordant read count%0A    if ins1%5B'INFO'%5D%5B'dr_count'%5D %3E ins2%5B'INFO'%5D%5B'dr_count'%5D:%0A        return True%0A%0A    return False%0A%0A%0Aif __name__ == '__main__':%0A    FORMAT = '%25(asctime)s %25(message)s'%0A    logging.basicConfig(format=FORMAT)%0A    logger.setLevel(logging.INFO)%0A%0A%0Aif len(sys.argv) %3E 2:%0A    assert not os.path.exists(sys.argv%5B1%5D), %22file exists: %25s, refusing to overwrite.%22 %25 sys.argv%5B1%5D%0A%0A    insertions = %5B%5D%0A%0A    for pfn in sys.argv%5B2:%5D:%0A        with open(pfn, 'r') as pickin:%0A            insertions += pickle.load(pickin)%0A%0A        logger.info('loaded %25s' %25 pfn)%0A%0A    logger.info('finished loading pickles, found %25d records, merging results...' %25 len(insertions))%0A%0A    insertions = resolve_duplicates(insertions)%0A%0A    logger.info('%25d records remain after merge.' %25 len(insertions))%0A%0A    with open(sys.argv%5B1%5D, 'w') as pickout:%0A        pickle.dump(insertions, pickout)%0A%0A    logger.info('wrote %25d records to %25s' %25 (len(insertions), sys.argv%5B1%5D))%0A%0Aelse:%0A    sys.exit(%22usage: %25s %3Coutput.pickle%3E %3Cinput1.pickle%3E %3Cinput2.pickle%3E ... %22)

add script