delta2js.py

#!/usr/bin/env python
"""
delta2js.py

Converts mummer .delta files to javascript database for webmummer.
https://github.com/sauloal/webmummer

Saulo Aflitos
Plant Research International - PRI
Cluster Bioinformatics
Wageningen University and Research Centre - WUR
Centre for BioSystems Genomics - CBSG

To Run:
    cd data
    rm <delta folder>/*.js
    delta2js.py <delta folder>/*.delta
    cp <delta folder>/*.js .

Documentation generated by:
    pydoc -w delta2js

License: GPL 2.0
"""

__author__  = "Saulo Aflitos"
__date__    = "Nov 2013"
__version__ = "2013_11_27_14_06"
__credits__ = "mummer, WUR, PRI, CBSG"

import sys, os
import subprocess
import re
import time
import datetime
from pprint import pprint as pp

#cd data
#rm ../../../out/*.js
#../delta2js.py ../../../out/*.delta
#cp -l ../../../out/*.js .
#rm ../webmummer.tar.xz; tar --exclude .git -ahcvf ../webmummer.tar.xz .


dry_run     = None
compulsory  = None
forbidden   = None
labelFields = None
titleFmt    = None
xlabelFmt   = None
ylabelFmt   = None
statusMatch = None

if __name__ == '__main__':
    if len(sys.argv) < 3:
        print "not enought arguments"
        sys.exit( 1 )

    config = sys.argv[ 1 ]

    if not os.path.exists( config ):
        print "config file %s does not exists" % config
        sys.exit( 1 )

    if not config.endswith( '.py' ):
        print "config file %s not a python script" % config
        sys.exit( 1 )



    execfile( config )



    for name, var in [
            ['dry_run'    , dry_run    ],
            ['compulsory' , compulsory ],
            ['forbidden'  , forbidden  ],
            ['labelFields', labelFields],
            ['titleFmt'   , titleFmt   ],
            ['xlabelFmt'  , xlabelFmt  ],
            ['ylabelFmt'  , ylabelFmt  ],
            ['statusMatch', statusMatch]
        ]:
        if var is None:
            print "variable %s not defined in config %s" % ( name, config )



    for k in labelFields:
        labelFields[k][0] = re.compile( labelFields[k][0] )
    """
    pre-compiles RE, replacing the original string
    """


    def statusMatcher( status ):
        """
        Given the extension of the file, return a verbose string describing it to be used in title.
        """
        try:
            return statusMatch[ status ]
        except:
            print "status %s does not have a name" % status
            sys.exit(1)





    labelFields['status'][2]= statusMatcher
    """
    Adds statusMatcher as a function to be called uppon parsing of file name.
    """



#pointFmt = "[{x:%d,y:%d},{x:%d,y:%d},{n:%d,s:%d,q:%.2f}]," # db format 1.0
#pointFmt = "[%d,%d,%d,%d,%d,%d,%.2f],"                     # db format 2.0
pointFmt = "%d,%d,%d,%d,%d,%d,%.2f"                         # db format 3.0
"""
All numbers (x1, x2, y1, y2, ref name, tgt name, quality) are printed in a single array.
By knowing the array size (7), each register can be recovered by simple arithmetics (k * 7).
"""

#solanum_lycopersicum_heinz_SL2.40ch12.fa_._solanum_pennellii_scaffold_final.assembly.fasta.delta.q.delta.fplot
#solanum_lycopersicum_heinz_SL2.40ch12.fa_._solanum_pennellii_scaffold_final.assembly.fasta.delta.q.delta.rplot

#-- reverse hits sorted by %sim
#0 0 0
#0 0 0


#18542978 36303894 83.8077999012671
#18549032 36297865 83.8077999012671


