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

#Copyright (C) 2013 by <NAME> (<EMAIL>)
#
#Released under the MIT license, see LICENSE.txtimport unittest

#Wed Apr 10 15:30:53 PDT 2013
#
#Generates necessary files to make assembly hub
#Input: 1/ Output directory (e.g ~nknguyen/public_html/ecoli/hub/testhubs)
#       2/ hal file of the multiple alignment
#       4/ (Optional: directory containing annotated bed files. e.g : genes)
#Output:
#   outdir/
#       hub.txt
#       genomes.txt
#

import os, sys, re
from optparse import OptionParser
from sonLib.bioio import system  #(may need to remove this dependency on sonLib)

def getLongestSeq(seq2len):
    seqs = sorted( [(seq, len) for seq, len in seq2len.iteritems()], key=lambda item:item[1], reverse=True )
    return seqs[0]

def getGenomeSequencesFromHal(halfile, genome):
    statsfile = "%s-seqStats.txt" %genome
    system("halStats --sequenceStats %s %s > %s" %(genome, halfile, statsfile))
    
    seq2len = {}
    f = open(statsfile, 'r')
    for line in f:
        if len(line) < 2 or re.search("SequenceName", line):
            continue
        items = line.strip().split(", ")
        seq = items[0].split('.')[-1]
        #seq = items[0]
        l = int(items[1])
        seq2len[seq] = l
    f.close()
    system("rm %s" %statsfile)

    return seq2len

def getGenomeSequences(halfile, genomes):
    genome2seq2len = {}
    for genome in genomes:
        seq2len = getGenomeSequencesFromHal(halfile, genome)
        if len(seq2len) == 0:
            sys.stderr.write("Warning: genome %s contains 0 sequence - no browser was made.\n" %genome)
        else:
            genome2seq2len[genome] = seq2len
    return genome2seq2len

def getChromSizes(halfile, genome2seq2len, outdir):
    #Get chrom.sizes for all genomes
    #Output to outdir/genomeChr.sizes for each genome
    for genome, seq2len in genome2seq2len.iteritems():
        outfile = os.path.join(outdir, genome)
        f = open(outfile, 'w')
        for s, l in seq2len.iteritems():
            if l > 0:
                f.write("%s\t%d\n" %(s, l))
        f.close()

def getGenomesFromHal(halfile):
    #Get a list of all genomes from the output of halStats
    statsfile = "halStats.txt"
    system("halStats --genomes %s > %s" %(halfile, statsfile))
    
    f = open(statsfile, 'r')
    genomes = f.readline().strip().split(',')
    f.close()

    #clean up
    system("rm %s" %statsfile)
    
    return genomes

def liftoverBedFiles(indir, halfile, genome2seq2len, outdir):
    #beddir has the hierachy: indir/genome/chr1.bed, chr2.bed...
    #for each genome in beddir, lifeover the bed records of that genome to the coordinate of all other genomes

    #Get chr.sizes for all genomes:
    chrsizedir = os.path.join(outdir, "chrsizes")
    system("mkdir -p %s" %chrsizedir)
    getChromSizes(halfile, genome2seq2len, chrsizedir)

    #liftover bed file of each genome with availabel beds to all genomes
    genomes = genome2seq2len.keys()
    for genome in os.listdir(indir):
        if genome not in genomes:
            continue
        genomeindir = os.path.join(indir, genome)
        assert os.path.isdir(genomeindir)

        #Create bed directory for current genome
        genomeoutdir = os.path.join(outdir, genome)
        system("mkdir -p %s" %genomeoutdir)

        #Concatenate all the input bed files and convert it into bigbed to outdir/genome/genome.bb
        tempbed = "%s-temp.bed" % os.path.join(genomeoutdir, genome)
        system( "cat %s/*bed > %s" %(genomeindir, tempbed) )
        system( "bedSort %s %s" % (tempbed, tempbed) )

        outbigbed = os.path.join(genomeoutdir, "%s.bb" %genome) 
        chrsizefile = os.path.join(chrsizedir, genome)
        system( "bedToBigBed %s %s %s" %(tempbed, chrsizefile, outbigbed) )

        #Liftover to all other genomes:
        for othergenome in genomes:
            if othergenome == genome:
                continue
            liftovertempbed = "%s.bed" % os.path.join(genomeoutdir, othergenome)
            system("halLiftover %s %s %s %s %s" %(halfile, genome, tempbed, othergenome, liftovertempbed))
            if re.search("Anc", othergenome): #HACK
                system("awk ' $0 !~ /#/ {split($1, arr, \".\"); print arr[2] \"\t\" $2 \"\t\" $3 \"\t\" $4} ' %s > %s-reformat" %(liftovertempbed, liftovertempbed))
                system("mv %s-reformat %s" %(liftovertempbed, liftovertempbed))

            system("bedSort %s %s" %(liftovertempbed, liftovertempbed))
            outbigbed = os.path.join(genomeoutdir, "%s.bb" %othergenome)
            chrsizefile = os.path.join(chrsizedir, othergenome)
            system("bedToBigBed %s %s %s" %(liftovertempbed, chrsizefile, outbigbed))
            #Cleanup:
            system("rm %s" % liftovertempbed)
        system("rm %s" %tempbed) #cleanup
    #Cleanup
    system("rm -Rf %s" %chrsizedir)

