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

#Copyright (C) 2014 by <NAME>
#
#Released under the MIT license, see LICENSE.txt
import sys
import os
import argparse
import logging
import numpy as np
import math
import copy

from teHmm.track import TrackList
from teHmm.trackIO import readTrackData
from teHmm.common import myLog, EPSILON, initBedTool, cleanBedTool
from teHmm.common import addLoggingOptions, setLoggingFromOptions, logger
from teHmm.common import runShellCommand, getLogLevelString, getLocalTempPath

"""

The HMM cannot use annotation tracks as output by several tools (ex repeatmasker)
without first doing some name-munging and setting some scaling parameters for binning.

This script takes as input a list of "raw" annotation tracks, and runs all the necessary scripts to produce a list of "clean" tracks that can be used by the HMM.

NOTE: This script is really hardcoded to run only on ./mustang_alyrata_tracks.xml at the moment, and needs to be updated to reflect changes to that file (maybe).  A more general workflow will need to be put in place later....
"""

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

    parser = argparse.ArgumentParser(
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
        description="Generate HMM-usable tracklist from raw tracklist. EX "
        "used to transform mustang_alyrata_tracks.xml -> "
        "mustang_alyrata_clean.xml.  Runs cleanChaux.py cleanLtrFinder.py and "
        " cleanTermini.py and addTsdTrack.py and setTrackScaling.py (also runs "
        " removeBedOverlaps.py before each of the clean scripts)")
    
    parser.add_argument("tracksInfo", help="Path of Tracks Info file "
                        "containing paths to genome annotation tracks")
    parser.add_argument("allBed", help="Bed file spanning entire genome")
    parser.add_argument("cleanTrackPath", help="Directory to write cleaned BED"
                        " tracks to")
    parser.add_argument("outTracksInfo", help="Path to write modified tracks XML"
                        " to.")
    parser.add_argument("--numBins", help="Maximum number of bins after scaling",
                        default=10, type=int)
    parser.add_argument("--scaleTracks", help="Comma-separated list of tracks "
                        "to process for scaling. If not set, all"
                        " tracks listed as having a multinomial distribution"
                        " (since this is the default value, this includes "
                        "tracks with no distribution attribute) will be"
                        " processed.", default=None)
    parser.add_argument("--skipScale", help="Comma-separated list of tracks to "
                        "skip for scaling.", default=None)
    parser.add_argument("--chaux", help="Name of chaux track", default="chaux")
    parser.add_argument("--ltrfinder", help="comma-sep Name(s) of ltrfinder track(s)",
                        default="ltr_finder,ltr_harvest")
    parser.add_argument("--termini", help="Name of termini track",
                        default="termini")
    parser.add_argument("--sequence", help="Name of fasta sequence track",
                        default="sequence")
    parser.add_argument("--tsd", help="Name of tsd track to generate",
                        default="tsd")
    parser.add_argument("--tir", help="Name of tir_termini track",
                        default="tir_termini")
    parser.add_argument("--hollister", help="Name of hollister track",
                        default="hollister")
    parser.add_argument("--repbase", help="Name of repbase track",
                        default="repbase")
    parser.add_argument("--repeat_modeler", help="Name of repeat_modeler track",
                        default="repeat_modeler")
    parser.add_argument("--noScale", help="Dont do any scaling", default=False,
                        action="store_true")
    parser.add_argument("--noTsd", help="Dont generate TSD track.  NOTE:"
                        " TSD track is hardcoded to be generated from "
                        "termini and (non-LTR elements of ) chaux",
                        default=False, action="store_true")
    
    addLoggingOptions(parser)
    args = parser.parse_args()
    setLoggingFromOptions(args)
    tempBedToolPath = initBedTool()
    args.logOpString = "--logLevel %s" % getLogLevelString()
    if args.logFile is not None:
        args.logOpString += " --logFile %s" % args.logFile

    try:
        os.makedirs(args.cleanTrackPath)
    except:
        pass
    if not os.path.isdir(args.cleanTrackPath):
        raise RuntimeError("Unable to find or create cleanTrack dir %s" %
                           args.cleanTrackPath)

    tempTracksInfo = getLocalTempPath("Temp_mustang_alyrata_clean", "xml")
    runCleaning(args, tempTracksInfo)
    assert os.path.isfile(tempTracksInfo)

    runTsd(args, tempTracksInfo)
    
    runScaling(args, tempTracksInfo)

    runShellCommand("rm -f %s" % tempTracksInfo)

    cleanBedTool(tempBedToolPath)

def cleanPath(args, track):
    """ path of cleaned track """
    oldPath = track.getPath()
    oldFile = os.path.basename(oldPath)
    oldName, oldExt = os.path.splitext(oldFile)
    return os.path.join(args.cleanTrackPath, oldName + "_clean" + ".bed")
    
