#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# iMembrane : Annotation database constructor
#
# Author: Sebastian Kelm
# Created: 10/05/2010
#


import sys
import os
import subprocess
from glob import glob

sys.path = [os.path.join(os.path.abspath(os.path.dirname(sys.argv[0]))), 
            os.path.join(os.path.abspath(os.path.dirname(sys.argv[0])), "..", "lib-python"),
            os.path.join(os.path.abspath(os.path.dirname(sys.argv[0])), "..", "..", "lib-python")] + sys.path

from prosci.util.pdb import Pdb
from prosci.util.ali import Ali
from prosci.shell import Params

from prosci.imembrane.builder import Classifier, Builder, DEFAULT_TIMELIMIT, DEFAULT_SURFACELIMIT, DEFAULT_GRANULARITY



params = Params()

if params.opts or len(params.args) < 2 or len(params.args) > 5:
    sys.stderr.write("""
iMembrane - database building component

Input: CGDB-style annotation
Output: JOY-style annotation, Z-coordinate list for membrane layers, b-factor annotated structures

JOY-style annotation types:
"membrane contact"   Residues are either in contact with the lipid tails
                 (T) or the polar head groups (H) or neither (N).
"membrane layer"     Residues are either in the lipid tail layer (T) or
                 in a peripheral polar head group layer (H) or not
                 within the membrane at all (N).

Residues that cannot be annotated are marked with question marks (?).


USAGE:
"""+params.scriptname+""" cgdb_directory imembranedb_directory [time [surface [granularity]]]
\n""")
    sys.exit(1)


DBDIR  = os.path.abspath(params.args[0])
OUTDIR = os.path.abspath(params.args[1])
timeLimit, surfaceLimit, granularity = params.getArgs(2, 3, [DEFAULT_TIMELIMIT, DEFAULT_SURFACELIMIT, DEFAULT_GRANULARITY], [float]*3)

ENTRYDIR = "%s/entries"%(OUTDIR)
SEQDIR   = "%s/seq"%(OUTDIR)



if not os.path.exists(DBDIR):
    sys.stderr.write("ERROR: Database directory does not exist: %s\n" % (DBDIR))
    sys.exit(2)

if not os.path.exists("%s"%(OUTDIR)):
    os.mkdir("%s"%(OUTDIR))



if not os.path.exists(ENTRYDIR):
    os.mkdir(ENTRYDIR)

for entrystructure in glob("%s/chains/*.atm"%(DBDIR)):
    entrycode = os.path.splitext(os.path.basename(entrystructure))[0]
    
    if os.path.exists("%s/%s/annotation.tem"%(ENTRYDIR, entrycode)):
      continue
    
    classifier = Classifier(DBDIR, entrycode)
    sequence = classifier.pdb_data.get_seq()
    ali_file = Ali(">%s\nstructure\n%s*\n"%(entrycode, sequence))
    
    a = Builder(ali_file, classifier)
    
    a.make_all_basic_annotation(timeLimit, surfaceLimit, granularity)
    
    a.write_output(ENTRYDIR)



if not os.path.exists(SEQDIR):
    os.mkdir(SEQDIR)

os.chdir(SEQDIR)
subprocess.call('cat "%s/"*/sequence.fasta > chains.fasta'%(ENTRYDIR), shell=True)
subprocess.call("grep '>' chains.fasta > chains.list", shell=True)
subprocess.call("makeblastdb -in chains.fasta -out chains -dbtype prot", shell=True)
subprocess.call('chmod a+rx "%s" "%s" "%s"' % (OUTDIR, ENTRYDIR, SEQDIR), shell=True)
subprocess.call('chmod a+r "%s/"*' % (SEQDIR), shell=True)
