<commit_before>#!/usr/bin/env python

#Copyright (C) 2013 by <NAME>
#
#Released under the MIT license, see LICENSE.txt
import sys
import os
import argparse

from pybedtools import BedTool, Interval
from teHmm.common import intersectSize, getLocalTempPath, runShellCommand
from teHmm.common import initBedTool, cleanBedTool

"""
Take the Bed output of ltr_finder, and make some different bed files that we can use for the model by changing the ids

output.bed : Uniuqe id's removed
output_sym.bed : and right and left removed
output_tsd_as_gap.bed : tsd states removed
output_sym_tsd_as_gap.bed : tsd states removed and right and left removed
output_tsd_as_ltr.bed : tsd states changed to ltr
output_sym_tsd_as_ltr.bed : tsd states changed to ltr right and left removed
output_single.bed : all annotated elems get same id 

"""

def main(argv=None):
    if argv is None:
        argv = sys.argv

    parser = argparse.ArgumentParser(
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
        description="Remove ltr_finder ids from 4th column")
    parser.add_argument("inBed", help="bed with ltr results to process")
    parser.add_argument("outBed", help="bed to write output to.  Will also "
                        "write outBed_sym.bed outBed_tsd_as_gap.bed etc.")
    parser.add_argument("--keepOl", help="by default, if LTR elements "
                        "overlap, the one with the highest score (length "
                        "in event of tie) is kept. This option disables"
                        " this logic.", action="store_true", default=False)
    
    args = parser.parse_args()
    tempBedToolPath = initBedTool()
    assert os.path.exists(args.inBed)
    baseOut, ext = os.path.splitext(args.outBed)

    inBed = args.inBed

    toRm = []
    if not args.keepOl:
        inBed = getLocalTempPath("Temp", ".bed")
        removeOverlaps(args.inBed, inBed)
        toRm.append(inBed)

    os.system("sed -e \"s/|LTR_TE|[0-9]*//g\" -e \"s/|-//g\" %s > %s" % (
        inBed, args.outBed))

    symBed = baseOut + "_sym" + ext
    os.system("sed -e \"s/|left//g\" -e \"s/|right//g\" %s > %s" % (args.outBed,
                                                                    symBed))

    tsd_as_gapsBed = baseOut + "_tsd_as_gap" + ext
    os.system("grep -v TSD %s > %s" % (args.outBed, tsd_as_gapsBed))

    sym_tsd_as_gapsBed = baseOut + "_sym_tsd_as_gap" + ext
    os.system("grep -v TSD %s > %s" % (symBed, sym_tsd_as_gapsBed))

    tsd_as_ltrBed = baseOut + "_tsd_as_ltr" + ext
    os.system("sed -e \"s/TSD/LTR/g\" %s > %s" % (args.outBed, tsd_as_ltrBed))

    sym_tsd_as_ltrBed = baseOut + "_sym_tsd_as_ltr" + ext
    os.system("sed -e \"s/TSD/LTR/g\" %s > %s" % (symBed, sym_tsd_as_ltrBed))

    singleBed = baseOut + "_single" + ext
    os.system("sed -e \"s/LTR/inside/g\" %s > %s" % (sym_tsd_as_ltrBed,
                                                     singleBed))

    for path in toRm:
        runShellCommand("rm -f %s" % path)
        
    cleanBedTool(tempBedToolPath)
        