class exp(object):
    """
    Exporter Class: exports database headers (title, axis labels and limits) and, uppon
    receiving each coordinates, append it to the array.
    """
    def __init__(self, outfile, refName, refChrom, tgtName, tgtChrom, status, title="title", xlabel="xlabel", ylabel="ylabel"):
        """
        Accepts output file, title and axis labels.
        Initialize limits.
        """
        self.title   = title
        self.xlabel  = xlabel
        self.ylabel  = ylabel
        self.outf    = outfile

        self.minX    = sys.maxint
        self.maxX    = 0

        self.minY    = sys.maxint
        self.maxY    = 0

        self.tgts    = {}

        self.fhd     = open(outfile, 'w')
        self.fhd2    = open(outfile+'.gff3', 'w')

        self.linec   = 0;

        self.dbreg   = "_filelist[ '%s' ][ '%s' ][ '%s' ][ '%s' ][ '%s' ]" % ( refName, refChrom, tgtName, tgtChrom, status )

        self.gffReg   = {}
        self.gffNames = []

        #outfiles[ refName ][ chromNumber ][ spp ][ status ]
        self.fhd.write("""\

%(dbreg)s[ 'title'  ] = '%(title)s';
%(dbreg)s[ 'xlabel' ] = '%(xlabel)s';
%(dbreg)s[ 'ylabel' ] = '%(ylabel)s';

%(dbreg)s[ 'points' ] = [""" % { 'dbreg': self.dbreg, 'title': self.title, 'xlabel': self.xlabel, 'ylabel': self.ylabel } )

        self.fhd2.write("##gff-version 3\n")


    def append(self, x1, y1, x2, y2, name, refName, sense, q):
        """
        Adds query name to database of names (a list of all names).
        Exports only species ID (index of the scaffold name in the database of names);

        Replaces sense for 0/1;

        Updates min and max positions;
        """

        self.linec += 1

        minx = x1 if (x1 < x2) else x2
        maxx = x1 if (x1 > x2) else x2

	print "minx %d  maxx %d x1 %d x2 %d" % (minx, maxx, x1, x2)
	assert minx < maxx, "minx (%d) > maxx(%d) for %s" % (minx, maxx, name)

        if name not in self.tgts:
            self.tgts[name] = len(self.tgts)

            self.gffReg[name] = []
            self.gffReg[name].append( [ refName, '.', 'mRNA', minx, maxx, '.', '+'                     , '.', 'ID=%s;Name=%s' % (refName+'_'+name,  name) ] )

            if name not in self.gffNames:
                self.gffNames.append( name )


        gffId   = "%s_frag_%05d" % ( refName+'_'+name, len( self.gffReg[ name ] ) )
        gffName = "%s_frag_%05d" % (             name, len( self.gffReg[ name ] ) )
        #self.gffReg[name].append( [refName, '.', 'CDS', minx, maxx, '.', '+' if sense == 'fwd' else '-', len( self.gffReg[ name ] )-1, 'ID=%s;Name=%s;Parent=%s' % (gffId, gffName, refName+'_'+name) ] )
        self.gffReg[name].append( [refName, '.', 'CDS', minx, maxx, '.', '+' if sense == 'fwd' else '-', 0, 'ID=%s;Name=%s;Parent=%s' % (gffId, gffName, refName+'_'+name) ] )
        self.gffReg[name][0][4] = maxx

        name  = self.tgts[name]

        sense = 0 if sense == 'fwd' else 1

        line  = pointFmt % ( x1, y1, x2, y2, name, sense, q )

        if self.linec != 1:
            line = ',' + line

        self.fhd.write( line )

        if x1 > self.maxX: self.maxX = x1
        if x2 > self.maxX: self.maxX = x2

        if x1 < self.minX: self.minX = x1
        if x2 < self.minX: self.minX = x2

        if y1 > self.maxY: self.maxY = y1
        if y2 > self.maxY: self.maxY = y2

        if y1 < self.minY: self.minY = y1
        if y2 < self.minY: self.minY = y2

    def add(self, sense, reg):
        """
        Filters from register only the necessary information.
        Calls self.append;
        """
        #print "adding", reg
        refStart, refEnd, tgtStart, tgtEnd, targetStart, targetEnd, refLen, tgtLen, refSub, tgtSub, idd, refName, tgtName = reg

        self.append( refStart, targetStart, refEnd, targetEnd, tgtName, refName, sense, idd )

    def close(self):
        """
        Exports min/max;
        Exports query species list;
        Closes filehandle;
        """
        line = """\
];

%(dbreg)s[ 'xmin'   ]  = %(minX)12d;
%(dbreg)s[ 'xmax'   ]  = %(maxX)12d;

%(dbreg)s[ 'ymin'   ]  = %(minY)12d;
%(dbreg)s[ 'ymax'   ]  = %(maxY)12d;
%(dbreg)s[ 'tgts'   ]  = [""" % { 'dbreg': self.dbreg, 'minX': self.minX, 'maxX': self.maxX, 'minY': self.minY, 'maxY': self.maxY }

        for scaf in sorted(self.tgts, key=lambda p: self.tgts[p]):
            line += "'%s'," % scaf

        line  = line.strip(',')
        line += '];'

        self.fhd.write( line )

        if len(self.gffReg) > 0:
            for key in self.gffNames:
                for line in self.gffReg[key]:
                    self.fhd2.write("\t".join( [ str(x) for x in line ] ) + "\n")


        self.fhd.close()
        self.fhd2.close()



