OpenStructure
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
align_3dcomb.py
Go to the documentation of this file.
1 """
2 3DCOMB module
3 
4 Author: Niklaus Johner
5 
6 This module is for structural alignment from OpenStructure using the external program 3DCOMB.
7 
8 How To Use This Module:
9  1. Import it (e.g. as "from ost.bindings import align_3dcomb")
10  2. Use it (e.g. as "alignment,transformation_list = align_3dcomb.AlignStructures(view_list)")
11 
12 Requirement:
13  - 3DCOMB installed
14 """
15 
16 from ost.bindings import utils
17 import subprocess,os
18 from ost import settings
19 from ost import io
20 import ost
21 import ost.geom
22 
23 def _GetExecutable(comb_exe, comb_env):
24  """
25  Function to check if 3DCOMB executable is present
26 
27  :param comb_exe: Explicit path to 3DCOMB executable
28  :param msms_env: Environment variable pointing to 3DCOMB executable
29  :returns: Path to the executable
30  :raises: :class:`~ost.FileNotFound` if executable is not found
31  """
32  return settings.Locate(['3DCOMB_linux','3DCOMB_win.exe'], explicit_file_name=comb_exe,
33  env_name=comb_env)
34 
35 def _SetupFiles(structure_list):
36  """
37  Setup files for MSMS calculation in temporary directory
38 
39  :param structure_list: A list of EntityView that will be aligned.\
40  each EntityView should contain a single chain and each residue needs to have a CA atom.
41  :returns: calss:settings.TempDir
42  :raises: class:`RuntimeError` if on of the Views is not valid
43  """
44 
45  #write out temporary pdb files
46  if not all([ev.IsValid() for ev in structure_list]):
47  raise RuntimeError, "Invalid EntityView in structure_list"
48  tpd=utils.TempDirWithFiles(structure_list)
49 
50  #Write out the file containing the list of all structures
51  outfile=open(os.path.join(tpd.dirname,'list'),'w')
52  outfile.write('\n'.join(tpd.files))
53  outfile.close()
54  return tpd
55 
56 
57 def _Run3DCOMB(command,tpd):
58  """
59  Run the 3DCOMB alignment command
60 
61  This functions starts the external 3DCOMB executable and returns the stdout of
62  3DCOMB.
63 
64  :param command: Command to execute
65  :returns: stdout of 3DCOMB
66  :raises: :class:`CalledProcessError` for non-zero return value
67  """
68  outname=os.path.join(tpd.dirname,'align.out')
69  outfile=open(outname,'w')
70  returncode=subprocess.call(command, shell=True, stdout=outfile)
71  outfile.close()
72  #check for successful completion of 3DCOMB
73  if returncode!=0:
74  print "WARNING: 3DCOMB error\n"
75  raise subprocess.CalledProcessError
76  return returncode
77 
78 def _ParseOutput(tpd):
79  #Read Alignment
80  ali=io.LoadAlignment(os.path.join(tpd.dirname,'list.ali'),'fasta')
81  for i in range(ali.GetCount()):
82  ali.SetSequenceName(i,'struc{0}'.format(i))
83  #Read Transformations
84  f=open(os.path.join(tpd.dirname,'list.rmt'),'r')
85  Tl=[]
86  for l in f:
87  if l.startswith('>'):
88  fl=ost.FloatList()
89  for i in range(3):
90  l=f.next()
91  sl=l.split(',')
92  fl.extend([float(el) for el in sl[0].split()+[sl[1]]])
93  fl.extend([0,0,0,1])
94  T=ost.geom.Transform()
95  M=ost.geom.Mat4(*fl)
96  T.SetMatrix(M)
97  Tl.append(T)
98  #Read standard output
99  outfile=open(os.path.join(tpd.dirname,'align.out'),'r')
100  results={}
101  for line in outfile:
102  if line.startswith('Objective function value is'):results['objective_function']=float(line.split()[-1])
103  if line.startswith('CORE_LEN'):
104  l=line.split(',')
105  for el in l:
106  s=el.split('=')
107  results[s[0]]=float(s[1])
108  return ali,Tl,results
109 
110 
111 def AlignStructures(structure_list,apply_transform=False,name_list=None,comb_exe=None,comb_env=None):
112  comb_executable=_GetExecutable(comb_exe, comb_env)
113  tpd=_SetupFiles(structure_list)
114  command=' '.join([comb_executable,os.path.join(tpd.dirname,'list')])
115  returncode=_Run3DCOMB(command,tpd)
116  try:ali,Tl,results=_ParseOutput(tpd)
117  except:
118  print 'could not parse output'
119  raise RuntimeError
120  if apply_transform:
121  for T,ev in zip(Tl,structure_list):
122  ev.handle.FixTransform()
123  ev.handle.SetTransform(T)
124  tpd.Cleanup()
125  return ali,Tl,results
126 
127 
128 
129