Source Code for Module csb.apps.bfite

  1  """ 
  2  Python application for robust structure superposition of an ensemble of structures. 
  3  bfite models non-rigid displacements in protein ensembles with outlier-tolerant 
  4  probability distributions. 
  5  """ 
  6  import numpy 
  7   
  8  import csb.apps 
  9  import csb.bio.structure 
 10   
 11  from csb.bio.io.wwpdb import LegacyStructureParser 
 12  from csb.bio.utils import average_structure, fit, wfit 
 13  from csb.statistics.scalemixture import ScaleMixture, GammaPrior 
 14   
 15   
 16 -class ExitCodes(csb.apps.ExitCodes): 
 17      IO_ERROR = 2 
 18   
 19 -class AppRunner(csb.apps.AppRunner): 
 20   
 21      @property 
 22 -    def target(self): 
 23          return BFitApp 
 24   
 25 -    def command_line(self): 
 26           
 27          cmd = csb.apps.ArgHandler(self.program, __doc__) 
 28   
 29          # Input structures 
 30          cmd.add_positional_argument('pdb', str, 
 31                                      'full path to the ensemble') 
 32   
 33          # Optional arguments  
 34          cmd.add_scalar_option('chain', 'c', str, 
 35                                'Chain', 
 36                                default='A') 
 37           
 38          cmd.add_scalar_option('scalemixture', 's', str, 
 39                                'Scale mixture distribution', 
 40                                default='student', 
 41                                choices=['student', 'k']) 
 42           
 43          cmd.add_scalar_option('alignment', 'a', str, 
 44                                'Alignment in fasta format defining equivalent positions\n' 
 45                                + 'Assumes that chain1 is the first sequence of ' 
 46                                + 'the alignment and chain2 the second sequence') 
 47   
 48          cmd.add_scalar_option('outfile', 'o', str, 
 49                                'file to which the rotated second ' +  
 50                                'structure will be written', 
 51                                default='bfit.pdb') 
 52   
 53          cmd.add_scalar_option('niter', 'n', int, 
 54                                'Number of optimization steps', 
 55                                default=200) 
 56           
 57          return cmd 
 58   
 59   
 60   
 61 -class BFitApp(csb.apps.Application): 
 62      """ 
 63      Python application for robust structure superposition of two protein structures 
 64      """ 
 65   
 66 -    def main(self): 
 67          try: 
 68              parser = LegacyStructureParser(self.args.pdb) 
 69              models = parser.models() 
 70   
 71          except IOError as e: 
 72              self.exit('PDB file parsing failed\n' + str(e.value), ExitCodes.IO_ERROR) 
 73   
 74          if len(models) < 2: 
 75              self.exit('PDB file contains only one model', ExitCodes.USAGE_ERROR) 
 76   
 77          ensemble = parser.parse_models(models) 
 78          X = numpy.array([model[self.args.chain].get_coordinates(['CA'], True) for model in ensemble]) 
 79          x_mu = average_structure(X) 
 80          #n = X.shape[1] 
 81          m = X.shape[0] 
 82          R = numpy.zeros((m, 3, 3)) 
 83          t = numpy.ones((m, 3)) 
 84   
 85   
 86          prior = GammaPrior() 
 87          mixture = ScaleMixture(scales=X.shape[1], 
 88                                 prior=prior, d=3) 
 89   
 90          for i in range(m): 
 91              R[i, :, :], t[i, :] = fit(x_mu, X[i]) 
 92           
 93          # gibbs sampling cycle 
 94          for j in range(self.args.niter): 
 95              # apply rotation 
 96              data = numpy.array([numpy.sum((x_mu - numpy.dot(X[i], numpy.transpose(R[i])) - t[i]) ** 2, -1) ** 0.5 
 97                                  for i in range(m)]).T 
 98              # sample scales 
 99              mixture.estimate(data) 
100              # sample rotations 
101              for i in range(m): 
102                  R[i, :, :], t[i, :] = wfit(x_mu, X[i], mixture.scales) 
103   
104   
105          out_ensemble = csb.bio.structure.Ensemble() 
106   
107          for i, model in enumerate(ensemble): 
108              model.transform(R[i], t[i]) 
109              out_ensemble.models.append(model) 
110   
111          out_ensemble.to_pdb(self.args.outfile) 
112   
113   
114 -def main(): 
115      AppRunner().run() 
116       
117       
118  if __name__ == '__main__': 
119      main() 
120