def parseDelta(delta):
    """
    Converts .delta to .coords (if needed) and parses .coords file.
    """
    coordslines = []
    scafOrder   = []
    scafLens    = {}
    lineCount   = 0

    if not delta.endswith('.coords'):
        show_coords = "show-coords -l -r -T %s" % (delta)
        print show_coords
        coords      = subprocess.Popen(show_coords, stdout=subprocess.PIPE, shell=True).communicate()[0]
        coordslines = coords.split("\n")
    else:
        coordslines = open(delta, 'r')



    for line in coordslines:
        lineCount += 1
        line       = line.strip()
        #print line

        if len(line) == 0  : continue
        if lineCount < 5   : continue
        if line[0] == '='  : continue

        cols = line.split("\t")
        #print cols
        refStart    , refEnd    , tgtStart, tgtEnd, refLen, tgtLen = [ int(  x)  for x in cols[ 0: 6] ]
        idd                                                        =   float(cols[6])
        referenceLen, targetLen                                    = [ int(  x)  for x in cols[ 7: 9] ]
        refName     , tgtName                                      = [ x.strip() for x in cols[ 9:11] ]

        refSub = refEnd - refStart
        tgtSub = tgtEnd - tgtStart

        #print "tgt name  %-15s start %9d end %9d len %9d size %9d sub %9d" % ( tgtName, tgtStart, tgtEnd, tgtLen, targetLen   , tgtSub )
        #print "ref name  %-15s start %9d end %9d len %9d size %9d sub %9d" % ( refName, refStart, refEnd, refLen, referenceLen, refSub )
        #print "identity  %.2f"                                             % ( idd                                                     )
        #print

        refMin = min( [refStart, refEnd] )
        refMax = max( [refStart, refEnd] )
        tgtMin = min( [tgtStart, tgtEnd] )
        tgtMax = max( [tgtStart, tgtEnd] )
        d      = [refStart, refEnd, tgtStart, tgtEnd, refLen, tgtLen, refSub, tgtSub, idd, refName, tgtName]
        #print d

        if tgtName in scafLens:
            prevRefStart = scafLens[tgtName][0]
            prevRefEnd   = scafLens[tgtName][1]

            if refMin < prevRefStart:
                scafLens[tgtName][0] = refMin

            if refMax > prevRefEnd:
                scafLens[tgtName][1] = refMax


            prevTgtStart = scafLens[tgtName][3]
            prevTgtEnd   = scafLens[tgtName][4]

            if tgtMin < prevTgtStart:
                scafLens[tgtName][3] = tgtMin

            if tgtMax > prevTgtEnd:
                scafLens[tgtName][4] = tgtMax

            scafLens[tgtName][6].append(d)

        else:
            #                    0       1       2             3       4       5          6
            scafLens[tgtName] = [refMin, refMax, referenceLen, tgtMin, tgtMax, targetLen, [d]]

    scafOrder = sorted(scafLens.keys(), key=lambda s: scafLens[s][0])
    #print scafOrder


    if delta.endswith('.coords'):
        coordslines.close()


    return ( scafOrder, scafLens )