def removeOverlaps(inBed, outBed):
    """ Little hack to get this script workign with different settings of ltr_finder
    where annotations can overlap.  To resolve overlaps, we choose the best element
    (by score, then length), and delete anything it touches.  TODO: incorporate Dougs
    script for his lastz stuff? """

    bedIntervals = [x for x in BedTool(inBed).sort()]
    outFile = open(outBed, "w")

    def getLtrID(interval):
        return interval.chrom + interval.name[interval.name.rfind("|") + 1:]
    
    # pass 1: element sizes
    sizes = dict()
    for interval in bedIntervals:
        id = getLtrID(interval)
        length = int(interval.end) - int(interval.start)
        if id in sizes:
            sizes[id] += length
        else:
            sizes[id] = length

    # pass 2: greedy kill (not optimal for all transitive cases)
    dead = set()
    for i, interval in enumerate(bedIntervals):
        id = getLtrID(interval)
        size = sizes[id]
        if id in dead:
            continue
        for j in xrange(i-1, -1, -1):
            if intersectSize((interval.chrom, interval.start, interval.end),
                             (bedIntervals[j].chrom, bedIntervals[j].start,
                             bedIntervals[j].end)) <= 0:
                break
            otherId = getLtrID(bedIntervals[j])
            if otherId not in dead and (
                    bedIntervals[j].score > interval.score or
                    (bedIntervals[j].score == interval.score and
                   sizes[otherId] > size)):
                dead.add(id)
                break
        if id in dead:
            continue
        for j in xrange(i+1, len(bedIntervals), 1):
            if intersectSize((interval.chrom, interval.start, interval.end),
                             (bedIntervals[j].chrom, bedIntervals[j].start,
                             bedIntervals[j].end)) <= 0:
                break
            otherId = getLtrID(bedIntervals[j])
            if otherId not in dead and (
                    bedIntervals[j].score > interval.score or
                    (bedIntervals[j].score == interval.score and
                   sizes[otherId] > size)):
                dead.add(id)
                break

    # pass 3: write non-killed
    for interval in bedIntervals:
        id = getLtrID(interval)
        if id not in dead:
            if interval.strand == "?":
                interval.strand = "."
            outFile.write(str(interval))
            
        
if __name__ == "__main__":
    sys.exit(main())<commit_msg>fix bug by making overlap test two sided<commit_after>#!/usr/bin/env python

#Copyright (C) 2013 by <NAME>
#
#Released under the MIT license, see LICENSE.txt
import sys
import os
import argparse

from pybedtools import BedTool, Interval
from teHmm.common import intersectSize, getLocalTempPath, runShellCommand
from teHmm.common import initBedTool, cleanBedTool

"""
Take the Bed output of ltr_finder, and make some different bed files that we can use for the model by changing the ids

output.bed : Uniuqe id's removed
output_sym.bed : and right and left removed
output_tsd_as_gap.bed : tsd states removed
output_sym_tsd_as_gap.bed : tsd states removed and right and left removed
output_tsd_as_ltr.bed : tsd states changed to ltr
output_sym_tsd_as_ltr.bed : tsd states changed to ltr right and left removed
output_single.bed : all annotated elems get same id 

"""

def main(argv=None):
    if argv is None:
        argv = sys.argv

    parser = argparse.ArgumentParser(
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
        description="Remove ltr_finder ids from 4th column")
    parser.add_argument("inBed", help="bed with ltr results to process")
    parser.add_argument("outBed", help="bed to write output to.  Will also "
                        "write outBed_sym.bed outBed_tsd_as_gap.bed etc.")
    parser.add_argument("--keepOl", help="by default, if LTR elements "
                        "overlap, the one with the highest score (length "
                        "in event of tie) is kept. This option disables"
                        " this logic.", action="store_true", default=False)
    parser.add_argument("--all", help="write _sym, _tsd_as_gap, etc. versions"
                        " of output", action="store_true", default=False)
    
    args = parser.parse_args()
    tempBedToolPath = initBedTool()
    assert os.path.exists(args.inBed)
    baseOut, ext = os.path.splitext(args.outBed)

    inBed = args.inBed

    toRm = []
    if not args.keepOl:
        inBed = getLocalTempPath("Temp", ".bed")
        removeOverlaps(args.inBed, inBed)
        toRm.append(inBed)

    os.system("sed -e \"s/|LTR_TE|[0-9]*//g\" -e \"s/|-//g\" %s > %s" % (
        inBed, args.outBed))

    if args.all:
        symBed = baseOut + "_sym" + ext
        os.system("sed -e \"s/|left//g\" -e \"s/|right//g\" %s > %s" % (args.outBed,
                                                                        symBed))

        tsd_as_gapsBed = baseOut + "_tsd_as_gap" + ext
        os.system("grep -v TSD %s > %s" % (args.outBed, tsd_as_gapsBed))

        sym_tsd_as_gapsBed = baseOut + "_sym_tsd_as_gap" + ext
        os.system("grep -v TSD %s > %s" % (symBed, sym_tsd_as_gapsBed))

        tsd_as_ltrBed = baseOut + "_tsd_as_ltr" + ext
        os.system("sed -e \"s/TSD/LTR/g\" %s > %s" % (args.outBed, tsd_as_ltrBed))

        sym_tsd_as_ltrBed = baseOut + "_sym_tsd_as_ltr" + ext
        os.system("sed -e \"s/TSD/LTR/g\" %s > %s" % (symBed, sym_tsd_as_ltrBed))

        singleBed = baseOut + "_single" + ext
        os.system("sed -e \"s/LTR/inside/g\" %s > %s" % (sym_tsd_as_ltrBed,
                                                         singleBed))

    for path in toRm:
        runShellCommand("rm -f %s" % path)
        
    cleanBedTool(tempBedToolPath)
        