def runCleaning(args, tempTracksInfo):
    """ run scripts for cleaning chaux, ltr_finder, and termini"""
    trackList = TrackList(args.tracksInfo)

    # run cleanChaux.py on chaux track
    chauxTrack = trackList.getTrackByName(args.chaux)
    if chauxTrack is not None:
        inFile = chauxTrack.getPath()
        outFile = cleanPath(args, chauxTrack)
        tempBed = getLocalTempPath("Temp_chaux", ".bed")
        runShellCommand("cleanChaux.py --keepUnderscore %s > %s" % (inFile,
                                                                    tempBed))
        runShellCommand("removeBedOverlaps.py %s > %s" % (tempBed, outFile))
        runShellCommand("rm -f %s" % tempBed)
        chauxTrack.setPath(outFile)
    else:
        logger.warning("Could not find chaux track")

    # run cleanChaux.py on hollister track
    hollisterTrack = trackList.getTrackByName(args.hollister)
    if hollisterTrack is not None:
        inFile = hollisterTrack.getPath()
        outFile = cleanPath(args, hollisterTrack)
        tempBed = getLocalTempPath("Temp_hollister", ".bed")
        runShellCommand("cleanChaux.py %s > %s" % (inFile, tempBed))
        runShellCommand("removeBedOverlaps.py %s > %s" % (tempBed, outFile)) 
        runShellCommand("rm -f %s" % tempBed)
        hollisterTrack.setPath(outFile)
    else:
        logger.warning("Could not find hollister track")

    # run cleanChaux.py on repbase track
    repbaseTrack = trackList.getTrackByName(args.repbase)
    if repbaseTrack is not None:
        inFile = repbaseTrack.getPath()
        outFile = cleanPath(args, repbaseTrack)
        tempBed = getLocalTempPath("Temp_repbase", ".bed")
        runShellCommand("cleanChaux.py %s > %s" % (inFile, tempBed))
        runShellCommand("removeBedOverlaps.py %s > %s" % (tempBed, outFile)) 
        runShellCommand("rm -f %s" % tempBed)
        repbaseTrack.setPath(outFile)
    else:
        logger.warning("Could not find repbase track")

    # run cleanChaux.py on repeat_modeler track
    repeat_modelerTrack = trackList.getTrackByName(args.repeat_modeler)
    if repeat_modelerTrack is not None:
        inFile = repeat_modelerTrack.getPath()
        outFile = cleanPath(args, repeat_modelerTrack)
        tempBed1 = None<commit_msg>add back repbase threshold so work on lyrata<commit_after>#!/usr/bin/env python

#Copyright (C) 2014 by <NAME>
#
#Released under the MIT license, see LICENSE.txt
import sys
import os
import argparse
import logging
import numpy as np
import math
import copy

from teHmm.track import TrackList
from teHmm.trackIO import readTrackData
from teHmm.common import myLog, EPSILON, initBedTool, cleanBedTool
from teHmm.common import addLoggingOptions, setLoggingFromOptions, logger
from teHmm.common import runShellCommand, getLogLevelString, getLocalTempPath

"""

The HMM cannot use annotation tracks as output by several tools (ex repeatmasker)
without first doing some name-munging and setting some scaling parameters for binning.

This script takes as input a list of "raw" annotation tracks, and runs all the necessary scripts to produce a list of "clean" tracks that can be used by the HMM.

NOTE: This script is really hardcoded to run only on ./mustang_alyrata_tracks.xml at the moment, and needs to be updated to reflect changes to that file (maybe).  A more general workflow will need to be put in place later....
"""

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

    parser = argparse.ArgumentParser(
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
        description="Generate HMM-usable tracklist from raw tracklist. EX "
        "used to transform mustang_alyrata_tracks.xml -> "
        "mustang_alyrata_clean.xml.  Runs cleanChaux.py cleanLtrFinder.py and "
        " cleanTermini.py and addTsdTrack.py and setTrackScaling.py (also runs "
        " removeBedOverlaps.py before each of the clean scripts)")
    
    parser.add_argument("tracksInfo", help="Path of Tracks Info file "
                        "containing paths to genome annotation tracks")
    parser.add_argument("allBed", help="Bed file spanning entire genome")
    parser.add_argument("cleanTrackPath", help="Directory to write cleaned BED"
                        " tracks to")
    parser.add_argument("outTracksInfo", help="Path to write modified tracks XML"
                        " to.")
    parser.add_argument("--numBins", help="Maximum number of bins after scaling",
                        default=10, type=int)
    parser.add_argument("--scaleTracks", help="Comma-separated list of tracks "
                        "to process for scaling. If not set, all"
                        " tracks listed as having a multinomial distribution"
                        " (since this is the default value, this includes "
                        "tracks with no distribution attribute) will be"
                        " processed.", default=None)
    parser.add_argument("--skipScale", help="Comma-separated list of tracks to "
                        "skip for scaling.", default=None)
    parser.add_argument("--chaux", help="Name of chaux track", default="chaux")
    parser.add_argument("--ltrfinder", help="comma-sep Name(s) of ltrfinder track(s)",
                        default="ltr_finder,ltr_harvest")
    parser.add_argument("--termini", help="Name of termini track",
                        default="termini")
    parser.add_argument("--sequence", help="Name of fasta sequence track",
                        default="sequence")
    parser.add_argument("--tsd", help="Name of tsd track to generate",
                        default="tsd")
    parser.add_argument("--tir", help="Name of tir_termini track",
                        default="tir_termini")
    parser.add_argument("--hollister", help="Name of hollister track",
                        default="hollister")
    parser.add_argument("--repbase", help="Name of repbase track",
                        default="repbase")
    parser.add_argument("--repeat_modeler", help="Name of repeat_modeler track",
                        default="repeat_modeler")
    parser.add_argument("--noScale", help="Dont do any scaling", default=False,
                        action="store_true")
    parser.add_argument("--noTsd", help="Dont generate TSD track.  NOTE:"
                        " TSD track is hardcoded to be generated from "
                        "termini and (non-LTR elements of ) chaux",
                        default=False, action="store_true")
    
    addLoggingOptions(parser)
    args = parser.parse_args()
    setLoggingFromOptions(args)
    tempBedToolPath = initBedTool()
    args.logOpString = "--logLevel %s" % getLogLevelString()
    if args.logFile is not None:
        args.logOpString += " --logFile %s" % args.logFile

    try:
        os.makedirs(args.cleanTrackPath)
    except:
        pass
    if not os.path.isdir(args.cleanTrackPath):
        raise RuntimeError("Unable to find or create cleanTrack dir %s" %
                           args.cleanTrackPath)

    tempTracksInfo = getLocalTempPath("Temp_mustang_alyrata_clean", "xml")
    runCleaning(args, tempTracksInfo)
    assert os.path.isfile(tempTracksInfo)

    runTsd(args, tempTracksInfo)
    
    runScaling(args, tempTracksInfo)

    runShellCommand("rm -f %s" % tempTracksInfo)

    cleanBedTool(tempBedToolPath)