def genCoords(exporter, scafOrder, scafLens):
    """
    Check sanity of data and adds data to exporter.
    Converts query coordnates to Y coordinates where scaffolds are ordered and their lengths sum
    to create a single Y axis.
    """

    targetPos = 1
    for tgtName in scafOrder:
        #print tgtName
        refMin, refMax, referenceLen, tgtMin, tgtMax, targetLen, coordsData = scafLens[tgtName]
        #refSize = refMax - refMin

        tgtSize = tgtMax - tgtMin

        if tgtSize > targetLen:
            print "error"

        alignmentEnd = targetPos + targetLen

        for line in coordsData:
            refStart, refEnd, tgtStart, tgtEnd, refLen, tgtLen, refSub, tgtSub, idd, refName, tgtName = line
            #print line
            targetStart = targetPos + tgtStart
            targetEnd   = targetPos + tgtEnd

            plotStr    = """%d %d %.2f\n%d %d %.2f\n\n\n""" %   (
                    targetStart, refStart, idd,
                    targetEnd  , refEnd  , idd,
                )

            #print targetEnd, targetStart, alignmentEnd

            if targetEnd   > alignmentEnd:
                print " error"

            if targetStart > alignmentEnd:
                print " error"

            if abs(tgtSub) > abs(refSub * 1.1):
                #print " target size %d > ref size %d"  % (tgtSub, refSub)
                pass

            if (abs(tgtSub) > targetLen) or (abs(refSub) > (targetLen*1.1)):
                print "incongruency"
                print "tgt name  %-15s start %9d end %9d len %9d size %9d sub %9d" % ( tgtName, tgtStart, tgtEnd, tgtLen, targetLen   , tgtSub )
                print "ref name  %-15s start %9d end %9d len %9d size %9d sub %9d" % ( refName, refStart, refEnd, refLen, referenceLen, refSub )
                print "identity  %.2f"                                             % ( idd                                                     )
                print "tgtSub", abs(tgtSub), "targetLen   ", targetLen
                print "refSub", abs(refSub), "referenceLen", referenceLen
                print plotStr
                print
                sys.exit(1)


            #print "tgt name  %-15s start %9d end %9d len %9d size %9d sub %9d" % ( tgtName, tgtStart, tgtEnd, tgtLen, targetLen   , tgtSub )
            #print "ref name  %-15s start %9d end %9d len %9d size %9d sub %9d" % ( refName, refStart, refEnd, refLen, referenceLen, refSub )
            #print "identity  %.2f"                                             % ( idd                                                     )
            #print "tgtSub", abs(tgtSub), "targetLen   ", targetLen
            #print "refSub", abs(refSub), "referenceLen", referenceLen
            #print plotStr

            reg = ( refStart, refEnd, tgtStart, tgtEnd, targetStart, targetEnd, refLen, tgtLen, refSub, tgtSub, idd, refName, tgtName )

            if tgtSub > 0: # fwd
                exporter.add('fwd', reg )

            else: # rev
                exporter.add('rev', reg )

        targetPos += targetLen + 1

    print "last pos", targetPos



