#!/usr/bin/python
# File created on Nov 27 Jan 2012
from __future__ import division

__author__ = "Kishori M Konwar"
__copyright__ = "Copyright 2013, MetaPathways"
__credits__ = ["r"]
__version__ = "1.0"
__maintainer__ = "Kishori M Konwar"
__status__ = "Release"

try:
    from os import makedirs, sys, remove, rename
    from sys import path
    import re, traceback
    from optparse import OptionParser, OptionGroup
    from glob import glob
    from libs.python_modules.utils.metapathways_utils  import parse_command_line_parameters, fprintf, printf, eprintf
    from libs.python_modules.utils.metapathways_utils  import strip_taxonomy, ShortenORFId, ShortenContigId
    from libs.python_modules.utils.sysutil import getstatusoutput

except:
    print """ Could not load some user defined  module functions"""
    print """ Make sure your typed 'source MetaPathwaysrc'"""
    print """ """
    sys.exit(3)


usage=  sys.argv[0] + """ -d dbname1 -b parsed_blastout_for_database1 -w weight_for_database1 [-d dbname2 -b parsed_blastout_for_database2 -w weight_for_database2 ] [ --rRNA_16S  16SrRNA-stats-table ] [ --tRNA tRNA-stats-table ] [ --compact_output ]"""
parser = None

def createParser():
     global parser
     epilog = """Based on the parsed homology search results against reference protein databases, 
               tRNA scan results and matches against the LSU and SSU rRNA sequence databases
               annotations for individual ORFs are created.
               The functional annotation is done by synthesizing  the annotations from the 
               reference databases, after removing redundant texts and also making sure 
               the annotations from different reference databases are in agreement with each other.
               The taxonomic annotation is done by applying the LCA (lowest common ancestor rule) 
               on the hits from the RefSeq NR protein database."""

     epilog = re.sub(r'\s+',' ', epilog)
     parser = OptionParser(usage = usage, epilog = epilog)

     parser.add_option("-b", "--blastoutput", dest="input_blastout", action='append', default=[],
                       help='blastout files in TSV format [at least 1 REQUIRED]')
     parser.add_option("-a", "--algorithm", dest="algorithm", default="BLAST", help="algorithm BLAST or LAST" )

     parser.add_option("-m", "--contig-map", dest="contig_map_file", default=None, help="contig map file" )
     
     parser.add_option("-d", "--dbasename", dest="database_name", action='append', default=[],
                       help='the database names [at least 1 REQUIRED]')

     parser.add_option("-D", "--blastdir", dest="blastdir",  default=None,
                       help='the blast dir where all the BLAST outputs are located')
     
     parser.add_option("-s", "--samplename", dest="sample_name",  default=None,
                       help='the sample name')
     
     parser.add_option("-w", "--weight_for_database", dest="weight_db", action='append', default=[], type='float',
                       help='the map file for the database  [at least 1 REQUIRED]')
     
     parser.add_option( "--rRNA_16S", dest="rRNA_16S", action="append", default=[], 
                       help='the 16s rRNA stats file [OPTIONAL]')
     
     parser.add_option( "--tRNA", dest="tRNA", action="append", default=[], 
                       help='the tRNA stats file [OPTIONAL]')
     
     cutoffs_group =  OptionGroup(parser, 'Cuttoff Related Options')
     
     cutoffs_group.add_option("--min_score", dest="min_score", type='float', default=20,
                       help='the minimum bit score cutoff [default = 20 ] ')
     
     cutoffs_group.add_option("--max_evalue", dest="max_evalue", type='float', default=1e-6,
                       help='the maximum E-value cutoff [ default = 1e-6 ] ')
     c
utoffs_group.add_option("--min_length", dest=
"min_length", type='float', default=30,
                       help='the minimum length of query cutoff [default = 30 ] ')
     cutoffs_group.add_option("--max_length", dest="max_length", type='float', default=10000,
                       help='the maximum length of query cutoff [default = 10000 ] ')
     
     cutoffs_group.add_option("--min_identity", dest="min_identity", type='float', default=20,
                       help='the minimum identity of query cutoff [default 30 ] ')
     cutoffs_group.add_option("--max_identity", dest="max_identity", type='float', default=100,
                       help='the maximum identity of query cutoff [default = 100 ] ')
     
     cutoffs_group.add_option("--limit", dest="limit", type='float', default=5,
                       help='max number of hits per query cutoff [default = 5 ] ')
     
     cutoffs_group.add_option("--min_bsr", dest="min_bsr", type='float', default=0.00,
                       help='minimum BIT SCORE RATIO [default = 0.00 ] ')
     parser.add_option_group(cutoffs_group)
     
     
     output_options_group =  OptionGroup(parser, 'Output Options')
     output_options_group.add_option("--tax", dest="taxonomy", action='store_true', default=False,
                       help='add the taxonomy info [useful for refseq] ')
     parser.add_option_group(output_options_group)
     
     parser.add_option('-o' , "--output_gff", dest="output_gff",
                      help='the output gff file [REQUIRED]')
     
     parser.add_option('--output-comparative-annotation', dest="output_comparative_annotation",
                      help='the comparative output table [REQUIRED]')

     parser.add_option( "--compact_output", dest="compact_output", action='store_true', default=False,
                      help='compact output [OPTIONAL]')

     
     parser.add_option('--input_gff', dest='input_gff',
                metavar='INPUT', help='Unannotated gff file [REQUIRED]')



def check_arguments(opts, args):

    return True 

    if len(opts.input_blastout) == 0:
         eprintf("There should be at least one blastoutput file\n")  
         return False

    if len(opts.database_name) == 0:
         eprintf("There should be at least one database name\n")  
         return False

    if len(opts.weight_db) == 0:
         eprint("There should be at least one weight\n")  
         return False

    if len(opts.input_blastout) != len(opts.database_name) or\
         len(opts.input_blastout) !=  len(opts.weight_db) :
         eprint("The num of database names, blastoutputs and database map file should be equal\n")
         return False

    if opts.output_gff == None:
       eprintf("Must specify the output gff file\n")
       return False

    if opts.output_comparative_annotation == None:
       eprintf("Must specify the output tables for comparative annotation\n")
       return False

    if opts.input_gff == None:
       eprintf("Must specify the input gff file\n")
       return False

    return True

def insert_attribute(attributes, attribStr):
     rawfields = re.split('=', attribStr)
     if len(rawfields) == 2:
       if rawfields[0].strip().lower()=='id':
           orfid = ShortenORFId(rawfields[1].strip())
           attributes[rawfields[0].strip().lower()] = orfid
       else:
           attributes[rawfields[0].strip().lower()] = rawfields[1].strip()

def split_attributes(str, attributes):
     rawattributes = re.split(';', str)
     for attribStr in rawattributes:
        insert_attribute(attributes, attribStr)

     return attributes

def insert_orf_into_dict(line, contig_dict):
     rawfields = re.split('\t', line)

     fields = []
     for field in rawfields:
        fields.append(field.strip());
     
     if( len(fields) != 9):
       return

     attributes = {}

  
     seqname = fields[0]
     try:
       seqname = ShortenContigId(fields[0])
     except:
       seqname = fields[0]


     attributes['seqname'] =  seqname   # this is a bit of a  duplication  
     attributes['source'] =  fields[1]
     attributes['feature'] =  fields[2]
     attributes['start'] =  int(fields[3])
     attributes['end'] =  int(fields[4])

     try:
        attributes['score'] =  float(fields[5])
     except:
        attributes['score'] =  fields[5]

     attributes['strand'] =  fields[6]