###########################################################################
#                                                                         #
# C O P Y R I G H T   N O T I C E                                         #
#  Copyright (c) 2003-10 by:                                              #
#    * California Institute of Technology                                 #
#                                                                         #
#    All Rights Reserved.                                                 #
#                                                                         #
# Permission is hereby granted, free of charge, to any person             #
# obtaining a copy of this software and associated documentation files    #
# (the "Software"), to deal in the Software without restriction,          #
# including without limitation the rights to use, copy, modify, merge,    #
# publish, distribute, sublicense, and/or sell copies of the Software,    #
# and to permit persons to whom the Software is furnished to do so,       #
# subject to the following conditions:                                    #
#                                                                         #
# The above copyright notice and this permission notice shall be          #
# included in all copies or substantial portions of the Software.         #
#                                                                         #
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #
# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #
# SOFTWARE.                                                               #
###########################################################################
#
# data for Monodelphis domestica
import string
from cistematic.genomes import Genome
from os import environ

if environ.get("CISTEMATIC_ROOT"):
    cisRoot = environ.get("CISTEMATIC_ROOT") 
else:
    cisRoot = "/proj/genome"

geneDB = "%s/M_domestica/mdomestica.genedb" % cisRoot


def loadChromosome(db, chromPath, chromOutPath):
    seqArray = []
    seqLen = 0
    mdGenome = Genome("mdomestica", dbFile=db)
    inFile = open(chromPath, "r")
    header = inFile.readline()
    while header != "":
        seqArray = []
        seqLen = 0
        chromID = header.strip()[1:]
        currentLine = inFile.readline()
        while currentLine != "" and currentLine[0] != ">":
            lineSeq = currentLine.strip()
            seqLen += len(lineSeq)
            seqArray.append(lineSeq)
            currentLine = inFile.readline()

        seq = string.join(seqArray, "")
        if seqLen < 500000:
            print "Added contig %s to database" % chromID
            mdGenome.addSequence(("mdomestica", chromID), seq, "chromosome", str(seqLen))
            mdGenome.addChromosomeEntry(chromID, chromID, "db")
        else:
            outFileName = "%s%s.bin" % (chromOutPath, chromID)
            outFile = open("%s%s" % (cisRoot, outFileName), "w")
            outFile.write(seq)
            outFile.close()
            print "Added contig file %s to database" % outFileName
            mdGenome.addChromosomeEntry(chromID, outFileName, "file")

        header = currentLine

    inFile.close()


def loadGeneEntries(db, gFile):
    """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES
    """
    geneEntries = []
    mdGenome = Genome("mdomestica", dbFile=db)
    geneFile = open(gFile, "r")
    for line in geneFile:
        cols = line.split("\t")
        gid = cols[0]
        start = int(cols[5])
        stop = int(cols[6])
        sense = cols[2]
        chrom = cols[1]
        if sense == "+":
            sense = "F"
        else:
            sense = "R"

        geneID = ("mdomestica", gid)
        gidVersion = 1
        geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion))

    print "Adding %d gene entries" % len(geneEntries)
    mdGenome.addGeneEntryBatch(geneEntries)


def loadGeneAnnotations(db, annotPath):
    geneAnnotations = []
    annotFile = open(annotPath, "r")  
    mdGenome = Genome("mdomestica", dbFile=db)
    for line in annotFile:
        try:
            cols = line.split("\t")
            locID = cols[0]
            geneDesc = cols[6]
            if len(locID) > 0:
                geneAnnotations.append((("mdomestica", locID), string.replace(geneDesc.strip(), "'", "p")))
        except:
            pass

    print "Adding %d annotations" % len(geneAnnotations)
    mdGenome.addAnnotationBatch(geneAnnotations)


def loadGeneFeatures(db, gfile):
    geneFile = open(gfile, "r")
    senseArray = {"+": "F",
                  "-": "R",
                  ".": "F"
    }

    seenArray = []
    insertArray = []
    for geneLine in geneFile:
        geneFields = geneLine.split("\t")
        exonNum = int(geneFields[7])
        exonStarts = geneFields[8].split(",")
        exonStops = geneFields[9].split(",")
        chrom = geneFields[1]
        sense = senseArray[geneFields[2]]
        gstop = int(geneFields[6]) - 1
        gstart = int(geneFields[5]) - 1
        geneid = geneFields[0]
        try:
            geneID = ("mdomestica", geneid)
        except:
            continue

        gidVersion = "1"
        if geneID in seenArray:
            gidVersion = "2" # doesn't deal with more than 2 refseq's for the same locus, yet.
        else:
            seenArray.append(geneID)

        for index in range(exonNum):
            estart = int(exonStarts[index]) - 1
            estop = int(exonStops[index]) - 1
            if estart >= gstart and estop <= gstop:
                insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, "CDS"))
            elif estop <= gstart:
                if sense == "F":
                    fType = "5UTR"
                else:
                    fType = "3UTR"

                insertArray.append((geneID, 1, chrom, estart, estop, sense, fType))
            elif estart >= gstop:
                if sense == "F":
                    fType = "3UTR"
                else:
                    fType = "5UTR"

                insertArray.append((geneID, 1, chrom, estart, estop, sense, fType))
            elif estart <= gstop:
                if sense == "F":
                    fType = "3UTR"
                else:
                    fType = "5UTR"

                insertArray.append((geneID, 1, chrom, estart, gstop, sense, "CDS"))
                insertArray.append((geneID, 1, chrom, gstop + 1, estop, sense, fType))
            else:
                if sense == "F":
                    fType = "5UTR"
                else:
                    fType = "3UTR"

                insertArray.append((geneID, 1, chrom, gstart, estop, sense, "CDS"))
                insertArray.append((geneID, 1, chrom, estart, gstart - 1, sense, fType))

    geneFile.close()
    mdGenome = Genome("mdomestica", dbFile=db)
    print "Adding %d features" % len(insertArray)
    mdGenome.addFeatureEntryBatch(insertArray)


def createDBFile(db):
    mdGenome = Genome("mdomestica",  dbFile=db)
    mdGenome.createGeneDB(db)


def createDBindices(db):
    mdGenome = Genome("mdomestica", dbFile=db)
    mdGenome.createIndices()


def buildMdomesticaDB(db=geneDB):
    genePath = "%s/download/mondom/genscan.txt" % cisRoot
    chromoPath = "%s/download/mondom/softMask.fa" % cisRoot
    chromoOutPath = "/M_domestica/"

    print "Creating database %s" % db
    createDBFile(db)

    print "Adding gene entries"
    loadGeneEntries(db, genePath)

    print "Adding gene features"
    loadGeneFeatures(db, genePath)

    print "Loading chromosomes"
    loadChromosome(db, chromoPath, chromoOutPath)

    print "Creating Indices"
    createDBindices(db)

    print "Finished creating database %s" % db