Source Code for Module csb.bio.io.noe

  1  """ 
  2  Simple XEASY and Sparky peak list parsers. 
  3  """ 
  4   
  5  from abc import ABCMeta, abstractmethod 
  6  from csb.bio.nmr import NOESpectrum 
  7   
  8   
  9 -class PeakListFormatError(ValueError): 
 10      pass 
 11   
 12 -class BasePeakListReader(object): 
 13       
 14      __metaclass__ = ABCMeta 
 15       
 16      @abstractmethod 
 17 -    def read(self, table): 
 18          """ 
 19          Parse a peak list table. 
 20           
 21          @param table: input peak list table 
 22          @type table: str 
 23          @rtype: L{NOESpectrum} 
 24          """ 
 25          pass 
 26       
 27 -    def read_file(self, filename): 
 28          """ 
 29          Parse a peak list file. 
 30           
 31          @param filename: input file name 
 32          @type filename: str 
 33          @rtype: L{NOESpectrum} 
 34          """ 
 35          with open(filename) as input: 
 36              return self.read(input.read()) 
 37           
 38 -    def read_all(self, filenames): 
 39          """ 
 40          Parse a list of peak list files and merge the resulting spectra. 
 41          All spectra must have identical dimensions. 
 42           
 43          @param filenames: input file names 
 44          @type filenames: iterable of str 
 45           
 46          @return: joint spectrum 
 47          @rtype: L{NOESpectrum} 
 48          """ 
 49          spectra = [self.read_file(f) for f in filenames] 
 50          return NOESpectrum.join(*spectra) 
 51           
 52 -class SparkyPeakListReader(BasePeakListReader): 
 53      """ 
 54      Sparky NOE peak list parser. 
 55       
 56      @param elements: list of element names for each dimension 
 57      @type elements: list of (str or L{EnumItem}) 
 58      @param connected: list of covalently connected dimension indices in the 
 59                        format: [(i1,i2),...] 
 60      @type connected: list of (int,int) tuples 
 61      """ 
 62       
 63 -    def __init__(self, elements, connected): 
 64           
 65          self._elements = list(elements) 
 66          self._connected = [(d1, d2) for d1, d2 in connected] 
 67           
 68          if len(self._elements) < 1: 
 69              raise ValueError("Can't parse a 0-dimensional peak list") 
 70       
 71 -    def read(self, table): 
 72          """ 
 73          Parse a Sparky peak list table. 
 74           
 75          @param table: input peak list 
 76          @type table: str 
 77          @rtype: L{NOESpectrum} 
 78          """ 
 79          offset = 0                 
 80          spectrum = NOESpectrum(self._elements) 
 81           
 82          for d1, d2 in self._connected: 
 83              spectrum.connect(d1, d2) 
 84           
 85          for l in table.splitlines(): 
 86              if not l.strip() or ('w1' in l and 'w2' in l): 
 87                  if l.lstrip().lower().startswith('assignment'): 
 88                      offset = 1 
 89                  continue 
 90               
 91              line = l.split()[offset:] 
 92              try: 
 93                  float(line[-1])             # last item may or may not be a comment 
 94              except ValueError: 
 95                  if len(line) > 0: 
 96                      line.pop() 
 97               
 98              items = list(map(float, line)) 
 99              intensity = items[-1] 