def makeTwoBitSeqFile(genome, halfile, outdir):
    fafile = os.path.join(outdir, "%s.fa" %genome)
    system("hal2fasta --outFaPath %s %s %s" %(fafile, halfile, genome))
    
    #if sequence headers have "." (e.g genome.chr), reformat the header to only have "chr"
    fafile2 = "%s2" %fafile
    cmd = "awk '{ if($0 ~/>/){split($1, arr, \".\"); if(length(arr) > 1 ){print \">\" arr[2]}else{print $0} }else{ print $0} }' %s > %s" %(fafile, fafile2)
    system(cmd)
    system("rm %s" %fafile)

    #convert to 2bit files
    twobitfile = os.path.join(outdir, "%s.2bit" %genome)
    system("faToTwoBit %s %s" %(fafile2, twobitfile))
    system("rm %s" %fafile2)

def writeTrackDb_bigbeds(f, bigbeddir, genomes, currgenome):
    for genome in os.listdir(bigbeddir):
        if genome == 'chrsizes':
            continue
        f.write("track gene%s\n" % genome)
        f.write("longLabel %s Liftovered Genes\n" % genome)
        f.write("shortLabel %sgenes\n" % genome)
        f.write("bigDataUrl ../%s\n" % os.path.join( os.path.basename(bigbeddir), genome, "%s.bb" % currgenome ) )
        f.write("type bigBed\n")
        if genome == currgenome: 
            f.write("visibility dense\n")
        else:
            f.write("visibility hide\n")
        f.write("\n")

def writeTrackDb_snakes(f, halfile, genomes, currgenome):
    for genome in genomes:
        if re.search(genome, currgenome): #current genome
            continue
        #SNAKE TRACKS
        f.write("track snake%s\n" %genome)
        f.write("longLabel %s Snake\n" %genome)
        f.write("shortLabel %s\n" %genome)
        f.write("otherSpecies %s\n" %genome)
        f.write("visibility full\n")
        f.write("bigDataUrl %s\n" % halfile)
        f.write("type halSnake\n")
        f.write("\n")

def writeTrackDbFile(genomes, halfile, outdir, bigbeddir):
    currgenome = outdir.rstrip('/').split("/")[-1]
    filename = os.path.join(outdir, "trackDb.txt")
    f = open(filename, 'w')
    
    if bigbeddir:
        writeTrackDb_bigbeds(f, bigbeddir, genomes, currgenome)
    
    writeTrackDb_snakes(f, halfile, genomes, currgenome)
    f.close()

def writeDescriptionFile(genome, outdir):
    filename = os.path.join(outdir, "description.html")
    f = open(filename, 'w')
    f.write("%s\n" %genome)
    f.close()
    return

def fixLodFilePath(lodtxtfile, localHalfile, outdir):
    #fix the path of the original hal file to point to the created
    #link relative to the output directory
    relPath = os.path.relpath(localHalfile, start=outdir)
    lodTxtBuf = ''
    for line in open(lodtxtfile):
        tokens = line.split()
        if len(tokens) == 2 and tokens[0] == '0':
            lodTxtBuf += '0 %s\n' % relPath
        else:
            lodTxtBuf += line
    with open(lodtxtfile, 'w') as lodFile:
        lodFile.write(lodTxtBuf)
    
