00001 """
00002 Interface to LGA
00003
00004 Author: Marco Biasini
00005 """
00006 import tempfile
00007 import os
00008 from ost import io
00009 import subprocess
00010 from ost import geom
00011 from ost import settings
00012
00013 def _FindLGABinary(lga_bin):
00014 return settings.Locate('lga', explicit_file_name=lga_bin,
00015 env_name='LGA_BINARY')
00016
00017 def _PrepareInput(pdb1, pdb2, output_pdb):
00018 """
00019 Deal with the brain-dead input requirements of LGA.
00020 """
00021 mol1=os.path.basename(os.path.splitext(pdb1)[0])
00022 mol2=os.path.basename(os.path.splitext(pdb2)[0])
00023 os.system('echo MOLECULE %s > %s' % (mol1, output_pdb))
00024 os.system('cat %s >> %s' % (pdb1, output_pdb))
00025 os.system('echo "\n" >> %s' % (output_pdb))
00026 os.system('echo MOLECULE %s >> %s' % (mol2, output_pdb))
00027 os.system('cat %s >> %s' % (pdb2, output_pdb))
00028
00029 class GDTResult:
00030 def __init__(self, rotation, shift, gdt_ts, gdt_ha):
00031 self.rotation=rotation
00032 self.shift=shift
00033 self.gdt_ts=gdt_ts
00034 self.gdt_ha=gdt_ha
00035 def GetTransform(self):
00036 transform=geom.Mat4()
00037 transform.PasteTranslation(self.shift)
00038 return transform*geom.Mat4(self.rotation)
00039
00040 def _ParseRotationAndShift(lines):
00041 t=[l.split() for l in lines]
00042 rot=geom.Mat3()
00043 shift=geom.Vec3()
00044 for i, x in enumerate(t):
00045 rot[(i, 0)]=+float(x[2])
00046 rot[(i, 1)]=+float(x[6])
00047 rot[(i, 2)]=+float(x[10])
00048 shift[i]=float(x[14])
00049 return rot, shift
00050
00051 def _ParseGDTSection(section, residue_count):
00052 cutoffs=[float(e) for e in section[0].split()[2:]]
00053 num_ca=[int(e) for e in section[1].split()[2:]]
00054 gdt_ts=[float(e) for e in section[2].split()[2:]]
00055 scores=dict(zip(cutoffs, gdt_ts))
00056 numbers=dict(zip(cutoffs, num_ca))
00057 factor=(1.0/(4*residue_count))*100
00058 ts_cutoffs=(1.0, 2.0, 4.0, 8.0)
00059 ha_cutoffs=(0.5, 1.0, 2.0, 4.0)
00060 gdt_ts=(sum([numbers[c] for c in ts_cutoffs]))*factor
00061 gdt_ha=(sum([numbers[c] for c in ha_cutoffs]))*factor
00062 return gdt_ts, gdt_ha
00063
00064 def _ParseLGAOutput(output, residue_count):
00065 result=GDTResult(geom.Mat3(), geom.Vec3(), 0.0, 0.0)
00066 found_gdt_section=False
00067 found_transform_section=False
00068 for index, line in enumerate(output):
00069 if line.startswith('GLOBAL_DISTANCE_TEST'):
00070 next_lines=output[index+1:index+5]
00071 result.gdt_ts, result.gdt_ha=_ParseGDTSection(next_lines, residue_count)
00072 found_gdt_section=True
00073 if line.startswith('Unitary ROTATION matrix'):
00074 next_lines=output[index+1:index+4]
00075 result.rotation, result.shift=_ParseRotationAndShift(next_lines)
00076 found_transform_section=True
00077 break
00078 assert found_transform_section and found_gdt_section
00079 return result
00080
00081 def GDT(pdb1, pdb2, chain1='', chain2='', reference_length=None, lga_bin=None):
00082 """
00083 Calculate GDT value between pdb1 and pdb2. It is assumed that the
00084 corresponding residues in pdb1 and pdb2 have the same residue numbers.
00085 """
00086 lga_bin=_FindLGABinary(lga_bin)
00087 temp_d=tempfile.mkdtemp(prefix='lga_gdt_ts_')
00088 for d in ('MOL2', 'TMP', 'RESULTS'):
00089 os.mkdir(os.path.join(temp_d, d))
00090
00091 if not chain1:
00092 chain1=pdb1.chains[0].name
00093 if not chain2:
00094 chain2=pdb2.chains[0].name
00095 pdb_one_name=os.path.join(temp_d, 'MOL2', 'one.pdb')
00096 pdb_two_name=os.path.join(temp_d, 'MOL2', 'two.pdb')
00097 io.SaveEntity(pdb1, pdb_one_name)
00098 io.SaveEntity(pdb2, pdb_two_name)
00099 _PrepareInput(pdb_one_name, pdb_two_name,
00100 os.path.join(temp_d, 'MOL2', 'input.pdb'))
00101 output_file=os.path.join(temp_d, 'out.txt')
00102 ch1, ch2=('', '')
00103 if len(chain1.strip()): ch1=' -ch1:%s ' % chain1
00104 if len(chain2.strip()): ch2=' -ch2:%s ' % chain2
00105
00106 params=(temp_d, lga_bin, 'input.pdb', ch1, ch2)
00107 command='cd %s; %s %s %s%s-ie -3 -d:4 -sda'
00108
00109 expanded_cmd=command % params
00110 lga_proc=subprocess.Popen(expanded_cmd, shell=True,
00111 stdout=subprocess.PIPE)
00112 length=reference_length or max(pdb1.residue_count, pdb2.residue_count)
00113 result=_ParseLGAOutput(lga_proc.stdout.readlines(), reference_length)
00114 os.system('rm -r %s' % temp_d)
00115 return result
