text: ###############################################################################
  #
  # tableaubuild.py - class for building protein tableaux in Python
  #
  # File:    tableaubuild.py
  # Author:  Alex Stivala
  # Created: May 2008 (moved from pytableaucreate.py)
  #
  # $Id$
  #
  ###############################################################################
  
  """
  Build a protein tableau.
  The implemntation is actually in pttableau.py which is used by ptgraph2.py
  (Pro-Origami),
  
  Also used to create SSE midpoint distance matrix.
  
  Tableaux are described by Kamat and Lesk 2007
  'Contact Patterns Between Helices and Strands of Sheet Define Protein
   Folding Patterns' Proteins 66:869-876
  and Lesk 2003 'From Electrons to Proteins and Back Again'
  Int. J. Quant. Chem. 95:678-682
  and Lesk 1995 'Systematic representation of folding patterns'
  J. Mol. Graph. 13:159-164.
  
  The implementation is based on Arun Konagurthu's TableauCreator program, see
  Konagurthu, Stuckey and Lesk 2008 'Structural search and retrieval using
  a tableau representation of protein folding patterns' Bioinformatics
  (advance access, to be published Jan 5 2008).
  
  Filenames may be either in the format above or the pdbq1lp.pdb format.
  Compressed pdb files are supported (gzip) (e.g. pdb1qlp.ent.gz).
  
  It is written in Python and depends on some Python libraries:
  
  . BioPython (including Bio.PDB)
    http://www.biopython.org
  
    Reference for Bio.PDB is:
    Hamelryck and Manderick 2003 "PDB parser and structure class implemented
    in Python" Bioinformatics 19:2308-2310
  
    which in turn depends on Numeric
    http://sourceforge.net/projects/numpy
  
  
  Developed on Linux 2.6.9 (x86_64) with Python 2.5.1
  and BioPython 1.43 with Numeric 24.2
  
  """
  
  import warnings # so we can suppress the annoying tempnam 'security' warning
  import sys,os
  import getopt
  import numpy.oldnumeric as Numeric
  from Bio.PDB import *
  
  from ptnode import *
  import pttableau
  import ptsecstruct
  from ptdomain import *
  from ptutils import cleanup_tmpdir,get_int_icode,biopdbresid_to_pdbresseq
  from domeval import *
  import getdomains
  from ptdistmatrix import compute_sse_midpoint_dist_matrix
  
  #-----------------------------------------------------------------------------
  #
  # Class definitions 
  #
  #-----------------------------------------------------------------------------
  
  
  #
  # Empty classes for exceptions
  #
  
  class NoSSE_Exception(Exception): # raised when no helices or strands found
      pass
  
  #
  # Real classes
  #
  class TableauBuild:
      """
      The protein representation consists of a sequence of structure
      (helix, strand) nodes with sequence edges in and out of them in
      sequence from N-terminus to C-terminus and adjacency edges for
      SSEs less than a threshold distance apart.
  
      Note there may be multiple such sequences (one for each
      chain).
  
      Also the nodes are all labelled with start and end residue
      sequence numbers, and node types etc. but this is not used at all
      in the here, it is only included because this code was reused from
      another program (ptraph2.py) which does require the node
      labelling.
      """
  
      #
      # member functions
      #
  
      def __init__(self, pdb_structure, pdbid,
                   include_310_helices = False, include_pi_helices = False,
                   add_loop_nodes = False):
          """
          Construct empty TableauBuild. To build the structure call
          build_graph_from_secstruct().
  
          Parameters:
              pdb_structure - parsed PDB structure from Bio.PDB
              pdbid - PDB identifier
              include_310_helices - include 3_10 helices in the graph if True
              include_pi_helices - include pi_helices in the graph if True
              add_loop_nodes - include nodes for loop regions between SSEs if True
  
          """
          self.pdb_struct = pdb_structure
          self.pdbid = pdbid
          self.chain_dict = None  # Each value of the chain_dict is a
                                  # List of nodes in order from N to C terminus
                                  # so chain_dict is { chainid : node_list }
          self.seqnum2node = {}   # dictionary of {  seqnum : PTNode }
                                  # maping int sequence numbers to PTNode objs
          self.tableau = None     # PTTableau build in build_tableau
          self.include_310_helices = include_310_helices
          self.include_pi_helices = include_pi_helices
          self.pdb_resid_dict = None # dict of { {chainid,pdb_resseq) : seqindx }
                                     # where chainid and pdb_resseq make up
                                     # the PDB residue identifier, the pdb_resseq
                                     # being string resnum+icode if any e.g.
                                     # '60' or '60A', seqindx is the indiex
                                     # into sequential list of all residues
                                     # residue_list.
          self.residue_list = None   # list of all residues (for all chains)
                                     # in sequence, built by get_residue_list()
  
  
          
      def iter_chains(self):
          """
          This generator function iterates over all chains in this PTGraph.
          A chain is just a list of nodes so it yields a node list for each
          chain.
  
          Parameters: Nonde.
          Return value: YIELDs a node list.
          Uses data members (readony):
              chain_dict - dict of {chainid:node_list}
          """
          # FIXME: can we just 'return self.chain_dict.itervalues()' here?
          for nodelist in self.chain_dict.itervalues():
              yield nodelist
              
  
      def iter_nodes(self):
          """
          This generator function iterates over all the node in this PTGraph.
  
          Parameters: None
          Return Value: YIELDs a node.
          Uses data members: (readonly):
               chain_dict - dict of {chainid_node_list}
          """
          for nodelist in self.iter_chains():
              for ptnode in nodelist:
                  yield ptnode
  
  
  
      def build_graph_from_secstruct(self, secstruct, domain, chainid=None,
                                     ignore_insertion_codes=False):
          """
          Build the list of nodes from the the supplied PTSecStruct
          object. 
  
  
          Parameters:
              secstruct - PTSecStruct (ptsecstruct.py) object to build from
              domain - PTDomain (ptdomain.py) object listing the segment(s)
                       that make up this domain (only one domain processed at a
                       time).
                       (in/out) NOTE: may be modified by having a segment
                       added if SSE is only partly in domain.
              chainid - chain identifier to build graph for only this chain,
                        or None for all chains (default)
              ignore_insertion_codes - If True, a hack to make it work with
                        PMML (only) which does not report insertion codes
                        unlike DSSP and STRIDE
  
          Uses member data (write):
              chain_dict - dict of { chainid : node_list } where node_list is
                            list of nodes in order, built in this function
              secstruct - keeps a pointer to the supplied secstruct
  
            (readonly):
              pdb_struct - The Bio.PDB parsed PDB struct (atomic co-ordinates)
                           for this protein.
              include_310_helices, include_pi_helices - if true, include
                           these kinds of helices.
  
          Raises exceptions:
             NoSSE_Exception if no helices or strands found
          
          Return value:
              None.
              
          """
  
          self.secstruct = secstruct
  
          helix_num = 1 
          strand_num = 1
  
          num_helices_in_domain = 0
          num_strands_in_domain = 0
  
          #
          # Build dictionary mapping (chainid, pdb_resid) to index in residue_list
          # for ALL residues, not just those in this domain.
          #
          self.residue_list = self.get_residue_list(self.pdb_struct,
                                                    PTDomain(None, None))
          self.pdb_resid_dict = {}
          seq_indx = 0
          while seq_indx < len(self.residue_list):
              residue = self.residue_list[seq_indx]
              self.pdb_resid_dict[( ptsecstruct.pdb_chainid_to_stride_chainid(
                                                  residue.get_full_id()[2]), 
                                    biopdbresid_to_pdbresseq(
                                                residue.get_id(),
                                                ignore_insertion_codes)
                                                                 )] = seq_indx
              seq_indx += 1
          
          # Note that now we are only adding elements in the supplied domain,
          # so the so-called 'chains' may really be segments, i.e. subsequences
          # of chains (rest of chain may be in other domain(s)
  
          self.chain_dict = {} # dict of {chainid : node_list}
  
          for (start_chainid, start_resnum, end_chainid, end_resnum, helixtype) \
                in secstruct.helix_list:
              assert(start_chainid == end_chainid) #helix must be same chain
              if chainid and chainid != start_chainid:
                  continue # chainid specified, skip ones not in that chain
              # will consider structures in domain if first residue is in domain
              if domain.is_in_domain(start_chainid,
                                     get_int_icode(start_resnum)[0]):
                  num_helices_in_domain += 1
                  if helixtype == "H":
                      idprefix = "ALPHAHELIX_"
                      htype = "ALPHA"
                      this_helix_num = helix_num
                      helix_num += 1
                  elif helixtype == "I":
                      if not self.include_pi_helices:
                          continue
                      idprefix = "PIHELIX_"
                      htype = "PI"
                      this_helix_num = helix_num
                      helix_num += 1
                  elif helixtype == "G":
                      if not self.include_310_helices:
                          continue
                      idprefix = "310HELIX_"
                      htype = "310"
                      this_helix_num = helix_num
                      helix_num += 1
                  else: # shouldn't happen
                      sys.stderr.write("ERROR: bad helix type " + helixtype+"\n")
                  ah_node = PTNodeHelix(htype,
                                        idprefix + start_chainid+"_" +\
                                        str(this_helix_num),
                                        this_helix_num,
                                        start_resnum, end_resnum, start_chainid,
                                        domain.domainid,
                                        self.residue_list, self.pdb_resid_dict)
                  if not self.chain_dict.has_key(start_chainid):
                      self.chain_dict[start_chainid] = []
                  self.chain_dict[start_chainid].append(ah_node)
  
                  # we must already have handled the case of SSEs that cross
                  # domain boundaries (by moving whole SSE to one of the domains)
                  assert( domain.is_in_domain(end_chainid,  get_int_icode(end_resnum)[0]) )
  
          for (start_chainid, start_resnum, end_chainid, end_resnum) \
                  in secstruct.strand_list:
              assert(start_chainid == end_chainid) # must be in same chain
              if chainid and chainid != start_chainid:
                  continue # chainid specified, skip ones not in that chain
              if domain.is_in_domain(start_chainid,
                                     get_int_icode(start_resnum)[0]):
                  num_strands_in_domain += 1
                  bs_node = PTNodeStrand("STRAND_"+start_chainid +"_"+\
                                         str(strand_num),
                                         strand_num,
                                         start_resnum, end_resnum, start_chainid,
                                         domain.domainid,
                                         self.residue_list,
                                         self.pdb_resid_dict)
                  strand_num += 1
                  if not self.chain_dict.has_key(start_chainid):
                      self.chain_dict[start_chainid] = []
  
                  # we must already have handled the case of SSEs that cross
                  # domain boundaries (by moving whole SSE to one of the domains)
                  assert( domain.is_in_domain(end_chainid,  get_int_icode(end_resnum)[0]) )
                  self.chain_dict[start_chainid].append(bs_node)
  
  
          # raise an exception if there are no SSEs at all in this domain
          if num_helices_in_domain == 0 and num_strands_in_domain == 0:
              raise NoSSE_Exception
  
          delete_chainid_list = [] # list of chainids to delete from chain_dict
          for (chainid, nodelist) in self.chain_dict.iteritems():
              # sort in order of start residue id ascending (all must be disjoint)
              nodelist.sort()
  
              if len(nodelist) < 1:
                  # There are no SSEs in this chain, get rid of it.
                  sys.stderr.write('WARNING: no SSEs in chain ' + chainid +
                                   '; chain ignored\n')
                  delete_chainid_list.append(chainid) # don't delete while in loop
                  continue
              else:
                  # Check for chain with only SSEs that will not be drawn
                  # (i.e. pi or 310 helices), and delete those too
                  found_useful_node = False
                  for ptnode in nodelist:
                      if isinstance(ptnode, PTNodeStrand):
                          found_useful_node = True
                          break
                      elif isinstance(ptnode, PTNodeHelix):
                          if ptnode.get_type() == "ALPHA":
                              found_useful_node = True
                              break
                          elif ((ptnode.get_type() == "310" and
                                   self.include_310_helices) or
                                  (ptnode.get_type() == "PI" and
                                   self.include_pi_helices)):
                              found_useful_node = True
                              break
                  if not found_useful_node:
                      sys.stderr.write('WARNING: only pi or 310 helices in chain '
                                       + chainid +
                                       '; chain ignored\n')
                      delete_chainid_list.append(chainid)
                      continue
                  
  
          # delete chains from chain_dict that were marked earlier for deletion
          for chainid in delete_chainid_list:
              self.chain_dict.pop(chainid)
  
          # -------------------------------------------------------------------
  
          # This is needed only for labelling sheets for HH and KK codes
          # (see dfs_strands() etc. below)
          
          # add edges for hydrogen bonds
          # uses secstruct and chainid member data
          # these are used for determining which side bridge partners are
          # on (and also for drawing a hydrogen bond graph if requested)
          self.add_hbond_edges_from_secstruct()
          
          # add edges for bridge partners
          # uses secstruct and chainid member data
          self.add_bridge_edges_from_secstruct()
  
          #---------------------------------------------------------------------
  
  
          # for sequential numbering, we'll build this dictionary mapping
          # sequential number (note NOT restarting for each chain)
          # to PTNode
          # so that sequential numbers as used in ptgraph2 -b sequential
          # option.
          # this is a dictionary of { seqnum : PTNode }
          self.seqnum2node = {}
          for (seqnum, node) in \
              enumerate([node for node in self.iter_nodes() if \
                         not ( (isinstance(node, PTNodeTerminus)) or
                                (isinstance(node, PTNodeHelix) and
                                 ( (node.get_type() == "310" and
                                    not self.include_310_helices) or
                                   (node.get_type() == "PI" and
                                    not self.include_pi_helices) ) ) ) ]):
                  self.seqnum2node[seqnum+1] = node # start at 1 not 0
  
      # ------------------------------------------------------------------------
  
      def get_residue_list(self, pdb_struct, domain, getchainid = None):
          """
          Return list of Bio.PDB Residue objects in this domain, and optionally
          in the specified chain.,
  
          Parameters:
               pdb_struct - Bio.PDB parsed PDB struct for the protein
               domain -  PTDomain (ptdomain.py) object listing the segment(s)
                           that make up this domain (only one domain processed at a
                           time).
               getchainid - chain identifier to get residues in (default None -
                         all chains).
  
          Return value:
               list of Bio.PDB Residue objects in the domain (and optionally chain).
          Raises exceptions:
             NoSSE_Exception for empty structure (happens eg on d1oayi_.ent)
  
          """
          residue_list = []
          try:
              pdb_model = self.pdb_struct[0] # TODO always using model 0 for now
          except KeyError:
              raise NoSSE_Exception
          
          for chain in pdb_model:
              chainid = ptsecstruct.pdb_chainid_to_stride_chainid(chain.get_id())
              if getchainid and getchainid != chainid:
                  continue # this is not the chain we want
  
              # Build a list of Bio.PDB Residue objects that are in this
              # domain.
              # id of a residue in Bio.PDB is tuple (hetatm, resseqnum, icode)
              # so we choose those where residue PDB number
              # (in the current chain) is in the domain.
              # TODO: maybe should use polypeptide builder for this instead
              # (and indeed should probably use it right from the beginning) -
              residue_list += [ residue for residue in chain.get_unpacked_list()
                                if is_aa(residue) and
                                domain.is_in_domain(chainid, residue.get_id()[1])
                              ]
              if getchainid:
                  break # if getchainid specified, we now have it so can quit
          return residue_list
  
  
      # iter_strands(), dfs_strands(),
      # find_connected_components() and label_sheets() are needed to assign
      # strands to sheets in order for HH and KK codes to be used for strands
      # only when they are in the same sheet.
      # TODO this code is cut&pasted from ptgraph2.py, probably should
      # have a base class that this at PTGraph2 both inherit from or somethibng
      # rather than all this duplication.
  
      def add_hbond_edges_from_secstruct(self):
          """
          Add edges between structural elements for hydrogen bonds between
          those nodes. Called by build_graph_from_secstruct().
  
          NB: adds bonds between STRANDs only, not between HELIXes (helices).
          
          Parameters: None.
          Return value: None.
          Uses data members:
             readonly:
                secstruct - PTSecStruct object to get hbonds from
                chainid - chainid of chain in PTSecStruct to use
             read/write:
                chain_dict - dict by chainid of
                             list of nodes (changes node data, not list as such)
  
          Precondition: each nodelist in chain_dict
                        is sorted (by start res seq ascending);
                        this is done by build_graph_from_secstruct()
                        before calling.
  
          """
          hbond_list = self.secstruct.hbond_list
          # TODO: do this more efficiently using presorting (ie how it used to
          # be done when only one chain)
          for (chainid1, resnum1, chainid2, resnum2, dist) in hbond_list:
              for ptnode in self.iter_strands():
                  if chainid1 == ptnode.get_chainid() and \
                     resnum1 >= ptnode.get_start_res_seq() and \
                     resnum1 <= ptnode.get_end_res_seq():
                      dest_node = self.find_node_containing_seqnum(resnum2,
                                                                   chainid2)
                      if dest_node != None and \
                             isinstance(dest_node, PTNodeStrand): # only STRANDs
                          ptnode.add_hbond(dest_node, resnum1, resnum2, dist)
  
      def add_bridge_edges_from_secstruct(self):
          """
          Add edges between strand nodes representing beta brdiges between
          those nodes (add just one edge between any two strands).
          Called by build_graph_from_secstruct().
  
          NB: adds bonds between STRANDs only, not between HELIXes (helices).
  
          Parameters: None.
          Return value: None.
          Uses data members:
             readonly:
                secstruct - PTSecStruct object to get hbonds from
                chainid - chainid of chain in PTSecStruct to use
             read/write:
                chain_dict - dict by chainid of
                              list of nodes (changes node data, not list as such)
  
          """
  
          bridge_list =  self.secstruct.bridgeres_list
          #         (chainid1, resnum1, chainid2, resnum2, bdir)
   
          # TODO: do this more efficiently using presorting (ie how it used to
          # be done when only one chain)
  
          for ptnode in self.iter_strands():
              for (chainid1, resnum1, chainid2, resnum2, bdir) in bridge_list:
                  if chainid1 == ptnode.get_chainid() and \
                     resnum1 >= ptnode.get_start_res_seq() and \
                     resnum1 <= ptnode.get_end_res_seq():
                      try:
                          dest_node = self.find_node_containing_seqnum(resnum2,
                                                                       chainid2)
                      except KeyError:
                          dest_node = None
                          sys.stderr.write('WARNING: chain ' + chainid2 + \
                                           ' involved in beta bridge not found.'+\
                                           '\n  Probably due to domain parsing' +\
                                           ' breaking a beta sheet.\n')
                      if dest_node != None and \
                             isinstance(dest_node, PTNodeStrand): # only STRANDs
                          if ptnode == dest_node:
                              sys.stderr.write('WARNING: ignoring self-bridge ' +
                                               ptnode.nodeid + '\n')
                          else:
                              ptnode.add_bridge(dest_node, bdir)
  
      
      def iter_strands(self):
          """
          This generator function iterates over all strands in this PTGraph
          object. I.e. it yields a strand for each strand in the 
          node lists.
  
          Parameters: None.
          Return value: YIELDs a strand.
          Uses data members (readonly):
             self.chain_dict - dict of { chainid : list of nodes }
          """
          for nodelist in self.iter_chains():
              for ptnode in nodelist:
                  if isinstance(ptnode, PTNodeStrand):
                      yield ptnode
  
      def find_node_containing_seqnum(self, res_seqnum, chainid):
          """
          Find and return node in node list for chain chainid
          containing supplied PDB residue
          sequence number.
  
          Parameters:
             res_seqnum - PDB residue sequence number to find node for
             chainid - chain identifier to find node in
  
          Return value:
             PTNode pointer of PTNode containing the supplied residue seq num
             in supplied chainid
             or None if the residue is not in a structural element PTNode
  
          Uses data members (readonly):
             chain_dict - chainid dict of list of PTNodes
          """
          # TODO: since node_list is sorted should use binary search here
          #       (maybe try the Python bisect module)
          if not self.chain_dict.has_key(chainid):
              return None # no such chain, can happen due to domain parsing
          for ptnode in self.chain_dict[chainid]:
              if ptnode.is_in_interval(res_seqnum):
                  return ptnode
          return None
  
      def dfs_strands(self, start_strand, visited, dfs_list, from_node,
                      back_edge_list,
                      sheet_id=None):
          """
          Make a depth-first search traversal of STRAND nodes
          using bridge (not sequence)
          edges starting at the specfied strand.
  
          Parameters:
             start_strand - STRAND node to start at
             visited - (in/out) dictionary of {ptnode:True} visited nodes
             dfs_list - (in/out) list of ptnodes visited in dfs order
             from_node - node from which we are being (recursively) called
             back_edge_list - list of (node, node) tuples representing an
                               edge between the two nodes, which is a back
                               edge, i.e. from a node to an ancestor of that
                                node in the spanning tree. The back edge
                                means there is a cycle of which the back
                                edge forms a part.
             sheet_id - identifier of this sheet (connected component) to mark
                        each strand in it with, or None to not mark at all
                        (default).
                        
  
          Recursive function. call initially as
              dfslist = []
              back_edge_list = []
              dfs_strands(startnode, {}, dfslist, None, back_edge_list)
  
          Return value:
              None. (Output is dfs_list, back_edge_list parameters)
          
          Uses members (readonly):
             chain_dict - dict by chainid of list of PTNodes
             
          """
          visited[start_strand] = True
          if sheet_id != None:
              start_strand.set_sheet_id(sheet_id)
              #print 'xxx',str(start_strand),sheet_id
          dfs_list.append(start_strand)
          for (node, bdir_unused, side_unused) in start_strand.get_bridge_list():
              if node not in visited:
                  self.dfs_strands(node, visited, dfs_list, start_strand,
                                   back_edge_list, sheet_id)
              elif node != from_node: #not parent of start_strand in spanning tree
                  # don't add duplicate back edges
                  # ((node1,node2) is same as (node2,node1))
                  duplicate = False
                  for (a,b) in back_edge_list:
                      if ((start_strand == a and node == b) or
                          (node == a and start_strand == b)):
                          duplicate = True
                          break
                  if not duplicate:
                      if verbose:
                          sys.stderr.write('dfs_strands back edge from ' +
                                           str(start_strand) + ' to ' +
                                           str(node) +
                                           '\n')
                      back_edge_list.append((start_strand, node))
  
  
  
      def find_connected_components(self):
          """
          Find the connected components (considering only STRAND nodes
          and bridge [not sequence] edges in the graph).
  
          This is done by a DFS traversal at every node in the graph
          (skipping already visited ones), giving us the partition of
          the graph into connected components.
          
          Parameters: None
  
          Uses member data: 
              chain_dict - dict by chainid of list
                          of PTNodes in the graph (modifies PTNodes not list)
  
              (WRITE):
  
               sheet_dict - 
                dictionary of { sheet_id : ptnode_list } where sheet_id is 'A',
                'B', etc. and ptnode_list is a list of PTNodeStrand instances
                 in that connected component (sheet).
  
                 self.sheet_backedges_dict  -
                   dict of {sheet_id : ((node1,node2))}
                   listing 'back edges' i.e. edges
                   to an ancestor in DFS spanning tree
                   in the connected component (sheet).
                   note (node1,node2) and (node2,node1)
                   are the same (undirected graph) and
                   only one of the two is present in the
  
  
          Labels each strand node with the sheet id it belongs to as it goes.
          """
  
          sheet_id = 'A' # sheet id is single alpha char A, B, etc.
                         # (will be a problem for more than 26 sheets... eg
                         # this actually happens on 2J28), wrap to lowercase
                         
          visited = {}   # dictionary of {ptnode : True} visited nodes
          back_edge_list = []  # list of (ptnode, ptnode) tuples for back edges
          self.sheet_dict = {} # dictionary of {sheet_id : nodelist}
          self.sheet_backedges_dict = {} # dict of {sheet_id : ((node1,node2))}
                                         # listing 'back edges' i.e. edges
                                         # to an ancestor in DFS spanning tree
                                         # in the connected component (sheet).
                                         # note (node1,node2) and (node2,node1)
                                         # are the same (undirected graph) and
                                         # only one of the two is present in the
                                         # list.
          for node in self.iter_strands():
              if node not in visited:
                  connected_node_list = []
                  back_edge_list = []
                  self.dfs_strands(node, visited, connected_node_list, None,
                                   back_edge_list,
                                   sheet_id)
                  self.sheet_dict[sheet_id] = list(connected_node_list)
                  self.sheet_backedges_dict[sheet_id] = list(back_edge_list)
                  sheet_id = chr(ord(sheet_id)+1)
                  if sheet_id == '[':
                      sheet_id = 'a' # if go past Z, wrap to lowercase
  
      
      def label_sheets(self):
          """
          Label strands with sheet id to which each belongs by finding
          connected components; strands in a connected componenent of
          the graph (considering nonly STRAND nodes and bridge edges)
          form a sheet.
  
          Parameters: None
  
          Uses member data:
             node_list - list of nodes. Modifies nodes by labelling them.
  
          Return value:
             Returns the sheet dictionary (dictionary of
             { sheet_id : ptnode_list }) from find_connected_components.
          """
          # ACtually don't do anything except call find_connected_components()
          # which does the labeling itself (more efficient since it knows
          # as each one is added which sheet it is added to)
          return self.find_connected_components()
              
  
      # -------------------------------------------------------------------------
      
      def build_tableau(self, pdbid, domain, ptnode_list = None,
                        use_hk = True):
          """
          Build the tableau data member (see PTTableau in pttableau.py)
          by calling function in  pttableau.py.
  
          Parameters:
             pdbid - PDB identifier of the strucutre
             domain - The PTDomain object for our current domain
             ptnode_list - list of PTNodes (in sequence order, but not
                           necessarily continguous) to build the tableau for,
                           or None to use all nodes in domain.
                           Default None.
          use_hk - If True, use the HH and KK codes for respectively
                   antiparallel and parallel strands. Default True.
  
          Return value: None
          Uses data members (WRITE):
             tableau - created by this function
             (readonly):
             chain_dict - dict { chainid : ptnode_list } of nodes in chains
             pdb_structure - Bio.PDB parsed PDB structure
          """
          if ptnode_list == None:
              # Build list of all helix and strand PTNodes
              ptnode_list = []
              for nodelist in self.iter_chains():
                  for node in nodelist: # these nodes are only those in our domain
                      if (not isinstance(node, PTNodeTerminus)): # not terminii
                          ptnode_list.append(node)
  
          self.tableau = pttableau.compute_tableau(ptnode_list, self.pdb_struct,
                                                   use_hk)
  
  
  
      def build_omega_matrix(self, pdbid, domain, ptnode_list = None):
          """
          Return the relative angles matrix by calling function in pttableau.py
  
          Parameters:
             pdbid - PDB identifier of the strucutre
             domain - The PTDomain object for our current domain
             ptnode_list - list of PTNodes (in sequence order, but not
                           necessarily continguous) to build the tableau for,
                           or None to use all nodes in domain.
                           Default None.
          Return value: Numeric.array Omega matrix.
          Uses data members:
             (readonly):
             chain_dict - dict { chainid : ptnode_list } of nodes in chains
             pdb_structure - Bio.PDB parsed PDB structure
          """
          if ptnode_list == None:
              # Build list of all helix and strand PTNodes
              ptnode_list = []
              for nodelist in self.iter_chains():
                  for node in nodelist: # these nodes are only those in our domain
                      if (not isinstance(node, PTNodeTerminus)): # not terminii
                          ptnode_list.append(node)
  
          return pttableau.compute_omega_matrix(ptnode_list, self.pdb_struct)
  
  
      def build_sse_dist_matrix(self, pdbid, domain, ptnode_list = None):
          """
          Return SSE axis midpoint distance matrix by calling function
          in ptdistmatrix.py
  
          Parameters:
             pdbid - PDB identifier of the strucutre
             domain - The PTDomain object for our current domain
             ptnode_list - list of PTNodes (in sequence order, but not
                           necessarily continguous) to build the matrix for,
                           or None to use all nodes in domain.
                           Default None.
          Return value: Numeric.array SSE midpoint distance matrix.
          Uses data members:
             (readonly):
             chain_dict - dict { chainid : ptnode_list } of nodes in chains
             pdb_structure - Bio.PDB parsed PDB structure
          """
          if ptnode_list == None:
              # Build list of all helix and strand PTNodes
              ptnode_list = []
              for nodelist in self.iter_chains():
                  for node in nodelist: # these nodes are only those in our domain
                      if (not isinstance(node, PTNodeTerminus)): # not terminii
                          ptnode_list.append(node)
  
          return compute_sse_midpoint_dist_matrix(ptnode_list, self.pdb_struct)
  
  
  
              
  #-----------------------------------------------------------------------------
  #
  # Function definitions
  #
  #-----------------------------------------------------------------------------
  
  
  
  def make_tableaux(pdb_filename,
                  pdb_struct,
                  secstruct_program,
                  domain_program,
                  include_310_helices = False,
                  include_pi_helices = False,
                  use_numeric = False,
                  sse_id_list = None,
                  use_hk = False,
                  min_sse_len = None,
                  build_dist_matrix = False,
                  chainid = None,
                  domainid = None):
      """
      For the supplied filemame, read PDB format data from that file
      and create tableaux or SSE distance matrix for that structre.
      This function is called by get_tableaux() (below), which handles
      compressed files etc.
  
      Note: for multi-domains, will be multiple output tableaux, one for
      each domain.
  
      Paramteters:
         pdb_filename - filename of PDB file to read
         pdb_struct - Bio.PDB parsed PDB structure
         secstruct_program - secondary structure definition program
                         ('stride' or 'dssp' or 'pdb' or 'pmml') to use.
         domain_progam - domain decompositino method ('ddomain','cath', etc.)
         include_310_helices - if True, include 3_10 helices in the graph
         include_pi_helices - if True, include pi helices in the graph
         use_numeric - If True, use numeric matrix rather than tableau
         sse_id_list - list of ints representing SSE sequential id numbers
                       to include in tableau. Default None.
                       When None, all SSEs are included.
         use_hk - If True, use HH and KK codes for strands.
         min_sse_len - if not None, the minimum length of SSE to include
                      in tableau.
         build_dist_matrix - If True, build SSE midpoint distance matrix
                     instead of tableau.
         chainid - If not None, only build tableau for that chain id.
         domainid - If note None, only build tableau for that domain id.
  
      Return value: tuple (tableaux_list, sse_string_list)
                     where tableaux_list is
                     list of tableaux (only one in list unless domain decomp
                     is used and finds multiple domains);
                     or list of omega matrices (Numeric.array) if use_numeric
                     is True
                     or list of SSE axis midpiont distance matrices
                     (Numeric.array) if build_dist_matrix is True
                     and 
                     sse_string_list is SSE string description e.g. 'EEHHE' etc.
      """
      (pdbid,suffix) = os.path.splitext(os.path.basename(pdb_filename))
      pdbid = pdbid.upper()
      if len(pdbid) >= 6 and pdbid[:3] == "PDB":
          pdbid = pdbid[3:7]
  
      if secstruct_program == "pdb":
          secstruct = ptsecstruct.read_secstruct_from_pdb_file(pdb_filename)
          if secstruct != None:
              secstruct.pdb_header = pdb_struct.header['head']
          else:
              secstruct_program = "dssp"
              sys.stderr.write('WARNING: error with HELIX or SHEET cards in PDB'
                               ': ' + secstruct_program +
                               ' will be used instead\n')
      else:
          secstruct = None
  
      if secstruct == None:
          # read secondary structure information from STRIDE or DSSP
          if secstruct_program == "stride":
              secstruct = ptsecstruct.read_secstruct_from_stride(pdb_filename)
          elif secstruct_program == "dssp":
              secstruct = ptsecstruct.read_secstruct_from_dssp(pdb_filename)
          elif secstruct_program == "pmml":
              secstruct = ptsecstruct.read_secstruct_from_pmml(pdb_filename)
          else:
              assert(False)
  
      if domain_program != None:
          domain_list = getdomains.get_domains(domain_program,
                                               pdbid, pdb_filename, pdb_struct)
      else:
          domain_list = [PTDomain(None, None)] # one-domain protein, no further info
  
  
      # for SSEs that cross domain boundaries, move whole SSE to one of the domains
      fixup_crossdomain_sses(secstruct, domain_list)
  
      tableaux_list = [] # NB may be list of PTTableau or list of Numeric.array
      sse_str_list = []
      for domain in domain_list:
          if domainid and domain.domainid != domainid:
              if verbose:
                  sys.stderr.write("skipped domainid " + domainid + "\n")
              continue
          
          ptg = TableauBuild(pdb_struct, pdbid,
                             include_310_helices, include_pi_helices)
          # build tableaubuild object from secondary structure
          if secstruct_program == "pmml":
              ignore_insertion_codes = True #  PMML does not report them
          else:
              ignore_insertion_codes = False
          try:
              ptg.build_graph_from_secstruct(secstruct, domain, chainid,
                                             ignore_insertion_codes)
          except NoSSE_Exception:
              if chainid:
                  sys.stderr.write('WARNING: No helices or strands found in ' +
                                   pdbid +
                                   ' chain ' + chainid +
                                   ': skipping\n')
              else:
                  sys.stderr.write('WARNING: No helices or strands found in ' +
                                   pdbid +
                                   ': skipping\n')
              continue
          except ValueError:
              # PTNode.__init__ raises ValueEror if end_res_seq < start_res_seq
              # This happens on e.g. d1dkia_.ent when using PMML (only)
              # as PMML does not handle insertion codes so end up with
              # residue numbers all wrong. Nothign we can do with this much of
              # a mess, jsut give up.
              sys.stderr.write('ERROR: out of order residue numbers in ' +
                               pdbid + ', probably due to insertino codes ' +
                               'not workin in pmml: skipping this structure\n')
              continue
                  
  
          if use_hk: # only need to know sheets if using HH and KK codes
              if secstruct_program == "pmml":
                  sys.stderr.write('WARNING: PMML has no hbond or bridge information, cannot label sheets for HH and KK codes\n')
              ptg.label_sheets()
  
          if verbose:
              for nodelist in ptg.iter_chains():
                  for node in nodelist:
                      sys.stderr.write(str(node) + '\n')
  
          # if list of int SSE sequential ids supplied, convert to list of
          # PTNode objects
          if sse_id_list:
              try:
                  ptnode_list = [ptg.seqnum2node[sse_id] for sse_id in sse_id_list]
              except KeyError,k:
                  sys.stderr.write("SSE sequential id " + str(k)
                                   + " does not exist\n")
                  sys.exit(1)
          else:
              ptnode_list = None
  
          if not ptnode_list:
              # Build list of all helix and strand PTNodes with len >= min_sse_len
              ptnode_list = []
              for nodelist in ptg.iter_chains():
                  for node in nodelist: # these nodes are only those in our domain
                      if (not isinstance(node, PTNodeTerminus)): # not terminii
                          ptnode_list.append(node)
  
          if min_sse_len:
              ptnode_list = [node for node in ptnode_list
                             if node.get_span() >= min_sse_len]
              
          if build_dist_matrix:
              dist_matrix = ptg.build_sse_dist_matrix(pdbid, domain, ptnode_list)
              tableaux_list.append(dist_matrix)
          elif use_numeric:
              Omega = ptg.build_omega_matrix(pdbid, domain, ptnode_list)
              tableaux_list.append(Omega)
          else:
              ptg.build_tableau(pdbid, domain, ptnode_list, use_hk)
              tableaux_list.append(ptg.tableau)
  
          sse_str = ""
          for node in ptnode_list:
              if isinstance(node, PTNodeStrand):
                  sse_str += 'E'
              elif isinstance(node, PTNodeHelix):
                  sse_str += 'H'
              else:
                  raise ValueError('bad node type ' + str(node))
          sse_str_list.append(sse_str)
          
          
      return (tableaux_list, sse_str_list)
  
  
  
  def get_tableaux(pdb_filename,
                   secstruct_program = 'dssp',
                   domain_program = 'none',
                   include_310_helices = True,
                   include_pi_helices = True,
                   sse_id_list = None,
                   min_sse_len = None,
                   use_numeric = False, 
                   use_hk = False,
                   build_dist_matrix = False):
      
      """
      Get a tableau for a single PDB or ASTRAL pdb-style file
      (compressed files e.g. pdb1qlp.ent.gz) or uncompressed
      or the ASTRAL pdb-style hierarchy
      (uncompressed files e.g. d1qlpa_.ent).
  
      Parameters:
         pdb_filename - filename of PDB or ASTRAL pdb-style file, as above.
         secstruct_program - secondary structure definition program
                         ('stride' or 'dssp' or 'pdb') to use.
         domain_progam - domain decompositino method ('ddomain','cath', etc.)
         include_310_helices - if True, include 3_10 helices in the graph
         include_pi_helices - if True, include pi helices in the graph
         sse_id_list - list of ints representing SSE sequential id numbers
                       to include in tableau. Default None.
                       When None, all SSEs are included.
         min_sse_len - min number of residues in SSE to be ncluded.
                        Default None (no min length).
         use_numeric - if True build Numeric.array Omega matrix not PTTableau
         use_hk - If True build tableaux with HH and KK codes for strands in
                  same sheet. default False.
         build_dist_matrix - If True, build SSE midpoint distance matrices
                     instead of tableaux.
  
  
      Return value:
          tuple (pdbid, tableaux_list, sse_string_list)
           from the pdb file, only one in lists unless
           domain decomposition is used and finds multidomains in input.
           tableaux_list is list of tableaux or omega matrices
           sse_string_list is SSE string description e.g. 'EEHHE' etc.
      """
      tableaux_list = []
      # check for compressed files. We only support gzip (.gz)
      # Note we are not using the zlib or GzipFile python modules
      # since we are calling to external programs which require the
      # file uncompressed themsevles anyway so we'll just run gzip
      # to uncompress the file to a temporary directory.
      pdb_file_basename = os.path.basename(pdb_filename)
      (name,extension) = os.path.splitext(pdb_file_basename)
      if extension == '.gz':
          TMPDIR = os.tempnam(None, "ptgz")
          os.mkdir(TMPDIR)
          tmp_pdbfilename = os.path.join(TMPDIR, name)
          os.system("gzip " + pdb_filename + " -d -c > " + tmp_pdbfilename)
          our_pdb_filename = tmp_pdbfilename
          used_tmp_file = True
      else:
          our_pdb_filename = pdb_filename
          used_tmp_file = False
  
      try:
          pdbid = name
          if len(pdbid) >= 6 and pdbid[:3].upper() == "PDB":
              pdbid = pdbid[3:7].upper()
          # parse PDB file
          pdb_parser = PDBParser()
          pdb_struct = pdb_parser.get_structure(pdbid, our_pdb_filename)
          # create the Tableaux and output them
          (tableaux_list, sse_string_list) = make_tableaux(our_pdb_filename,
                      pdb_struct,
                      secstruct_program,
                      domain_program,
                      include_310_helices,
                      include_pi_helices,
                      use_numeric,
                      sse_id_list,
                      use_hk,
                      min_sse_len,
                      build_dist_matrix)
                                     
      finally:
          if used_tmp_file:
              cleanup_tmpdir(TMPDIR)
      return (pdbid, tableaux_list, sse_string_list)
  
  
  