def getLodFiles(localHalfile, options, outdir):
    lodtxtfile = os.path.join(outdir, "lod.txt") #outdir/lod.txt
    loddir = os.path.join(outdir, "lod") #outdir/lod
    if options.lodtxtfile and options.loddir: #if lod files were given, then just make soft links to them
        if os.path.exists(lodtxtfile):
            if os.path.abspath(lodtxtfile) != os.path.abspath(options.lodtxtfile):
                system("rm %s" %lodtxtfile)
                system("ln -s %s %s" %(os.path.abspath(options.lodtxtfile), lodtxtfile))
        else:
            system("ln -s %s %s" %(os.path.abspath(options.lodtxtfile), lodtxtfile))

        if os.path.exists(loddir):
            if os.path.abspath(loddir) != os.path.abspath(options.loddir):
                if os.path.islink(loddir):
                    system("rm %s" %loddir)
                else:
                    system("rm -Rf %s" %loddir)
                loddir = os.path.join(outdir, os.path.basename(options.loddir))
                system("ln -s %s %s" %(os.path.abspath(options.loddir), loddir))
        else:
            system("ln -s %s %s" %(os.path.abspath(options.loddir), loddir))
    else: #if lod files were not given, create them using halLodInterpolate.py
        system("halLodInterpolate.py %s %s --outHalDir %s %s" %(localHalfile, lodtxtfile, loddir, options.lodOpts))
        fixLodFilePath(lodtxtfile, localHalfile, outdir)
    return lodtxtfile, loddir

def writeGenomesFile(genome2seq2len, halfile, options, outdir):
    '''Write genome for all samples in hal file
    '''

    localHalfile = os.path.join(outdir, os.path.basename(halfile))
    if os.path.abspath(localHalfile) != os.path.abspath(halfile):
        if os.path.exists(localHalfile):
            system("rm %s" %localHalfile)
        system("ln -s %s %s" %(os.path.abspath(halfile), localHalfile))

    #Create lod files if useLod is specified
    lodtxtfile = ''
    loddir = ''
    options.lodOpts = ''
    if options.lodMaxBlock is not None:
        options.lodOpts += '--maxBlock %d ' % options.lodMaxBlock
    if options.lodScale is not None:
        options.lodOpts += '--scale %f ' % options.lodScale
    if options.lodMaxDNA is not None:
        options.lodOpts += '--maxDNA %d ' % options.lodMaxDNA
    if options.lodInMemory is True:
        options.lodOpts += '--inMemory'
    if len(options.lodOpts) > 0:
        options.lod = True
    if options.lod:
        lodtxtfile, loddir = getLodFiles(localHalfile, options, outdir)
    <commit_msg>always use --trans for halLodInterpolate.py<commit_after>#!/usr/bin/env python

#Copyright (C) 2013 by <NAME> (<EMAIL>)
#
#Released under the MIT license, see LICENSE.txtimport unittest

#Wed Apr 10 15:30:53 PDT 2013
#
#Generates necessary files to make assembly hub
#Input: 1/ Output directory (e.g ~nknguyen/public_html/ecoli/hub/testhubs)
#       2/ hal file of the multiple alignment
#       4/ (Optional: directory containing annotated bed files. e.g : genes)
#Output:
#   outdir/
#       hub.txt
#       genomes.txt
#

import os, sys, re
from optparse import OptionParser
from sonLib.bioio import system  #(may need to remove this dependency on sonLib)

def getLongestSeq(seq2len):
    seqs = sorted( [(seq, len) for seq, len in seq2len.iteritems()], key=lambda item:item[1], reverse=True )
    return seqs[0]

def getGenomeSequencesFromHal(halfile, genome):
    statsfile = "%s-seqStats.txt" %genome
    system("halStats --sequenceStats %s %s > %s" %(genome, halfile, statsfile))
    
    seq2len = {}
    f = open(statsfile, 'r')
    for line in f:
        if len(line) < 2 or re.search("SequenceName", line):
            continue
        items = line.strip().split(", ")
        seq = items[0].split('.')[-1]
        #seq = items[0]
        l = int(items[1])
        seq2len[seq] = l
    f.close()
    system("rm %s" %statsfile)

    return seq2len

def getGenomeSequences(halfile, genomes):
    genome2seq2len = {}
    for genome in genomes:
        seq2len = getGenomeSequencesFromHal(halfile, genome)
        if len(seq2len) == 0:
            sys.stderr.write("Warning: genome %s contains 0 sequence - no browser was made.\n" %genome)
        else:
            genome2seq2len[genome] = seq2len
    return genome2seq2len

def getChromSizes(halfile, genome2seq2len, outdir):
    #Get chrom.sizes for all genomes
    #Output to outdir/genomeChr.sizes for each genome
    for genome, seq2len in genome2seq2len.iteritems():
        outfile = os.path.join(outdir, genome)
        f = open(outfile, 'w')
        for s, l in seq2len.iteritems():
            if l > 0:
                f.write("%s\t%d\n" %(s, l))
        f.close()