def parseFN(infile):
    """
    Parses file name to extract species information
    """
    bn = os.path.basename( infile )
    #solanum_lycopersicum_heinz_SL2.40ch00.fa_._solanum_arcanum_scaffold_final.assembly.fasta.delta

    print "parsing %s" % infile

    res = {}

    for label in sorted( labelFields ):
        #print 'parsing %s label %s' % ( infile, label )
        fmt = labelFields[label][0]
        bfr = labelFields[label][1]
        aft = labelFields[label][2]
        #print 'parsing %s label %s fmt %s' % ( infile, label, fmt.pattern )
        try:
            val = fmt.search(infile).group(1)
            #print 'parsing %s label %s val %s' % ( infile, label, str(val) )

            if bfr is not None:
                val = bfr( val )
                #print 'parsing %s label %s val %s b' % ( infile, label, str(val) )

            val = val.replace('_', ' ')
            #print 'parsing %s label %s val %s u' % ( infile, label, str(val) )

            if aft is not None:
                val = aft( val )
                #print 'parsing %s label %s val %s a' % ( infile, label, str(val) )

            #print 'parsing %s label %s val %s GOT' % ( infile, label, str(val) )

            res[label] = val
        except:
            print "error parsing %s in file name: %s using format %s" % (label, infile, fmt.pattern)
            sys.exit(1)

    return res