def cleanPath(args, track):
    """ path of cleaned track """
    oldPath = track.getPath()
    oldFile = os.path.basename(oldPath)
    oldName, oldExt = os.path.splitext(oldFile)
    return os.path.join(args.cleanTrackPath, oldName + "_clean" + ".bed")
    
def runCleaning(args, tempTracksInfo):
    """ run scripts for cleaning chaux, ltr_finder, and termini"""
    trackList = TrackList(args.tracksInfo)

    # run cleanChaux.py on chaux track
    chauxTrack = trackList.getTrackByName(args.chaux)
    if chauxTrack is not None:
        inFile = chauxTrack.getPath()
        outFile = cleanPath(args, chauxTrack)
        tempBed = getLocalTempPath("Temp_chaux", ".bed")
        runShellCommand("cleanChaux.py --keepUnderscore %s > %s" % (inFile,
                                                                    tempBed))
        runShellCommand("removeBedOverlaps.py %s > %s" % (tempBed, outFile))
        runShellCommand("rm -f %s" % tempBed)
        chauxTrack.setPath(outFile)
    else:
        logger.warning("Could not find chaux track")

    # run cleanChaux.py on hollister track
    hollisterTrack = trackList.getTrackByName(args.hollister)
    if hollisterTrack is not None:
        inFile = hollisterTrack.getPath()
        outFile = cleanPath(args, hollisterTrack)
        tempBed = getLocalTempPath("Temp_hollister", ".bed")
        runShellCommand("cleanChaux.py %s > %s" % (inFile, tempBed))
        runShellCommand("removeBedOverlaps.py %s > %s" % (tempBed, outFile)) 
        runShellCommand("rm -f %s" % tempBed)
        hollisterTrack.setPath(outFile)
    else:
        logger.warning("Could not find hollister track")

    # run cleanChaux.py on repbase track
    repbaseTrack = trackList.getTrackByName(args.repbase)
    if repbaseTrack is not None:
        inFile = repbaseTrack.getPath()
        outFile = cleanPath(args, repbaseTrack)
        tempBed = getLocalTempPath("Temp_repbase", ".bed")
        runShellCommand("cleanChaux.py %s --minScore 0 > %s" % (inFile, tempBed))
        runShellCommand("removeBedOverlaps.py %s > %s" % (tempBed, outFile)) 
        runShellCommand("rm -f %s" % tempBed)
        repbaseTrack.setPath(outFile)
    else:
        logger.warning("Could not find repbase track")

    # run cleanChaux.py on repeat_modeler track
    repeat_modelerTrack = trackList.getTrackByName(args.repeat_modeler)
    if repeat_modelerTrack is not None:
        inFile = repeat_modelerTrack.getPath()
        outFile = cleanPath(args, repeat_modelerTrack)
        tempBed1 = None
{'fraction_non_alphanumeric': 0.09720390398311791, 'fraction_numerical': 0.0009891849116328145, 'mean_word_length': 2.726044226044226, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 8}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '334997', 'n_tokens_mistral': 4468, 'n_tokens_neox': 4278, 'n_words': 1202}