def getGenomesFromHal(halfile):
    #Get a list of all genomes from the output of halStats
    statsfile = "halStats.txt"
    system("halStats --genomes %s > %s" %(halfile, statsfile))
    
    f = open(statsfile, 'r')
    genomes = f.readline().strip().split(',')
    f.close()

    #clean up
    system("rm %s" %statsfile)
    
    return genomes

def liftoverBedFiles(indir, halfile, genome2seq2len, outdir):
    #beddir has the hierachy: indir/genome/chr1.bed, chr2.bed...
    #for each genome in beddir, lifeover the bed records of that genome to the coordinate of all other genomes

    #Get chr.sizes for all genomes:
    chrsizedir = os.path.join(outdir, "chrsizes")
    system("mkdir -p %s" %chrsizedir)
    getChromSizes(halfile, genome2seq2len, chrsizedir)

    #liftover bed file of each genome with availabel beds to all genomes
    genomes = genome2seq2len.keys()
    for genome in os.listdir(indir):
        if genome not in genomes:
            continue
        genomeindir = os.path.join(indir, genome)
        assert os.path.isdir(genomeindir)

        #Create bed directory for current genome
        genomeoutdir = os.path.join(outdir, genome)
        system("mkdir -p %s" %genomeoutdir)

        #Concatenate all the input bed files and convert it into bigbed to outdir/genome/genome.bb
        tempbed = "%s-temp.bed" % os.path.join(genomeoutdir, genome)
        system( "cat %s/*bed > %s" %(genomeindir, tempbed) )
        system( "bedSort %s %s" % (tempbed, tempbed) )

        outbigbed = os.path.join(genomeoutdir, "%s.bb" %genome) 
        chrsizefile = os.path.join(chrsizedir, genome)
        system( "bedToBigBed %s %s %s" %(tempbed, chrsizefile, outbigbed) )

        #Liftover to all other genomes:
        for othergenome in genomes:
            if othergenome == genome:
                continue
            liftovertempbed = "%s.bed" % os.path.join(genomeoutdir, othergenome)
            system("halLiftover %s %s %s %s %s" %(halfile, genome, tempbed, othergenome, liftovertempbed))
            if re.search("Anc", othergenome): #HACK
                system("awk ' $0 !~ /#/ {split($1, arr, \".\"); print arr[2] \"\t\" $2 \"\t\" $3 \"\t\" $4} ' %s > %s-reformat" %(liftovertempbed, liftovertempbed))
                system("mv %s-reformat %s" %(liftovertempbed, liftovertempbed))

            system("bedSort %s %s" %(liftovertempbed, liftovertempbed))
            outbigbed = os.path.join(genomeoutdir, "%s.bb" %othergenome)
            chrsizefile = os.path.join(chrsizedir, othergenome)
            system("bedToBigBed %s %s %s" %(liftovertempbed, chrsizefile, outbigbed))
            #Cleanup:
            system("rm %s" % liftovertempbed)
        system("rm %s" %tempbed) #cleanup
    #Cleanup
    system("rm -Rf %s" %chrsizedir)

def makeTwoBitSeqFile(genome, halfile, outdir):
    fafile = os.path.join(outdir, "%s.fa" %genome)
    system("hal2fasta --outFaPath %s %s %s" %(fafile, halfile, genome))
    
    #if sequence headers have "." (e.g genome.chr), reformat the header to only have "chr"
    fafile2 = "%s2" %fafile
    cmd = "awk '{ if($0 ~/>/){split($1, arr, \".\"); if(length(arr) > 1 ){print \">\" arr[2]}else{print $0} }else{ print $0} }' %s > %s" %(fafile, fafile2)
    system(cmd)
    system("rm %s" %fafile)

    #convert to 2bit files
    twobitfile = os.path.join(outdir, "%s.2bit" %genome)
    system("faToTwoBit %s %s" %(fafile2, twobitfile))
    system("rm %s" %fafile2)

def writeTrackDb_bigbeds(f, bigbeddir, genomes, currgenome):
    for genome in os.listdir(bigbeddir):
        if genome == 'chrsizes':
            continue
        f.write("track gene%s\n" % genome)
        f.write("longLabel %s Liftovered Genes\n" % genome)
        f.write("shortLabel %sgenes\n" % genome)
        f.write("bigDataUrl ../%s\n" % os.path.join( os.path.basename(bigbeddir), genome, "%s.bb" % currgenome ) )
        f.write("type bigBed\n")
        if genome == currgenome: 
            f.write("visibility dense\n")
        else:
            f.write("visibility hide\n")
        f.write("\n")