def main():
    """
    Parse all input files names, exports them to databases and creates index of all files.
    """
    outfiles   = {}

    refsNames  = {}
    refsChroms = {}
    tgtsNames  = {}
    tgtsChroms = {}
    statuses   = {}

    if len(sys.argv) == 1:
        print "Usage:"
        print sys.argv[0], ' [input delta files +]'
        sys.exit(0)



    filecount = 0
    numfiles  = len(sys.argv[2:])
    for infile in sys.argv[2:]:
        filecount += 1
        print "INFILE : %s [%3d/%3d]" % ( infile, filecount, numfiles )

        hasC = False
        for compul in compulsory:
            if infile.endswith( compul ):
                hasC = True

        if not hasC:
            print "  does not have compulsory extension. skipping", compulsory
            continue

        hasF = False
        for forb in forbidden:
            if infile.endswith( forb ):
                hasF = True

        if hasF:
            print "  has compulsory extension", forbidden
            continue


        outfile  = infile + '.js'
        print "OUTFILE: %s" % outfile

        labels = parseFN( infile )


        refName  = labels['refName' ]
        refChrom = labels['refChrom']

        tgtName  = labels['tgtName' ]
        tgtChrom = labels['tgtChrom']

        status   = labels['status'  ]



        refsNames  [ refName  ] = 1
        refsChroms [ refChrom ] = 1

        tgtsNames  [ tgtName  ] = 1
        tgtsChroms [ tgtChrom ] = 1

        statuses   [ status   ] = 1



        title  = titleFmt  % labels
        xlabel = xlabelFmt % labels
        ylabel = ylabelFmt % labels

        print labels

        if not dry_run:
            if not os.path.exists( outfile ):
                print "parsing delta"
                scafOrder, scafLens = parseDelta(infile)


                print "parsing coords"
                exporter = exp(outfile, refName, refChrom, tgtName, tgtChrom, status, title=title, xlabel=xlabel, ylabel=ylabel)
                genCoords(exporter, scafOrder, scafLens)
                exporter.close()



        if refName not in outfiles:
            outfiles[ refName ] = {}

        if refChrom not in outfiles[ refName ]:
            outfiles[ refName ][ refChrom ] = {}

        if tgtName not in outfiles[ refName ][ refChrom ]:
            outfiles[ refName ][ refChrom ][ tgtName ] = {}

        if tgtChrom not in outfiles[ refName ][ refChrom ][ tgtName ]:
            outfiles[ refName ][ refChrom ][ tgtName ][ tgtChrom ] = {}

        if status not in outfiles[ refName ][ refChrom ][ tgtName ][ tgtChrom ]:
            outfiles[ refName ][ refChrom ][ tgtName ][ tgtChrom ][ status ] = [
                [os.path.abspath(infile ), os.path.basename(infile )],
                [os.path.abspath(outfile), os.path.basename(outfile)],
            ]
        else:
            print "data read twice"
            print "ref %s ref chrom %s tgt %s tgt chrom %s status %s" % ( refName, refChrom, tgtName, tgtChrom, status )
            sys.exit( 1 )

        #break

    print
    print "RESUME"
    print " refs"
    pp(refsNames )

    print " refs chroms"
    pp(refsChroms)

    print " tgts"
    pp(tgtsNames )

    print " tgts chroms"
    pp(tgtsChroms)

    print " statuses"
    pp(statuses  )


    if dry_run:
        #return
        pass


    with open( 'list.js', 'w' ) as fhd:
        """
        Creates a list of all available status, query species, chromosomes and references for quick listing.
        Creates a database of the available combitations and the .js files containing the data.
        outfiles[ refName ][ refChrom ][ tgtName ][ tgtChrom ][ status ] = <javascript database base name>
        """

        datestr = datetime.datetime.fromtimestamp( time.time() ).strftime('%Y_%m_%d_%H_%M_%S')

        fhd.write( 'var _db_domain   = "%s";\n' % project_id);

        fhd.write( 'var _db_version  = "%s";\n' % datestr);

        refStr       = ', '.join( [ "'%s'" % x for x in sorted(refsNames ) ] )
        fhd.write( 'var _refsNames   = [%s];\n' % refStr);

        refChromStr  = ', '.join( [ "'%s'" % x for x in sorted(refsChroms) ] )
        fhd.write( 'var _refsChroms  = [%s];\n' % refChromStr);

        tgtStr       = ', '.join( [ "'%s'" % x for x in sorted(tgtsNames ) ] )
        fhd.write( 'var _tgtsNames   = [%s];\n' % tgtStr);

        tgtChromStr  = ', '.join( [ "'%s'" % x for x in sorted(tgtsChroms) ] )
        fhd.write( 'var _tgtsChroms  = [%s];\n' % tgtChromStr);

        statusStr    = ', '.join( [ "'%s'" % x for x in sorted(statuses  ) ] )
        fhd.write( 'var _statuses    = [%s];\n' % statusStr);

        fhd.write( 'var _filelist    = {\n');


        refNames = sorted(outfiles)
        for v, refName in enumerate(refNames):
            fhd.write("  '%s': {\n" % refName )

            refChroms = sorted(outfiles[ refName ])
            for w, refChrom in enumerate(refChroms):
                fhd.write("    '%s': {\n" % refChrom )

                tgtNames = sorted(outfiles[ refName ][ refChrom ])
                for x, tgtName in enumerate(tgtNames):
                    fhd.write("      '%s': {\n" % tgtName )

                    tgtChroms = sorted(outfiles[ refName ][ refChrom ][ tgtName ])
                    for y, tgtChrom in enumerate(tgtChroms):
                        fhd.write("        '%s': {\n" % tgtChrom )

                        statuses = sorted(outfiles[ refName ][ refChrom ][ tgtName ][ tgtChrom ])
                        for z, status in enumerate(statuses):
                            filedata = outfiles[ refName ][ refChrom ][ tgtName ][ tgtChrom ][ status ]
                            fhd.write("          '%s': {\n" % status)
                            fhd.write("                  'filename': '%s'\n" % ( filedata[1][1] ) )
                            fhd.write('          }')
                            if z < len(statuses)-1:
                                fhd.write(',')
                            fhd.write('\n')

                        fhd.write('        }')
                        if y < len(tgtChroms)-1:
                            fhd.write(',')
                        fhd.write('\n')

                    fhd.write('      }')
                    if x < len(tgtNames)-1:
                        fhd.write(',')
                    fhd.write('\n')

                fhd.write('    }')
                if w < len(refChroms)-1:
                    fhd.write(',')
                fhd.write('\n')

            fhd.write('  }')
            if v < len(refNames)-1:
                fhd.write(',')
            fhd.write('\n')

        fhd.write('};')

    print "done"



if __name__ == '__main__': main()