100              dimensions = items[:-1] 
101               
102              if len(dimensions) != len(self._elements): 
103                  raise PeakListFormatError("Expected {0} dimensional spectrum, got {1}".format( 
104                                                                      len(self._elements), len(dimensions))) 
105               
106              spectrum.add(intensity, dimensions) 
107           
108          return spectrum        
109                       
110 -class XeasyPeakListReader(BasePeakListReader): 
111      """ 
112      XEASY NOE peak list parser. 
113      """ 
114           
115 -    def __init__(self): 
116          pass 
117               
118 -    def read(self, table): 
119          """ 
120          Parse an XEASY peak list table. 
121           
122          @param table: input peak list 
123          @type table: str 
124          @rtype: L{NOESpectrum} 
125          """       
126          lines = table.splitlines()   
127          spectrum = self._read_header(lines) 
128           
129          for l in lines: 
130              if not l.strip() or l.startswith('#'): 
131                  continue 
132               
133              parts = l.split()[1:] 
134              peak = parts[:spectrum.num_dimensions] 
135              height = parts[spectrum.num_dimensions + 2]  
136               
137              intensity = float(height) 
138              dimensions = map(float, peak) 
139               
140              spectrum.add(intensity, dimensions) 
141           
142          return spectrum            
143       
144       
145 -    def _read_header(self, lines): 
146           
147          num = 0 
148          dim = {} 
149          el = {} 
150          el2 = {} 
151          connectivity = None 
152           
153          for l in lines: 
154              if l.startswith('#'): 
155                  if l[1:].lstrip().lower().startswith('number of dimensions'): 
156                      num = int(l.split()[-1]) 
157                   
158                  if l.startswith('#INAME'): 
159                      parts = l.split()[1:] 
160                      if len(parts) != 2: 
161                          raise PeakListFormatError("Invalid Xeasy header")                     
162                       
163                      index = int(parts[0]) - 1 
164                      if index < 0: 
165                          raise PeakListFormatError("Invalid Xeasy header: dimension index < 1") 
166                       
167                      element = ''.join(i for i in parts[1] if i.isalpha()) 
168                      el[parts[1]] = index 
169                      el2[element] = index 
170                        
171                      dim[index] = element 
172                       
173                  if l.startswith('#CYANAFORMAT'): 
174                      connectivity = l.split()[1] 
175         
176          if len(dim) != num or num == 0: 
177              raise PeakListFormatError("Invalid Xeasy header") 
178           
179          elements = tuple(dim[i] for i in sorted(dim)) 
180          spectrum = NOESpectrum(elements) 
181           
182          if connectivity: 
183              # cyanaformat - explicitly defines connected dimensions: 
184              # upper case dimensions are connected, e.g. "#CYANAFORMAT hCH" => 2-3  
185              if connectivity.upper() != ''.join(elements).upper(): 
186                  raise ValueError("Invalid XEASY/CYANA header")  
187              for i1 in range(len(connectivity)): 
188                  for i2 in range(len(connectivity)): 
189                      e1, e2 = connectivity[i1], connectivity[i2]  
190                      if i1 != i2 and e1.isupper() and e2.isupper(): 
191                          spectrum.connect(i1, i2)                         
192          else: 
193              # dimension labels starting with a number are connected, e.g. "1A B2 3C" => 1-3 
194              if len(el) != num: 
195                  raise PeakListFormatError("Invalid XEASY header")                 
196              for e1 in el: 
197                  for e2 in el: 
198                      if e1 != e2: 
199                          element1 = dim[el[e1]] 
200                          element2 = dim[el[e2]] 
201                           
202                          num1 = e1.replace(element1, '') 
203                          num2 = e2.replace(element2, '') 
204                           
205                          if e1.startswith(num1) and e2.startswith(num2): 
206                              spectrum.connect(el[e1], el[e2]) 
207                           
208          return spectrum 
209   
210   
211 -class XeasyFileBuilder(object): 
212      """ 
213      XEASY output format builder. 
214       
215      @param stream: destination stream, were the output is written 
216      @type stream: file 
217      """ 
218       
219 -    def __init__(self, stream): 
220          self._out = stream 
221           
222 -    def add_spectrum(self, spectrum): 
223           
224          self.add_header(spectrum) 
225          self.add_peaks(spectrum) 
226           
227 -    def add_header(self, spectrum): 
228          """ 
229          Write the XEASY header. 
230           
231          @param spectrum: NOE spectrum 
232          @type spectrum: L{NOESpectrum}  
233          """ 
234           
235          self._out.write( 
236              '# Number of dimensions {0}\n'.format(spectrum.num_dimensions)) 
237           
238          conn = '' 
239           
240          for en, e in enumerate(spectrum.dimensions, start=1): 
241              element = repr(e).upper() 
242              self._out.write('#INAME {0} {1}{0}\n'.format(en, element)) 
243               
244              if spectrum.has_connected_dimensions(en - 1): 
245                  conn += element.upper() 
246              else: 
247                  conn += element.lower() 
248                   
249          self._out.write( 
250              '#CYANAFORMAT {0}\n'.format(conn)) 
251       
252 -    def add_peaks(self, spectrum): 
253          """ 
254          Write all peaks from C{spectrum}. 
255           
256          @param spectrum: NOE spectrum 
257          @type spectrum: L{NOESpectrum}  
258          """ 
259           
260          for pn, peak in enumerate(spectrum, start=1): 
261              self._out.write("{0:5} ".format(pn)) 
262               
263              for dim in range(spectrum.num_dimensions): 
264                  data = "{0:7.3f} ".format(peak.get(dim)) 
265                  self._out.write(data) 
266                               
267              self._out.write("2 U ") 
268              self._out.write("{0:18e} ".format(peak.intensity)) 
269              self._out.write("0.00e+00 m   0    0    0    0 0\n") 
270