def writeTrackDb_snakes(f, halfile, genomes, currgenome):
    for genome in genomes:
        if re.search(genome, currgenome): #current genome
            continue
        #SNAKE TRACKS
        f.write("track snake%s\n" %genome)
        f.write("longLabel %s Snake\n" %genome)
        f.write("shortLabel %s\n" %genome)
        f.write("otherSpecies %s\n" %genome)
        f.write("visibility full\n")
        f.write("bigDataUrl %s\n" % halfile)
        f.write("type halSnake\n")
        f.write("\n")

def writeTrackDbFile(genomes, halfile, outdir, bigbeddir):
    currgenome = outdir.rstrip('/').split("/")[-1]
    filename = os.path.join(outdir, "trackDb.txt")
    f = open(filename, 'w')
    
    if bigbeddir:
        writeTrackDb_bigbeds(f, bigbeddir, genomes, currgenome)
    
    writeTrackDb_snakes(f, halfile, genomes, currgenome)
    f.close()

def writeDescriptionFile(genome, outdir):
    filename = os.path.join(outdir, "description.html")
    f = open(filename, 'w')
    f.write("%s\n" %genome)
    f.close()
    return

def fixLodFilePath(lodtxtfile, localHalfile, outdir):
    #fix the path of the original hal file to point to the created
    #link relative to the output directory
    relPath = os.path.relpath(localHalfile, start=outdir)
    lodTxtBuf = ''
    for line in open(lodtxtfile):
        tokens = line.split()
        if len(tokens) == 2 and tokens[0] == '0':
            lodTxtBuf += '0 %s\n' % relPath
        else:
            lodTxtBuf += line
    with open(lodtxtfile, 'w') as lodFile:
        lodFile.write(lodTxtBuf)
    
def getLodFiles(localHalfile, options, outdir):
    lodtxtfile = os.path.join(outdir, "lod.txt") #outdir/lod.txt
    loddir = os.path.join(outdir, "lod") #outdir/lod
    if options.lodtxtfile and options.loddir: #if lod files were given, then just make soft links to them
        if os.path.exists(lodtxtfile):
            if os.path.abspath(lodtxtfile) != os.path.abspath(options.lodtxtfile):
                system("rm %s" %lodtxtfile)
                system("ln -s %s %s" %(os.path.abspath(options.lodtxtfile), lodtxtfile))
        else:
            system("ln -s %s %s" %(os.path.abspath(options.lodtxtfile), lodtxtfile))

        if os.path.exists(loddir):
            if os.path.abspath(loddir) != os.path.abspath(options.loddir):
                if os.path.islink(loddir):
                    system("rm %s" %loddir)
                else:
                    system("rm -Rf %s" %loddir)
                loddir = os.path.join(outdir, os.path.basename(options.loddir))
                system("ln -s %s %s" %(os.path.abspath(options.loddir), loddir))
        else:
            system("ln -s %s %s" %(os.path.abspath(options.loddir), loddir))
    else: #if lod files were not given, create them using halLodInterpolate.py
        system("halLodInterpolate.py %s %s --outHalDir %s %s" %(localHalfile, lodtxtfile, loddir, options.lodOpts))
        fixLodFilePath(lodtxtfile, localHalfile, outdir)
    return lodtxtfile, loddir

def writeGenomesFile(genome2seq2len, halfile, options, outdir):
    '''Write genome for all samples in hal file
    '''

    localHalfile = os.path.join(outdir, os.path.basename(halfile))
    if os.path.abspath(localHalfile) != os.path.abspath(halfile):
        if os.path.exists(localHalfile):
            system("rm %s" %localHalfile)
        system("ln -s %s %s" %(os.path.abspath(halfile), localHalfile))

    #Create lod files if useLod is specified
    lodtxtfile = ''
    loddir = ''
    options.lodOpts = ''
    if options.lodMaxBlock is not None:
        options.lodOpts += '--maxBlock %d ' % options.lodMaxBlock
    if options.lodScale is not None:
        options.lodOpts += '--scale %f ' % options.lodScale
    if options.lodMaxDNA is not None:
        options.lodOpts += '--maxDNA %d ' % options.lodMaxDNA
    if options.lodInMemory is True:
        options.lodOpts += '--inMemory '
    if len(options.lodOpts) > 0:
        options.lodOpts += '--trans '
        options.lod = True
    if options.lod:
        lodtxtfile, loddir = getLodFiles(localHalfile, options, outdir)
    
{'fraction_non_alphanumeric': 0.13523688456461155, 'fraction_numerical': 0.008311739109177138, 'mean_word_length': 2.6721723518850986, 'pattern_counts': {'":': 0, '<': 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': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1969139', 'n_tokens_mistral': 7373, 'n_tokens_neox': 6996, 'n_words': 1774}