def removeOverlaps(inBed, outBed):
    """ Little hack to get this script workign with different settings of ltr_finder
    where annotations can overlap.  To resolve overlaps, we choose the best element
    (by score, then length), and delete anything it touches.  TODO: incorporate Dougs
    script for his lastz stuff? """

    bedIntervals = [x for x in BedTool(inBed).sort()]
    outFile = open(outBed, "w")

    def getLtrID(interval):
        return interval.chrom + interval.name[interval.name.rfind("|") + 1:]
    
    # pass 1: element sizes
    sizes = dict()
    for interval in bedIntervals:
        id = getLtrID(interval)
        length = int(interval.end) - int(interval.start)
        if id in sizes:
            sizes[id] += length
        else:
            sizes[id] = length

    # pass 2: greedy kill (not optimal for all transitive cases)
    # strategy: any pairwise overlap will be detected in either
    # the left or right scan of at least one of the overlapping
    # elements. 
    dead = set()
    for i, interval in enumerate(bedIntervals):
        id = getLtrID(interval)
        size = sizes[id]
        if id in dead:
            continue
        for j in xrange(i-1, -1, -1):
            if intersectSize((interval.chrom, interval.start, interval.end),
                             (bedIntervals[j].chrom, bedIntervals[j].start,
                             bedIntervals[j].end)) <= 0:
                break
            otherId = getLtrID(bedIntervals[j])
            if otherId in dead:
                continue
            if (bedIntervals[j].score > interval.score or
                (bedIntervals[j].score == interval.score and
                 sizes[otherId] > size)):
                dead.add(id)
                break
            else:
                dead.add(otherId)
        if id in dead:
            continue
        for j in xrange(i+1, len(bedIntervals), 1):
            if intersectSize((interval.chrom, interval.start, interval.end),
                             (bedIntervals[j].chrom, bedIntervals[j].start,
                             bedIntervals[j].end)) <= 0:
                break
            otherId = getLtrID(bedIntervals[j])
            if otherId in dead:
                continue
            if (bedIntervals[j].score > interval.score or
                (bedIntervals[j].score == interval.score and
                 sizes[otherId] > size)):
                dead.add(id)
                break
            else:
                dead.add(otherId)
        if id in dead:
            continue

    # pass 3: write non-killed
    for interval in bedIntervals:
        id = getLtrID(interval)
        if id not in dead:
            if interval.strand == "?":
                interval.strand = "."
            outFile.write(str(interval))
            
        
if __name__ == "__main__":
    sys.exit(main())
{'fraction_non_alphanumeric': 0.1062383137854836, 'fraction_numerical': 0.0030596634370219276, 'mean_word_length': 2.0140881147540983, 'pattern_counts': {'":': 4, '<': 9, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '776033', 'n_tokens_mistral': 3837, 'n_tokens_neox': 3574, 'n_words': 1005}