###############################################################################
#
# ptdomain.py - object to represent protein domains and functions to
#               parse into domains from external programs
#
# File:    ptdomain.py
# Author:  Alex Stivala
# Created: September 2007
#
# $Id: ptdomain.py 2011 2008-10-30 01:54:20Z astivala $
#
# PTDomain is a class representing a protein domain. A domain is represented
# by a list of segments, which are contiguo
us subsequences of a chain.
#
# Functions are provided to parse domains using different domain parsing
# programs and return the corresponding list of PTDomain objects.
# Supported so far is:
#    . DDOMAIN (Zhou et al 2007 Protein Science 16:947-955) program
#    . CATH (CATH Domall File (CDF) 2.0) file
#
###############################################################################

import os,sys

from Bio.PDB import * # only needed for DDomain when segment spans chai
ns
from ptutils import cleanup_tmpdir,get_int_icode
from ptsecstruct import PTSecStruct

#-----------------------------------------------------------------------------
#
# Module globals
#
#-----------------------------------------------------------------------------

verbose = False

#-----------------------------------------------------------------------------
#
# Class definitions 
#
#-----------------------------------------------------------------------------


#
# Empty classes for exceptions
#

class NotInCATH_Exception(Exception): # identifier not found in CATH CDF file
    pass

#
# Real classes
#

class PTSegment:
    """
    The PTSegment object represents a segment (contiguous subsequence of
    a chain) by the chain identifier, and start and end residue sequence
    numbers

    Note residue numbers here are integers, not PDB residue
    'numbers' which are actually strings and may have insertion
    codes, are not sequential (may have gaps, run backwards, etc.);
    unlike in PTNode, PTGraph2 etc. where we store the PDB residue numbers,
    (but use dictionary to get sequence integers when necesary).
    Using integers makes everything much
    simpler, can do calculations easily in domeval.py etc.
    For this to work strictly correctly, should use proper purely sequential
    numbering in these classes (PTSegment, PTDomain), such as that assigned
    by DSSP or STRIDE, or equivalently by using the index into list of
    residues built by Polypeptide Builder or similar from Bio.PDB.
    """
    def __init__(self, chainid, start_resnum, end_resnum):
        """
        Construct segment with supplied chain identifier and start and
        end residue sequence numbers.

        Parameters:
             chainid - PDB chain identifier (may be '-' for none)
             start_resnum - start residue sequence number
             end_resnum - end residue sequence number

        Exceptions:
           raises ValueEror if end_resnum < start_resnum

        """
        if end_resnum < start_resnum:
            raise ValueError("end residue seqnum " + str(end_resnum) + " < " +
                             "start residue seqnum " + str(start_resnum))
        self.chainid = chainid
        self.start_resnum = start_resnum
        self.end_resnum = end_resnum


    def __str__(self):
        """
        Return string representation of segment as 'chainid:start-chainid:end'
        """
        return self.chainid + ":" + str(self.start_resnum) + \
               "-" + self.chainid + ":" + str(self.end_resnum)
    

    def is_in_segment(self, res_seq_num):
        """
        Return True iff the supplied residue sequence number is in the
        interval spanned by this segment (assumed to be same chainid)

        Parameters:
           res_seq_num - PDB residue sequence number to test

        Return value:
           True if res_seq_num is in >=start_resnum and <=end_resnum else False
        """
        if res_seq_num >= self.start_resnum and \
           res_seq_num <= self.end_resnum:
            return True
        else:
            return False

    # we will define only the rich comparison operators __eq__ and __ne__
    # (not __le__, __gt__, etc.) to test for equality or non-equality
    # of segments only.

    def __eq__(self, other):
        """
        Two segments are equal if they have same chain and same start and
        end residues
        """
        if (self.chainid == other.chainid and
            self.start_resnum == other.start_resnum and
            self.end_resnum == other.end_resnum):
            return True
        else:
            return False

    def __ne__(self,other):
        """
        Two segments are '!=' (or '<>') exactly when they are not equal.
        """
        return not self.__eq__(other)



class PTDomain:
    """
    The PTDomain object represents a protein domain by a list of segments.
    Segments are contiguous subsequences of a chain, and so each is
    represented by a chain identifier, and a start and end residue sequence
    number.

    The domain consisting of domainid == None and segment_list == None
    is a special domain signifying a single-domain protein. This is used
    because we don't want to have to specify multiple segments for multiple
    chains in a single domain - a single domain protein should be treated
    just as one unit without worrying about dividing anything up.
    """
    def __init__(self, domainid, segment_list):
        """
        Create new PTDomain with supplied domain identifier and segment list.

        Parameters:
           domainid -domain identifier, string
           segment_list - list of PTSegment objects

        NOTE if both parameters are None this marks this PTDomain as the
        one used as single element in domain list to signify single-domain
        protein with no further information.
        """
        self.domainid = domainid
        self.segment_list = segment_list


    def __str__(self):
        """
        Return a representation of the domain as list of segments
        separated by ';'
        """
        if self.domainid == None and self.segment_list == None:
            return "SINGLE-DOMAIN"
        
        s = ""
        for i in range(len(self.segment_list)):
            s += str(self.segment_list[i])
            if i < len(self.segment_list) - 1:
                s += ';'
        return s

    def is_single(self):
        """
        Return True iff this is the domain with no information representing
        a single-domain protein
        """
        if self.domainid == None and self.segment_list == None:
            return True
        else:
            return False

    def is_in_domain(self, chainid, res_seq_num):
        """
        Return True iff the supplied residue specified by chainid and
        residue sequence number is in this domain.

        Parameters:
           chainid - chainid of the residue to test
           res_seq_num - PDB residue sequence number of the residue to test

        Return value:
           True if the residue is in this domain, else False.
        """
        # Note we just do a linear search, no dictionaries or anything,
        # as there is only have a maximum of maybe 5 domains, and usually
        # only 1 or 2.
        if self.is_single():
            return True # always in the special 'single domain'
        else:
            for segment in self.segment_list:
                if segment.chainid == chainid and \
                   segment.is_in_segment(res_seq_num):
                    return True
            return False

    def get_segments_for_chainid(self, chainid):
        """
        Return a list of segments of the supplied chain in this domain

        Parameters:
           chainid - id of chain to find segments of
        Return value:
           list of PTSegment objects with supplied chain id
        """
        if self.is_single():
            return []
        else:
            return [f for f in self.segment_list if f.chainid == chainid]

    def get_minmax_res_seq_in_chain(self, chainid):
        """
        Return a tuple with the lowest and highest residue sequence numbers
        in the supplied chainid in this domai