Package csb :: Package apps :: Module promix |
[frames] | no frames]

# Source Code for Module csb.apps.promix

      1  """
      2  ProMix: Take a protein structure ensemble and find a mixture of rigid
      3  segments or a mixture of conformers. Writes K copies of the ensemble
      4  (for segments) or K subsets of the ensemble (for conformers) as PDB
      5  files, each superposed on different components.
      6
      7  Reference: Hirsch M, Habeck M. - Bioinformatics. 2008 Oct 1;24(19):2184-92
      8  """
      9
     10  import numpy
     11
     12  import csb.apps
     13  import csb.bio.structure
     14
     15  from csb.bio.io.wwpdb import LegacyStructureParser
     16  from csb.statistics import mixtures
     17
     18
     19 -class ExitCodes(csb.apps.ExitCodes):
     20      IO_ERROR = 2
     21
     22 -class AppRunner(csb.apps.AppRunner):
     23
     24      @property
     25 -    def target(self):
     26          return ProMixApp
     27
     28 -    def command_line(self):
     29          cmd = csb.apps.ArgHandler(self.program, __doc__)
     30
     31          cmd.add_scalar_option('components', 'K', int, 'Number of components', -1)
     32          cmd.add_scalar_option('type', 't', str, 'Type of mixture', 'segments', ('segments', 'conformers'))
     33          cmd.add_positional_argument('infile', str, 'input PDB file')
     34
     35          return cmd
     36
     37 -    def initapp(self, args):
     38          app = self.target
     39          return app(args)
     40
     41 -class ProMixApp(csb.apps.Application):
     42
     43 -    def main(self):
     44          try:
     45              parser = LegacyStructureParser(self.args.infile)
     46              models = parser.models()
     47          except:
     48              self.exit('PDB file parsing failed', ExitCodes.IO_ERROR)
     49
     50          if len(models) < 2:
     51              self.exit('PDB file contains only one model', ExitCodes.USAGE_ERROR)
     52
     53          ensemble = parser.parse_models(models)
     54          X = numpy.array([model.get_coordinates(['CA'], True) for model in ensemble])
     55
     56          if self.args.type == 'segments':
     57              self.main_segments(ensemble, X)
     58          elif self.args.type == 'conformers':
     59              self.main_conformers(ensemble, X)
     60          else:
     61              raise ValueError('type must be "segments" or "conformers"')
     62
     63 -    def main_segments(self, ensemble, X):
     64
     65          mixture = mixtures.SegmentMixture.new(X, self.args.components)
     66          self.log('Number of segments: {0}'.format(mixture.K))
     67
     68          for k,(sigma,w) in enumerate(zip(mixture.sigma, mixture.w)):
     69              outfile = 'promix_segment_{0}.pdb'.format(k+1)
     70              self.log('  {0}: sigma = {1:6.3f}, w = {2:.3f}, file = {3}'.format(k+1, sigma, w, outfile))
     71
     72              for model, R, t in zip(ensemble, mixture.R, mixture.t):
     73                  if k > 0:
     74                      model.transform(R[k-1], t[k-1])
     75                  R = R[k].T
     76                  t = -numpy.dot(R, t[k])
     77                  model.transform(R, t)
     78
     79              ensemble.to_pdb(outfile)
     80
     81 -    def main_conformers(self, ensemble, X):
     82
     83          mixture = mixtures.ConformerMixture.new(X, self.args.components)
     84          self.log('Number of conformers: {0}'.format(mixture.K))
     85
     86          membership = mixture.membership
     87
     88          for k,(sigma,w) in enumerate(zip(mixture.sigma, mixture.w)):
     89              outfile = 'promix_conformer_{0}.pdb'.format(k+1)
     90              self.log('  {0}: sigma = {1:6.3f}, w = {2:.3f}, file = {3}'.format(k+1, sigma, w, outfile))
     91
     92              ek = csb.bio.structure.Ensemble()
     93
     94              for model, R, t, mk in zip(ensemble, mixture.R, mixture.t, membership):
     95                  if mk != k:
     96                      continue
     97                  R = R[k].T
     98                  t = -numpy.dot(R, t[k])
     99                  model.transform(R, t)
    100                  ek.models.append(model)
    101
    102              ek.to_pdb(outfile)
    103
    104
    105 -def main():
    106      AppRunner().run()
    107
    108
    109  if __name__ == '__main__':
    110      main()
    111