> although I see that there is still there are separate isotope issues in
> the XML.
Attached are the latest updates to the CCPN:ARIA2 interface (for src/py/).
I have added the passing back of rejected peaks. These are suitably marked
in the details field and have a merit of 0.0, so would not be used again
in ARIA or in Analysis restraint list generation by default.
The updates also include further modifications to prevent the isotope
mismatch issue. As I may have mentioned before, the bigger problem due to
different dimension orders was sorted with the early updates. However,
occasionally it seems that an assignment possibility with a coincidentally
matching shift but wrong isotope type seems to be getting through. For
example a Phe Hd/Cd peak may match backbone H/N resonances. I can't see
any obvious problems with the CCPN interface so it might be coming from
ARIA. In any case, there is extra, explicit isotope checking in place to
stop any aberrant resonance assignments getting through to the CCPN
project, whatever their source.
T.
-------------------------------------------------------------------------------
Dr Tim Stevens Email: [log in to unmask]
Department of Biochemistry [log in to unmask]
University of Cambridge Phone: +44 1223 766018 (office)
80 Tennis Court Road +44 7816 338275 (mobile)
Old Addenbrooke's Site +44 1223 364613 (home)
Cambridge CB2 1GA WWWeb: http://www.bio.cam.ac.uk/~tjs23
United Kingdom http://www.pantonia.co.uk
-------------------------------------------------------------------------------
------ +NH3CH(CH(CH3)OH)C(O)NHCH(CH(CH3)CH2CH3)C(O)NHCH(CH2CH2SCH3)CO2- -------
-------------------------------------------------------------------------------
"""
ARIA -- Ambiguous Restraints for Iterative Assignment
A software for automated NOE assignment
Version 2.2
Copyright (C) Benjamin Bardiaux, Michael Habeck, Therese Malliavin,
Wolfgang Rieping, and Michael Nilges
All rights reserved.
NO WARRANTY. This software package is provided 'as is' without warranty of
any kind, expressed or implied, including, but not limited to the implied
warranties of merchantability and fitness for a particular purpose or
a warranty of non-infringement.
Distribution of substantively modified versions of this module is
prohibited without the explicit permission of the copyright holders.
$Author: bardiaux $
$Revision: 1.14 $
$Date: 2007/11/26 14:23:53 $
"""
from aria import *
from Settings import Settings
from xmlutils import XMLElement, XMLBasePickler
from N import *
from time import clock
from time import ctime
## BARDIAUX 280105
REPORT_UPDATED_SPECTRUM = '.assigned'
REPORT_SUMMARY = 'report'
REPORT_NOE_RESTRAINTS = 'noe_restraints'
MISSING_STRUCTURES = 'missing_structures'
KEYBOARD_INTERRUPT = 'keyboard_interrupt'
## local path where MOLMOL restraints files are stored
PATH_MOLMOL = 'molmol'
def files_exist(file_list):
import os
for fn in file_list:
fn = os.path.expanduser(fn)
if not os.path.exists(fn):
return 0
return 1
def dump_peaks_as_text(peak_list, filename, gzip, options = None):
check_type(peak_list, LIST, TUPLE)
check_string(filename)
check_int(gzip)
check_type(options, NONE, DICT)
import AriaPeak
## write contributions as text file (both, ambiguous and
## unambiguous in one file)
## TODO: so far, text-pickler options are hard-coded.
settings = AriaPeak.AriaPeakListTextPicklerSettings()
## set options
if options is not None:
settings.update(options)
pickler = AriaPeak.AriaPeakListTextPickler(settings)
## sort peak list such, that active restraints come first,
## inactive last.
active = [p for p in peak_list if p.isActive()]
inactive = [p for p in peak_list if not p.isActive()]
peak_list = active + inactive
## Dump assignments
pickler.dump_assignments(peak_list, filename + '.assignments')
## decompose peak_list into ambiguous and unambiguous
## restraints
ambig = [p for p in peak_list if p.isAmbiguous()]
unambig = [p for p in peak_list if not p.isAmbiguous()]
pickler.dump_ambiguous(ambig, filename + '.ambig', gzip)
pickler.dump_unambiguous(unambig, filename + '.unambig', gzip)
## Dump restraint list, sorted by violations.
## Active restraints come first, inactive last
violated = [p for p in peak_list if p.analysis.isViolated()]
if violated:
d_lower = []
d_upper = []
lower = []
upper = []
for r in violated:
## In case of a lower-bound violation,
## put restraint in list of lower-bound-violating
## restraints
d = r.analysis.getLowerBoundViolation()[0]
if d is not None and d > 0.:
d_lower.append(d)
lower.append(r)
else:
d = r.analysis.getUpperBoundViolation()[0]
if d is not None and d > 0.:
d_upper.append(d)
upper.append(r)
lower_indices = argsort(d_lower).tolist()
lower_indices.reverse()
upper_indices = argsort(d_upper).tolist()
upper_indices.reverse()
a = [lower[i] for i in lower_indices]
b = [upper[i] for i in upper_indices]
viol = a + b
else:
viol = []
active = [p for p in viol if p.isActive()]
inactive = [p for p in viol if not p.isActive()]
viol = active + inactive
pickler.dump_violations(viol, filename + '.violations', gzip)
class IterationSettings(Settings):
def create(self):
from Settings import TypeEntity, ChoiceEntity, NonNegativeInt
kk = {}
kk['number'] = NonNegativeInt()
kk['number_of_structures'] = NonNegativeInt()
kk['number_of_best_structures'] = NonNegativeInt()
choices = ['total_energy', 'restraint_energy','noe_violations', 'restraint_violations'] # BARDIAUX 2.2
msg = 'Sort structures according to one of the following' + \
'choices: %s' % ', '.join(choices)
entity = ChoiceEntity(choices, error_message = msg)
kk['sort_criterion'] = entity
msg = 'Argument has to be a %s instance.'
d = {'calibrator_settings': 'CalibratorSettings',
'peak_assigner_settings': 'PeakAssignerSettings',
'violation_analyser_settings': 'ViolationAnalyserSettings',
'contribution_assigner_settings': 'ContributionAssignerSettings',
'merger_settings': 'MergerSettings',
'network_anchoring_settings' : 'NetworkSettings'}
for name, type in d.items():
entity = TypeEntity(type, error_message = msg % type)
kk[name] = entity
return kk
def create_default_values(self):
import Calibrator, PeakAssigner
import ViolationAnalyser as VA
import ContributionAssigner as CA
import Merger
## BARDIAUX 2.2
import Network
d = {}
d['calibrator_settings'] = Calibrator.CalibratorSettings()
d['peak_assigner_settings'] = PeakAssigner.PeakAssignerSettings()
d['violation_analyser_settings'] = VA.ViolationAnalyserSettings()
d['contribution_assigner_settings'] = CA.ContributionAssignerSettings()
d['merger_settings'] = Merger.MergerSettings()
## BARDIAUX 2.2
d['network_anchoring_settings'] = Network.NetworkSettings()
for value in d.values():
value.reset()
d['number'] = 0
d['number_of_structures'] = 20
d['number_of_best_structures'] = 7
d['sort_criterion'] = 'total_energy'
return d
class IterationSettingsXMLPickler(XMLBasePickler):
order = ['number', 'n_structures', 'sort_criterion',
'n_best_structures', 'assignment', 'merging',
'calibration', 'violation_analysis',
'partial_assignment', 'network_anchoring']
def _xml_state(self, x):
e = XMLElement(tag_order = self.order)
e.number = x['number']
e.n_structures = x['number_of_structures']
e.sort_criterion = x['sort_criterion']
e.n_best_structures = x['number_of_best_structures']
e.assignment = x['peak_assigner_settings']
e.merging = x['merger_settings']
e.calibration = x['calibrator_settings']
e.violation_analysis = x['violation_analyser_settings']
e.partial_assignment = x['contribution_assigner_settings']
# BARDIAUX 2.2
e.network_anchoring = x['network_anchoring_settings']
return e
def load_from_element(self, e):
s = IterationSettings()
s['number'] = int(e.number)
s['number_of_structures'] = int(e.n_structures)
s['sort_criterion'] = str(e.sort_criterion)
s['number_of_best_structures'] = int(e.n_best_structures)
s['peak_assigner_settings'] = e.assignment
s['merger_settings'] = e.merging
s['calibrator_settings'] = e.calibration
s['violation_analyser_settings'] = e.violation_analysis
s['contribution_assigner_settings'] = e.partial_assignment
## BARDIAUX 2.2
if hasattr(e, 'network_anchoring'):
s['network_anchoring_settings'] = e.network_anchoring
else:
import Network
nas = Network.NetworkSettings()
nas.reset()
s['network_anchoring_settings'] = nas
return s
class ProtocolSettings(Settings):
def create(self):
from Settings import TypeEntity, ChoiceEntity
d = {}
## water-refinement parameters
data_type = 'WaterRefinementParameters'
msg = 'Argument has to be a %s instance.' % data_type
entity = TypeEntity(data_type, error_message = msg)
d['water_refinement'] = entity
## iteration settings
entity = TypeEntity(DICT)
entity.set({})
d['iteration_settings'] = entity
## floating assignment flag
choices = [YES, NO]
msg = 'Floating assignment: %s?' % ' / '.join(choices)
entity = ChoiceEntity(choices, error_message = msg)
d['floating_assignment'] = entity
return d
def create_default_values(self):
return {'floating_assignment': YES}
def addIterationSettings(self, settings):
check_type(settings, 'IterationSettings')
number = settings['number']
it_settings = self['iteration_settings']
if not it_settings.has_key(number):
it_settings[number] = settings
else:
msg = 'IterationSettings instance with same number ' + \
'(%d) already stored in protocol settings.' % number
self.error(KeyError, msg)
def delIterationSettings(self, s):
check_type(s, 'IterationSettings')
number = s['number']
it_settings = self['iteration_settings']
if it_settings.has_key(number):
del it_settings[number]
else:
msg = 'IterationSettings (number %d) d not exist.' % number
self.error(KeyError, msg)
class Protocol(AriaBaseClass):
def __init__(self, settings):
check_type(settings, 'ProtocolSettings')
AriaBaseClass.__init__(self, name = 'Protocol')
self.setSettings(settings)
self.__reports_written = 0
self.constraintSetSerial = None
## self.iterations = {}
def getIterationSettings(self, number):
check_int(number)
return self.getSettings()['iteration_settings'][number]
def getInfrastructure(self):
from Singleton import ProjectSingleton
project = ProjectSingleton()
return project.getInfrastructure()
def initialize(self):
"""
must be called after unpickling a protocol in
order to finalize its initialization.
create objects which implement the actual aria-method
set parameters of first iteration.
"""
from NOEModel import NOEModel, ISPA
import Calibrator
import ViolationAnalyser as VA
import ContributionAssigner as CA
import Merger
## BARDIAUX 2.2
import Network
## NOE model
## TODO: not nessecary to store model!
#self.model = NOEModel()
self.model = ISPA()
## calibrator
cs = Calibrator.CalibratorSettings()
self.calibrator = Calibrator.Calibrator(cs, self.model)
## violation-analyser
vas = VA.ViolationAnalyserSettings()
self.violation_analyser = VA.ViolationAnalyser(vas)
## contribution assigner
cas = CA.ContributionAssignerSettings()
self.contribution_assigner = CA.ContributionAssigner(cas)
## merger
ms = Merger.MergerSettings()
self.merger = Merger.Merger(ms)
## BARDIAUX 2.2
## network
nas = Network.NetworkSettings()
self.network_anchoring = Network.NetworkAnchoring(nas)
def _updateIterationSettings(self, iteration):
"""
applies the parameters-settings in 'iteration' to
sub-modules
"""
check_type(iteration, 'Iteration')
settings = self.getIterationSettings(iteration.getNumber())
cs = settings['calibrator_settings']
self.calibrator.getSettings().update(cs)
vas = settings['violation_analyser_settings']
self.violation_analyser.getSettings().update(vas)
cas = settings['contribution_assigner_settings']
self.contribution_assigner.getSettings().update(cas)
ms = settings['merger_settings']
self.merger.getSettings().update(ms)
## BARDIAUX 2.2
nas = settings['network_anchoring_settings']
self.network_anchoring.getSettings().update(nas)
def _updateSpectrumSettings(self, spectrum):
spectrum_data = spectrum.getDataSource()
self.calibrator.getSettings().update(spectrum_data)
self.violation_analyser.getSettings().update(spectrum_data)
def mergePeakLists(self, iteration):
## BARDIAUX 2.2
# remove previous combined peaks
iteration.setCombinedRestraints([])
if self.merger.getSettings()['method'] == 'no_merging':
return
peaks = [p for p in iteration.getPeakList() if p.isActive()]
## merge peaks. peaks which are merged (away), i.e.
## are marked via peak.isMerged(1)
n, combined = self.merger(peaks)
# add combined restraints
iteration.setCombinedRestraints(combined)
if n:
msg = 'Peaks merged: %d / %d peaks (%.1f %%) (see report for details).'
self.message(msg % (n, len(peaks), 100. * n / len(peaks)))
if combined:
msg = 'Peaks combined: %d / %d peaks (%.1f %%).'
self.message(msg % (len(combined), len(peaks), 100. * len(combined) / len(peaks)))
## self.message('old %d new %d.' % (len(peaks), len([p for p in iteration.getPeakList() if p.isActive()])))
def compileFileList(self, number, **kw):
"""
Compiles a tuple of filenames which must exist after
iteration 'number' has been successfully completed.
If the method is called with keywords, the returned tuple
contains only a sub-set of all filenames. The sub-set depends
on the particular keyword(s)
"""
check_int(number)
STRUCTURES = 'structures'
FLOAT_FILES = 'float_files'
from os.path import join
keywords = (STRUCTURES, FLOAT_FILES)
unknown = [key for key in kw if key not in keywords]
if unknown:
s = 'Unknown sub-list identifier ("%s"). Known identifier ' + \
'are: %s'
self.error(s % (str(unknown), str(keywords)))
if kw:
keywords = kw.keys()
infra = self.getInfrastructure()
## structures (pdb-files)
it_settings = self.getIterationSettings(number)
n_structures = it_settings['number_of_structures']
path = infra.get_iteration_path(number)
l = []
if STRUCTURES in keywords:
name_template = infra.get_file_root() + '_%d.pdb'
l += [join(path, name_template % j) \
for j in range(1, n_structures + 1)]
if FLOAT_FILES in keywords:
name_template = infra.get_file_root() + '_%d.float'
l += [join(path, name_template % j) \
for j in range(1, n_structures + 1)]
return tuple(l)
def findFirstIteration(self):
"""
For every iteration, this method checks if the required
PDB-files exist. If so, it seeks to the next iteration until
it ends up at an iteration with missing files. If files are
present for all iterations, None is returned. Otherwise the
number of the iteration with missing files is returned.
"""
infra = self.getInfrastructure()
n_iterations = infra.getSettings()['n_iterations']
for i in range(n_iterations):
## Get filenames of all structures for iteration i
pdb_files = self.compileFileList(i, structures=1)
if not files_exist(pdb_files):
return i
return None
def _writePDBFileList(self, iteration):
import os
infra = self.getInfrastructure()
iteration_path = infra.get_iteration_path(iteration.getNumber())
molmol_path = os.path.join(iteration_path, PATH_MOLMOL)
ok = 1
if not os.path.exists(molmol_path):
try:
os.makedirs(molmol_path)
except Exception, msg:
import tools
self.warning(tools.last_traceback())
msg = 'Could not create directory for MOLMOL output.'
self.warning(msg)
ok = 0
if not ok:
return None
## file.nam
## files are sorted according to 'sort_criterion'
## which is part of the <structure_generation>
## section of project xml.
ensemble = iteration.getStructureEnsemble()
file_list = ensemble.getFiles()
if file_list is None:
s = 'Could not write MolMol file.nam: No PDB-files' + \
'in structure ensemble.'
self.warning(s)
return None
filename = os.path.join(molmol_path, 'file.nam')
try:
f = open(filename, 'w')
except:
self.warning('Could not create %s.' % filename)
return None
s = '\n'.join(file_list)
try:
f.write(s)
f.close()
except:
self.warning('Could not write PDB-file list (%s).' % filename)
return None
return molmol_path
def structure_calc_done(self, engine, molecule, iteration):
## Check if there are missing structures
if engine.missingStructures():
self.__done = MISSING_STRUCTURES
return
self.message('Structure calculation done.')
it_settings = self.getIterationSettings(iteration.getNumber())
ensemble = engine.getEnsemble(it_settings, molecule)
## store structure-ensemble in iteration
iteration.setStructureEnsemble(ensemble)
self.__done = 1
def startStructureCalculation(self, iteration, molecule):
from Singleton import ProjectSingleton
## get list of all restraints
peaks = iteration.getPeakList()
## BARDIAUX 2.2 Add user-distance CCPN restraints
restraints = iteration.getDistanceRestraintsList()
peaks += restraints
# add also combined restraints
combined = iteration.getCombinedRestraints()
peaks += combined
project = ProjectSingleton()
engine = project.getStructureEngine()
## callback that will be called when pdb-files have been generated.
f = lambda engine = engine, molecule = molecule, it = iteration, \
c = self.structure_calc_done: c(engine, molecule, it)
engine.set_callback(f)
engine.go(peaks, self.getSettings(), iteration.getNumber())
def doSolventRefinement(self, iteration, molecule):
from Singleton import ProjectSingleton
import StructureEnsemble as SE
## If the given iteration lacks a structure ensemble,
## load PDB-files from last iteration.
ensemble = iteration.getStructureEnsemble()
if ensemble is None:
## create structure-ensemble from given 'filenames'
number = iteration.getNumber()
it_settings = self.getIterationSettings(number)
se_settings = SE.StructureEnsembleSettings()
## Set iteration-specific settings
se_settings.update(it_settings)
ensemble = SE.StructureEnsemble(se_settings)
filenames = self.compileFileList(number, structures=1)
naming_convention = 'cns'
ensemble.read(filenames, molecule, naming_convention)
iteration.setStructureEnsemble(ensemble)
project = ProjectSingleton()
engine = project.getStructureEngine()
file_list = engine.solvent_refine(self.getSettings(), ensemble)
## cleanup (for solvent refinement and last iteration)
infra = self.getInfrastructure()
last_it = infra.getSettings()['n_iterations'] - 1
engine.cleanup(infra.get_refinement_path())
engine.cleanup(infra.get_iteration_path(last_it))
def doViolationAnalysis(self, restraints, ensemble, store_analysis = 0):
"""
'restraints': list of AriaPeaks. The bounds of every
restraint in that list is checked against distances found
in the 'ensemble'.
'targets': list of AriaPeaks. The violationAnalyser will
store all intermediate results in their analysis-section.
Note: we assume, that peaks[i] corresponds to results[i]
for all i !. If a restraint has been violated, the
corresponding 'target'_restraint will be marked as violated.
"""
f = self.violation_analyser.analysePeak
violated = []
non_violated = []
## get theshold for current iteration
settings = self.violation_analyser.getSettings()
threshold = settings['violation_threshold']
for restraint in restraints:
R_viol = f(ensemble, restraint, store_analysis)
##
## If a restraint has been violated in too many structures
## (according to 'threshold'), mark is a violated.
##
if R_viol > threshold:
restraint.analysis.isViolated(1)
violated.append(restraint)
else:
restraint.analysis.isViolated(0)
non_violated.append(restraint)
## For violated restraints: if bound-correction is enabled,
## repeat violation-analysis with modified bounds.
if settings['lower_bound_correction']['enabled'] == YES:
new_lower = settings['lower_bound_correction']['value']
else:
new_lower = None
if settings['upper_bound_correction']['enabled'] == YES:
new_upper = settings['upper_bound_correction']['value']
else:
new_upper = None
if new_lower is not None or new_upper is not None:
## We forget 'store_analysis' here, since it has already
## been stored (if set).
R_viol = [f(ensemble, r, lower_correction = new_lower,
upper_correction = new_upper) for r in violated]
## List of restraint-indices which are no longer
## violated after bound modification.
indices = nonzero(less(R_viol, threshold))
new_non_violated = [violated[i] for i in indices]
[r.analysis.isViolated(0) for r in new_non_violated]
else:
new_non_violated = None
return violated, non_violated, new_non_violated
def done(self):
self.message('finished.')
def _dump_peak_list(self, peaks, path):
check_list(peaks)
check_string(path)
import os
from Singleton import ProjectSingleton
import Merger
project = ProjectSingleton()
s = project.getReporter()['noe_restraint_list']
## Dump only non-merged restraitns
not_merged = [r for r in peaks if not r.isMerged()]
if s['xml_output'] in (YES, GZIP):
from AriaXML import AriaXMLPickler as pickler
filename = os.path.join(path, REPORT_NOE_RESTRAINTS + '.xml')
t = clock()
p = pickler()
# BARDIAUX 2.2: It's not possible to dump
# CCPN restraints as XML
not_merged = [r for r in not_merged if not is_type(r, 'DistanceRestraint')]
gzip = {YES: 0, GZIP: 1}[s['xml_output']]
p.dump(not_merged, filename, gzip = gzip)
self.message('NOE-restraint list (xml) written (%s).' % filename)
self.debug('Time: %ss' % str(clock()-t))
if s['text_output'] in (YES, GZIP):
t = clock()
filename = os.path.join(path, REPORT_NOE_RESTRAINTS)
gzip = {YES: 0, GZIP: 1}[s['text_output']]
dump_peaks_as_text(not_merged, filename, gzip)
self.message('NOE-restaint list (text) written (%s).' % filename)
self.debug('Time: %ss' % str(clock()-t))
if s['pickle_output'] in (YES, GZIP):
from tools import Dump
filename = os.path.join(path, REPORT_NOE_RESTRAINTS + '.pickle')
t = clock()
gzip = {YES: 0, GZIP: 1}[s['pickle_output']]
Dump(not_merged, filename, gzip = gzip)
self.message('NOE-restraint list (python pickle) ' + \
'written (%s).' % filename)
self.debug('Time: %ss' % str(clock()-t))
## Write Merger report
t = clock()
merged = [r for r in peaks if r.isMerged()]
filename = os.path.join(path, REPORT_NOE_RESTRAINTS + '.merged')
p = Merger.MergerTextPickler()
p.dump(merged, filename)
self.message('Report for Merging-step written (%s).' % filename)
self.debug('Time: %ss' % str(clock() - t))
def _dump_spectra(self, peaks, spectrum, filename):
check_list(peaks)
check_type(spectrum, 'NOESYSpectrum')
check_string(filename)
import os
from Singleton import ProjectSingleton
from NOESYSpectrum import NOESYSpectrum
import Assignment
from copy import deepcopy
from AriaXML import AriaXMLPickler as pickler
project = ProjectSingleton()
s = project.getReporter()['spectra']
if s['write_unambiguous_only'] == YES:
peaks = [pk for pk in peaks if not pk.isAmbiguous()]
if peaks and s['write_assigned'] == YES:
spectrum = deepcopy(spectrum)
for pk in peaks:
ref_pk = spectrum.findPeak(pk.getReferencePeak().getNumber())
if ref_pk is None:
self.error('Inconsistency: Could not find reference peak %d in copy of spectrum "%s"' %
(ref_pk.getNumber(), spectrum.getName()))
active_contribs = [c for c in pk.getContributions() \
if c.getWeight() > 0.]
if not ref_pk.isAssigned() or \
(ref_pk.isAssigned() and s['write_assigned_force'] == YES) and\
active_contribs:
h = active_contribs[0]
ss = h.getSpinSystems()
for sp_sys in ss:
dim = ref_pk.getDimension(sp_sys)
AUTO = Assignment.ASSIGNMENT_TYPE_AUTOMATIC
A = Assignment.Assignment(sp_sys.getAtoms(), AUTO)
if dim == 1:
ref_pk.setProton1Assignments((A,))
elif dim == 2:
ref_pk.setProton2Assignments((A,))
for h in active_contribs[1:]:
ss = h.getSpinSystems()
for sp_sys in ss:
dim = ref_pk.getDimension(sp_sys)
AUTO = Assignment.ASSIGNMENT_TYPE_AUTOMATIC
A = Assignment.Assignment(sp_sys.getAtoms(), AUTO)
if dim == 1:
atoms_assi = [a.getAtoms() for a in ref_pk.getProton1Assignments()]
if not A.getAtoms() in atoms_assi:
ref_pk.addProton1Assignment(A)
elif dim == 2:
atoms_assi = [a.getAtoms() for a in ref_pk.getProton2Assignments()]
if not A.getAtoms() in atoms_assi:
ref_pk.addProton2Assignment(A)
t = clock()
p = pickler()
p.dump(spectrum, filename, gzip = 0)
self.debug('Time: %ss' % str(clock()-t))
self.message('Assigned spectrum "%s" written to file %s"' \
% (spectrum.getName(), filename))
def _dump_ccpn(self, iteration, path, is_water_refinement=False):
from Singleton import ProjectSingleton
project = ProjectSingleton()
if project.ccpn_project_instance is None:
return
import exportToCcpn as ccpnExport
from importFromCcpn import getKeysFromString
project_settings = project.getSettings()
run_name = project_settings['run'].strip().replace(' ', '_')
file_root = project_settings['file_root'].strip().replace(' ', '_')
ccpn_project = project.ccpn_project_instance
nmrProject = ccpn_project.currentNmrProject or ccpn_project.findFirstNmrProject()
## TBD: proper handling of NMR projects, e.g. the case in which
## multiple NMR projects exist.
if not nmrProject:
self.message('CCPN export: No NMR project found, creating new one.')
name = 'ARIA2_run_%s' % run_name
nmrProject = ccpn_project.newNmrProject(name=name)
s = project.getReporter()['ccpn']
export = 0
# TJS: Includes inactive ones as these are
# now passed back to CCPN (and marked)
restraints = iteration.getPeakList()
## TBD: multiple chains
aria_chain = project.getMolecule().get_chains()[0]
if s['export_structures'] == YES and \
(iteration.getNumber() == project_settings['n_iterations']-1 or \
is_water_refinement):
export = 1
# BARDIAUX: Update for multiple chains
chains = ccpnExport.getChainsFromAria2(restraints, ccpn_project,
aria_chain=aria_chain)
if not chains:
self.warning(StandardError, 'CCPN export: No molecular system found.')
return
structures = ccpnExport.getStructuresFromAria2Dir(path, chains[0])
if not structures:
self.warning('CCPN export: Unable to load any structures from iteration directory %s' % path)
else:
if self.constraintSetSerial is not None:
constraintSet = nmrProject.findFirstNmrConstraintStore(serial=self.constraintSetSerial)
self.message('CCPN export: Using existing constraint set.')
else:
self.message('CCPN export: Creating new constraint set.')
constraintSet = ccpnExport.makeNmrConstraintStore(nmrProject)
self.constraintSetSerial = constraintSet.serial
self.message('CCPN export: PDB files exported.')
struct_gen = ccpnExport.makeStructureGeneration(structures, constraintSet)
if not is_water_refinement:
struct_gen.name = 'ARIA2_%s_run_%s_iteration_%d' % (file_root, run_name, iteration.getNumber())
struct_gen.details = 'Structures created by ARIA2, run "%s", iteration %d.' % \
(project.getSettings()['run'], iteration.getNumber())
else:
struct_gen.name = 'ARIA2_%s_run_%s_water_refinement' % (file_root, run_name)
struct_gen.details = 'Structures created by ARIA2, run "%s", water refinement.'
else:
structures = None
if not is_water_refinement and (s['export_noe_restraint_list'] == 'all' or \
(s['export_noe_restraint_list'] == 'last' and \
iteration.getNumber() == project.getSettings()['n_iterations']-1)):
export = 1
if self.constraintSetSerial is not None:
constraintSet = nmrProject.findFirstNmrConstraintStore(serial=self.constraintSetSerial)
self.message('CCPN export: Using existing constraint set.')
else:
self.message('CCPN export: Creating new constraint set.')
constraintSet = ccpnExport.makeNmrConstraintStore(nmrProject)
self.constraintSetSerial = constraintSet.serial
# BARDIAUX: Update for multiple chains
chains = ccpnExport.getChainsFromAria2(restraints, ccpn_project,
aria_chain=aria_chain)
if restraints:
# TJS expert rejected (not active) restraints back to CCPN for analysis
constrList, rejectList, violList = ccpnExport.getConstraintsFromAria2(restraints, chains, constraintSet,
structures)
constrList.name = 'ARIA2_%s_run%s_NOEs_it%d' % (file_root, run_name, iteration.getNumber())
constrList.details = 'ARIA2 NOE restraint list, project "%s", run "%s", iteration %d.' % \
(project.getSettings()['name'], project.getSettings()['run'], iteration.getNumber())
# TJS constraint list contaning the inactive restraints
if rejectList is not None:
rejectList.name = 'ARIA2_REJECT_%s_run%s_NOEs_it%d' % (file_root, run_name, iteration.getNumber())
rejectList.details = 'ARIA2 NOE *INACTIVE* restraint list, project "%s", run "%s", iteration %d.' % \
(project.getSettings()['name'], project.getSettings()['run'], iteration.getNumber())
self.message('CCPN export: NOE restraint list exported.')
# BARDIAUX export Constraints from CCPN
for orig_list, distance_constraints in iteration.getDistanceRestraints().items():
# TJS adjust; we want the inactive now. ;-)
# distance_constraints = [r for r in distance_constraints if r.isActive()]
constrList, rejectList, violList = ccpnExport.getConstraintsFromAria2(distance_constraints, chains, constraintSet,
structures)
# TJS adjust
ccpnKeys = getKeysFromString(orig_list.getDataSource().get('ccpn_id'))
print "ccpnKeys", ccpnKeys
list_id = '%s|%s' % (ccpnKeys[-2],ccpnKeys[-1])# Project id is too verbose and redundant inside same project
constrList.name = 'ARIA2_%s_run%s_Dists_%s_it%d' % (file_root, run_name, list_id, iteration.getNumber())
constrList.details = 'ARIA2 DistanceConstraints list %s, project "%s", run "%s", iteration %d.' % \
(orig_list.getName(), project.getSettings()['name'], project.getSettings()['run'], iteration.getNumber())
# TJS constraint list contaning the inactive restraints
if rejectList is not None:
rejectList.name = 'ARIA2_REJECT_%s_run%s_Dists_%s_it%d' % (file_root, run_name, list_id, iteration.getNumber())
rejectList.details = 'ARIA2 *INACTIVE* DistanceConstraints list %s, project "%s", run "%s", iteration %d.' % \
(orig_list.getName(), project.getSettings()['name'], project.getSettings()['run'], iteration.getNumber())
self.message('CCPN export: DistanceConstraints list exported.')
if not is_water_refinement and (s['export_assignments'] == YES and \
iteration.getNumber() == project.getSettings()['n_iterations']-1):
import DataContainer as DC
export = 1
## compile dict that maps peaklist keys to name of new peaklist
ccpn_spectra = [x for x in project.ccpn_data_sources.get(DC.DATA_SPECTRUM, ())]
names_dict = {}
template = 'ARIA2_NOE_Peaks_%s_run%s_it%d_%s'
for x in ccpn_spectra:
ccpn_id = x['peaks']['ccpn_id']
args = file_root, run_name, iteration.getNumber(), ccpn_id
names_dict[ccpn_id] = template % args
peak_list_keys = ccpnExport.getPeakAssignmentsFromAria2(ccpn_project, restraints, namesDict=names_dict,
aria_chain=aria_chain)
self.message('CCPN export: NOE assignments exported.')
if export:
try:
ccpn_project.saveModified()
except Exception, msg:
self.warning('CCPN export: Could not save project. Error message was: %s' % msg) #by AWSS, msg is an object not a string
## unload objects to free some memory
if structures:
for s in structures:
#by AWSS begin
# CCPN API changed the way how unload() is called
try:
s.unload()
except Exception, msg:
self.warning("Couldn't unload structure: %s" % msg)
#by AWSS end
if constraintSet:
#by AWSS begin
# CCPN API changed the way how unload() is called
try:
constraintSet.unload()
except Exception, msg:
self.warning("Couldn't unload constraintSet: %s" % msg)
#by AWSS end
self.message('CCPN project saved.')
def _dump_iteration(self, iteration, path):
import os, Iteration
filename = os.path.join(path, REPORT_SUMMARY)
pickler = Iteration.IterationTextPickler()
pickler.dump(iteration, filename)
## BARDIAUX 2.2
def _dump_rms(self, peaks, iteration_number):
import RmsReport
infra = self.getInfrastructure()
text_path = infra.get_iteration_path(iteration_number)
graphics_path = infra.get_iteration_graphics_path(iteration_number)
rp = RmsReport.RmsReport(peaks, iteration_number, text_path, graphics_path)
rp.go()
def dumpIteration(self, iteration):
import os
from Singleton import ProjectSingleton
project = ProjectSingleton()
check_type(iteration, 'Iteration')
infra = self.getInfrastructure()
path = infra.get_iteration_path(iteration.getNumber())
# BARDIAUX 2.2
# dump also CCPN DistanceRestraints
# all_peaks = iteration.getPeaks()
all_peaks = {}
all_peaks.update(iteration.getPeaks())
all_peaks.update(iteration.getDistanceRestraints())
order = all_peaks.keys()
order.sort()
peak_list = []
for spectrum in order:
peaks = all_peaks[spectrum]
## sort peak_list wrt ref_peak number
peaks.sort(lambda a, b, c = cmp: \
c(a.getReferencePeak().getNumber(), \
b.getReferencePeak().getNumber()))
peak_list += peaks
## BARDIAUX 25/04/05
# Dump updated spectrum
if is_type(spectrum, 'ConstraintList'):
# Not for CCPN ConstraintList
continue
spec_name = spectrum.getName()
spec_name = spec_name.replace(' ','_')
filename = os.path.join(path, spec_name + \
REPORT_UPDATED_SPECTRUM +'.xml')
s = project.getReporter()['spectra']
if s['iteration'] == 'all' or \
(s['iteration'] == 'last' and \
iteration.getNumber() == project.getSettings()['n_iterations']-1):
self._dump_spectra(peaks, spectrum, filename)
self._dump_peak_list(peak_list, path)
self._dump_iteration(iteration, path)
self._dump_ccpn(iteration, path)
## BARDIAUX 2.2
## RMS report
if peak_list and iteration.getNumber() > 0:
try:
self._dump_rms(peak_list, iteration.getNumber())
except:
self.message('Error during RMS analysis/graphics generation.')
self.__reports_written = 1
def reportsWritten(self):
"""
returns 1 if any report files have been written.
"""
return self.__reports_written
def _get_peak_sizes(self, peaks):
"""
depending on the calibrator-setting 'volume_or_intensity',
the method returns volumes or intensities
"""
ref_peaks = [p.getReferencePeak() for p in peaks]
if self.calibrator.getSettings()['volume_or_intensity'] == 'volume':
peak_sizes = [p.getVolume()[0] for p in ref_peaks]
else:
peak_sizes = [p.getIntensity()[0] for p in ref_peaks]
return peak_sizes
# Bardiaux rMat
def _get_peak_theoric_volumes(self, peaks):
"""
get the theoric volume from the spin diffusin correction.
"""
return array([p.getTheoricVolume() for p in peaks])
# BARDIAUX 2.2
# get settings of DistanceRestraint sourec list
def __getListSource(self, p):
return p.getReferencePeak().getSpectrum().getListSource()
def setModelIntensities(self, restraints, ensemble, calibration_factor):
from Datum import Datum
f = self.model.calculatePeaksize
model_peak_sizes = array([f(r, ensemble) for r in restraints])
calculated_peak_sizes = model_peak_sizes * calibration_factor
for i in range(len(restraints)):
d = Datum(calculated_peak_sizes[i], None)
restraints[i].analysis.setCalculatedPeaksize(d)
def calculateBounds(self, factor, peaks, bound_corrected = None, ensemble = None):
"""
calculate lower- and upper bounds for every peak using
the calibration 'factor'. values are stored.
'bound_corrected': list of restraints which are classified
as correct after bound-modification.
"""
# BARDIAUX 2.2
# ConstraintList can bypass the calibration
if is_type(peaks[0], 'DistanceRestraint'):
calibrate = self.__getListSource(peaks[0])['calibrate']
if calibrate == NO or \
(self.findFirstIteration() == 0 and calibrate == 'all_iterations_except_first'):
return
factor = power(factor, 1./6)
peak_sizes = self._get_peak_sizes(peaks)
# Malliavin/Bardiaux rMat
cs = self.calibrator.getSettings()
if cs['relaxation_matrix'] == YES and ensemble is not None:
## from NOEModel import ISPA
## ispa = ISPA()
## f = ispa.calculatePeaksize
## ispa_peak_sizes = array([f(r, ensemble) for r in peaks])
ispa_peak_sizes = array([p.getIspa() for p in peaks])
peak_theoric_vol = self._get_peak_theoric_volumes(peaks)
ratio = ispa_peak_sizes / peak_theoric_vol
distances = factor * power(peak_sizes * ratio, -1. / 6)
else:
distances = factor * power(peak_sizes, -1. / 6)
## TODO: hard-coded 0.125
if cs['error_estimator'] == 'intensity':
errors = 0.125 * power((factor * power(peak_sizes, -1. / 6)), 2.)
else:
errors = 0.125 * power(distances, 2.)
## # BARDIAUX : To be implemented
## if cs['error_estimator'] == 'no_errors':
## errors = [0.] * len(distances)
# errors = [0.] * len(distances)
## distances = factor * power(peak_sizes, -1. / 6)
## ## TODO: hard-coded 0.125
## errors = 0.125 * power(distances, 2.)
## lower bounds are >= 0.
lower_bounds = clip(distances - errors, 0., 1.e10)
upper_bounds = distances + errors
for i in range(len(peaks)):
peak = peaks[i]
peak.setDistance(distances[i])
peak.setLowerBound(lower_bounds[i])
peak.setUpperBound(upper_bounds[i])
## Set new (fixed) bounds for bound-corrected restraints
if bound_corrected:
va_settings = self.violation_analyser.getSettings()
if va_settings['lower_bound_correction']['enabled'] == YES:
new_bound = va_settings['lower_bound_correction']['value']
[r.setLowerBound(new_bound) for r in bound_corrected]
if va_settings['upper_bound_correction']['enabled'] == YES:
new_bound = va_settings['upper_bound_correction']['value']
[r.setUpperBound(new_bound) for r in bound_corrected]
def doDumboCalibration(self, peaks):
# BARDIAUX 2.2
# ConstraintList can bypass calibration
if is_type(peaks[0], 'DistanceRestraint'):
do_cal = self.__getListSource(peaks[0])['calibrate'] == 'all_iterations_except_first'
if do_cal:
return 1.
peak_sizes = self._get_peak_sizes(peaks)
## Assume that average distance of atoms
## causing an NOE is 3.0A
d_calib = 3.0
sum_noe_calc = len(peaks) * (d_calib ** -6)
factor = sum(peak_sizes) / sum_noe_calc
return factor
def doCalibration(self, restraints, ensemble, store_analysis = 0):
do_cal = 1
# BARDIAUX 2.2
# ConstraintList can bypass calibration
if is_type(restraints[0], 'DistanceRestraint'):
do_cal = self.__getListSource(restraints[0])['calibrate'] == NO
if do_cal:
return 1.
if ensemble is None:
factor = self.doDumboCalibration(restraints)
else:
f = self.calibrator.calculateEstimator
factor = f(restraints, ensemble, store_analysis = store_analysis)
## if factor is None, i.e. if no NOEs stronger than a
## certain cutoff were found, set the cutoff to zero
## and calibrate again.
if factor is None:
s = self.calibrator.getSettings()
d_cutoff = s['distance_cutoff']
s = 'Could not perform 1st calibration, since ' + \
'no distances less than %.1f A were found in the ' + \
'ensemble. Omitting distance-cutoff and ' + \
'calibrating again...'
self.message(s % d_cutoff)
factor = f(restraints, ensemble, \
store_analysis = store_analysis, use_cutoff = 0)
return factor
def finalize_engine(self, molecule):
from Singleton import ProjectSingleton
project = ProjectSingleton()
engine = project.getStructureEngine()
engine.prepare(self.getSettings(), molecule)
def doAnalysis(self, restraints, ensemble):
##
## 1st calibration
##
## Store results.
##
factor = self.doCalibration(restraints, ensemble, store_analysis = 1)
##
## Violation analysis
##
## Store results:
##
## - average violation distance
## - whether a restraint is violated
## - degree of violation etc.
##
violated, non_violated, new_non_violated = \
self.doViolationAnalysis(restraints, ensemble,
store_analysis = 1)
## If possible, enlarge set of non-violated restraints
if new_non_violated:
non_violated += new_non_violated
##
## 2nd calibration
##
if non_violated:
factor = self.doCalibration(non_violated, ensemble)
## Calculate model peak-sizes
self.setModelIntensities(restraints, ensemble, factor)
def export_water_refined_structures(self, iteration):
if self.getSettings()['water_refinement']['enabled'] == NO:
return
infra = self.getInfrastructure()
path = infra.get_refinement_path()
self._dump_ccpn(iteration, path, is_water_refinement=True)
def writeReports(self, iteration):
t = clock()
s = 'Performing analysis on calculated structures...'
self.message(s)
ensemble = iteration.getStructureEnsemble()
# BARDIAUX 2.2
# add CCPN DistanceRestraints
all_restraints = {}
all_restraints.update(iteration.getPeaks())
all_restraints.update(iteration.getDistanceRestraints())
for spectrum, restraints in all_restraints.items():
#for spectrum, restraints in iteration.getPeaks().items():
## Set spectrum-specific settings
self._updateSpectrumSettings(spectrum)
self.doAnalysis(restraints, ensemble)
self.message('Analysis done.')
self.debug('Time: %ss' % str(clock() - t))
self.dumpIteration(iteration)
def readPDBFiles(self, filenames, molecule, float_files = None,
format = 'cns'):
"""
'filenames' is a dict or a tuple. In case of a dict,
keys are filenames, values are format-strings. In case of
a tuple, the argument 'format' is used as format-string.
Reads PDB-files 'filenames' and the respective
'float_files'. Returns a StructureEnsemble.
"""
check_type(filenames, TUPLE, DICT)
check_type(molecule, 'Molecule')
check_type(float_files, TUPLE, NONE)
check_string(format)
import StructureEnsemble as SE
if float_files == ():
float_files = None
if float_files is not None:
from FloatFile import FloatFile
parser = FloatFile()
swapped_atoms = [parser.parse(f) for f in float_files]
else:
swapped_atoms = None
## create structure-ensemble from given 'filenames'
se_settings = SE.StructureEnsembleSettings()
ensemble = SE.StructureEnsemble(se_settings)
ensemble.read(filenames, molecule, format, swapped_atoms)
return ensemble
def start(self, iteration, molecule):
"""
setup the first iteration. calculate distance-estimates
using either a template structure or by doing a dumbo-
calibration.
"""
check_type(molecule, 'Molecule')
check_type(iteration, 'Iteration')
## set-up first iteration
import os
from Singleton import ProjectSingleton
from DataContainer import DATA_TEMPLATE_STRUCTURE
t_iteration = clock()
## create seed-assignment
## if a template-structure has been specified,
## use that structure to create a better seed-assignment.
## otherwise we use an extended chain (in case of CNS,
## it is created by using generate_template.inp
infra = self.getInfrastructure()
project = ProjectSingleton()
templates = project.getData(DATA_TEMPLATE_STRUCTURE)
templates = [t for t in templates if t['enabled'] == YES]
## If template structures have been specified, we
## create an initial structure ensemble which will
## be used in the following steps to calculate
## an initial calibration factor etc.
if templates:
filenames = [t['filename'] for t in templates]
formats = [t['format'] for t in templates]
s = 'Template structure(s) specified (%s). Using template(s)' + \
' for calibration / seed-assignment generation.'
self.message(s % ', '.join(map(os.path.basename, filenames)))
## Compile local filenames (template PDB-files are
## read from local data-directory)
filenames = [infra.get_local_filename(f, DATA_TEMPLATE_STRUCTURE) \
for f in filenames]
## Build dict with filenames as keys and formats as
## values.
d = {}
for i in range(len(formats)):
d[filenames[i]] = formats[i]
ensemble = self.readPDBFiles(d, molecule)
ensemble.getSettings()['number_of_best_structures'] = 'all'
iteration.setStructureEnsemble(ensemble)
self.message('Initial iteration created.', verbose_level = VL_LOW)
self.debug('Time: %ss' % str(clock() -t_iteration))
## register initial iteration and start protocol
## self.iterations[iteration.getNumber()] = iteration
## Return last iteration
return self.run_protocol(molecule, iteration)
def go(self, iteration, molecule):
check_type(molecule, 'Molecule')
check_type(iteration, 'Iteration')
## Check which iterations have already been
## calculated, so that we can skip them.
first_iteration = self.findFirstIteration()
## If no iteration has been calculated yet, run the
## full protocol.
if first_iteration == 0:
iteration = self.start(iteration, molecule)
## If some iterations already exist, calculate the remaing ones.
elif first_iteration is not None:
t = clock()
s = 'Existing iterations found. Resuming with iteration ' + \
'%d.'
self.message(s % first_iteration)
iteration._setNumber(first_iteration - 1)
## Get filenames of structures which ought exist
## after parent-iteration has been successfully completed.
CFL = self.compileFileList
filenames = CFL(first_iteration - 1, structures = 1)
float_files = CFL(first_iteration - 1, float_files = 1)
it_settings = self.getIterationSettings(iteration.getNumber())
## Default PDB-format is 'cns'
ensemble = self.readPDBFiles(filenames, molecule,
float_files)
## Apply structure-ensemble settings of current
## iteration.
se_settings = ensemble.getSettings()
se_settings.update(it_settings)
ensemble.settingsChanged()
iteration.setStructureEnsemble(ensemble)
s = 'PDB-files / FLOAT-files for iteration %d read: %ss'
self.debug(s % (first_iteration, str(clock() - t)))
iteration = self.run_protocol(molecule, iteration)
else:
from Iteration import Iteration
self.message('Iterations are already complete.')
## TODO: this is a bit fiddled!
## Since the solvent-refinement step does not refer to any
## restraint, we simply create an empty iteration
## and set its number to the last iteration's.
infra = self.getInfrastructure()
last_iteration = infra.getSettings()['n_iterations'] - 1
iteration = Iteration(last_iteration)
## Start water refinement
self.doSolventRefinement(iteration, molecule)
## Export structures to CCPN
self.export_water_refined_structures(iteration)
## End of ARIA protocol.
self.done()
## return last iteration
return iteration
## BARDIAUX 2.2
def _checkSpectrumData(self, spectrum):
k = spectrum.getExperimentData().values()
return len([v for v in k if v <> 0.0]) == len(k)
def run_iteration(self, molecule, last_iteration, iteration,
is_first_iteration = 0):
"""
last_iteration is introduced, to provide time-shifted
pickling of restraint lists (i.e. while another
iteration is on its way)
"""
from time import sleep
check_type(molecule, 'Molecule')
check_type(last_iteration, 'Iteration')
## if all iterations shall be stored in memory ...
## self.iterations[iteration.getNumber()] = iteration
## apply settings of target iteration to all sub-modules
self._updateIterationSettings(iteration)
## get structure ensemble calculated in last iteration
ensemble = last_iteration.getStructureEnsemble()
## We now analyse the structure ensemble of the last
## iteration to generate a new set of restraints which
## will then be used to calculate a new ensemble of
## structures for the target iteration.
spectra = last_iteration.getPeaks().keys()
spectra.sort()
# BARDIAUX 2.2
# add CCPN DistanceRestraints
constraint_lists = last_iteration.getDistanceRestraints().keys()
constraint_lists.sort()
# Malliavin/Bardiaux rMat
cs = self.calibrator.getSettings()
for s in spectra:
if cs['relaxation_matrix'] <> NO and not self._checkSpectrumData(s):
self.warning("Missing experimental data for spectrum \"%s\". " % s.getName())
self.warning("Spin-diffusion correction will be disabled.\n")
cs['relaxation_matrix'] = NO
continue
if cs['relaxation_matrix'] == YES and ensemble is not None:
from NOEModel import SpinDiffusionCorrection
self.model = SpinDiffusionCorrection()
self.calibrator.setModel(self.model)
self.model.prepare(molecule, ensemble)
# initialize matrix
for spectrum in spectra:
self.model.setIntensityMatrix(spectrum)
spectra += constraint_lists
for spectrum in spectra:
self.message('Calibrating spectrum "%s"...' % spectrum.getName())
## get peaks stored in current iteration
## BARDIAUX 2.2 ConstraintList
if is_type(spectrum, 'ConstraintList'):
peaks = iteration.getDistanceRestraints(spectrum)
else:
peaks = iteration.getPeaks(spectrum)
##
## Calibration
##
## If we do not have a seed structure ensemble, we
## perform a very simple calibration
##
## set spectrum-specific parametes for Calibrator
self._updateSpectrumSettings(spectrum)
if ensemble is None:
factor = self.doCalibration(peaks, None)
new_non_violated = None
else:
## Calculate initial calibraton factor.
factor = self.doCalibration(peaks, ensemble)
## Calculate upper/lower bounds for restraints of current
## iteration.
t = clock()
self.calculateBounds(factor, peaks, ensemble = ensemble)
s = '1st calibration and calculation of new ' + \
'distance-bounds done (calibration factor: %e)'
self.message(s % factor, verbose_level = VL_LOW)
self.debug('Time: %ss' % str(clock() - t))
##
## Violation Analysis
##
## Assess every restraint regarding its degree of violation.
##
## Violated restraints will be disabled for the
## current iteration and thus will not be used during
## structure calculation.
##
t = clock()
violated, non_violated, new_non_violated = \
self.doViolationAnalysis(peaks, ensemble)
## Augment set of non-violated restraints
if new_non_violated:
non_violated += new_non_violated
n = len(peaks)
n_viol = len(violated)
p_viol = n_viol * 100. / n
s = 'Violation analysis done: %d / %d restraints ' + \
'(%.1f %%) violated.'
self.message(s % (n_viol, n, p_viol))
if new_non_violated:
s = 'Number of valid restraints has been increased ' + \
'by %d (%.1f%%) after applying a bound-correction.'
p_new = len(new_non_violated) * 100. / n
self.message(s % (len(new_non_violated), p_new))
self.debug('Time: %ss' % str(clock()-t))
##
## 2nd calibration - wrt to non-violated restraints.
## If no restraints have been violated, we use the
## 1st calibration factor.
## Again, we do not store results.
##
if non_violated:
factor = self.doCalibration(non_violated, ensemble)
##
## Activate restraints explicitly.
## We consider a restraint as active, if it has
## not been violated or if its reference cross-peak is
## 'reliable'.
##
for r in peaks:
if not r.analysis.isViolated() or \
r.getReferencePeak().isReliable():
r.isActive(1)
else:
r.isActive(0)
## Store final calibration factor for current iteration.
iteration.setCalibrationFactor(spectrum, factor)
## Calculate upper/lower bounds for restraint-list
## used in the current iteration. I.e. these bounds
## will be used to calculated the structures.
t = clock()
self.calculateBounds(factor, peaks, new_non_violated, ensemble = ensemble)
s = 'Final calibration and calculation of new distance-bounds' + \
' done (calibration factor: %e).' % factor
self.message(s, verbose_level = VL_LOW)
self.debug('Time: %ss' % str(clock() - t))
##
## Partial assignment for restraint-list used in
## current iteration, i.e. for all restraints:
## Calculate weight for every contribution
## and (depends possibly on partial analyser
## settings) throw away 'unlikely' contributions.
##
## If we do not have an ensemble, all contributions
## are activated.
##
t = clock()
# BARDIAUX 2.2
# if we are using the NA, we do not filter the contributions yet
# if the spectrum is a ConstraintList, check if needs to be filtered
NA = self.network_anchoring
ns = NA.getSettings()
do_filter = ns['enabled'] == NO
if is_type(spectrum, 'ConstraintList'):
do_filter = do_filter or \
spectrum.getListSource()['run_network_anchoring'] == NO
do_filter = do_filter and \
spectrum.getListSource()['filter_contributions'] == YES
self.contribution_assigner.assign(peaks, ensemble, filter_contributions = do_filter)
self.message('Partial assignment done.', verbose_level = VL_LOW)
self.debug('Time: %ss' % str(clock() - t))
# BARDIAUX 2.2 NA
NA = self.network_anchoring
ns = NA.getSettings()
if ns['enabled'] == YES:
NA.run(iteration)
# filter contributions
for spectrum in spectra:
# CCPN ConstraintList may not be filtered
if is_type(spectrum, 'ConstraintList'):
if spectrum.getListSource()['filter_contributions'] == NO:
continue
peaks = iteration.getDistanceRestraints(spectrum)
else:
peaks = iteration.getPeaks(spectrum)
self.contribution_assigner.filter_contributions(peaks)
##
## Merge spectra
##
self.mergePeakLists(iteration)
## to known when the structure generation has been completed.
self.__done = 0
## Display some cache statistics
if AriaBaseClass.cache and ensemble is not None:
cache = ensemble._StructureEnsemble__cache
hit_rate = float(cache['hit']) * 100. / cache['total']
self.debug('StructureEnsemble cache hit-rate: %.1f%%' % hit_rate)
##
## Run structure calculation
##
self.startStructureCalculation(iteration, molecule)
##
## Perform analysis wrt to last iteration
## Dump report files
##
if not is_first_iteration:
## apply settings of last iteration to sub-modules
self._updateIterationSettings(last_iteration)
self.writeReports(last_iteration)
## Restore old settings
self._updateIterationSettings(iteration)
## Wait for completion of structure generation.
self.message('Waiting for completion of structure calculation...')
try:
while not self.__done:
sleep(1.)
except KeyboardInterrupt:
self.__done = KEYBOARD_INTERRUPT
if self.__done in (MISSING_STRUCTURES, KEYBOARD_INTERRUPT):
## shut down job-manager
from Singleton import ProjectSingleton
project = ProjectSingleton()
engine = project.getStructureEngine()
engine.getJobScheduler().shutdown()
if self.__done == MISSING_STRUCTURES:
## BARDIAUX 2.2
# get list of missing structures
missing = engine.missingStructures()
msg = 'Structure calculation failed for structure %d.\n' #Please check your setup and the CNS output files for errors.'
err_msg = ''
for m in missing:
err_msg += msg % m
err_msg += 'Please check your setup and the CNS output files for errors.'
## err_msg = 'Structure calculation failed. ' + \
## 'Some structures are missing. ' + \
## 'Please check your setup and/or CNS output files.'
else:
err_msg = 'Interrupted by user.'
self.error(StandardError, err_msg)
## Return current, most recently calculated iteration
return iteration
def run_protocol(self, molecule, iteration):
check_type(molecule, 'Molecule')
check_type(iteration, 'Iteration')
from copy import copy
from Singleton import ProjectSingleton
infra = self.getInfrastructure()
project = ProjectSingleton()
engine = project.getStructureEngine()
msg = '---------------------- Iteration %d -----------------------\n'
is_first_it = 1
n_iterations = len(self.getSettings()['iteration_settings'])
t_main_protocol = clock()
n_first = iteration.getNumber()
for it_number in range(n_first, n_iterations - 1):
## create new iteration
t = clock()
target = copy(iteration)
target._setNumber(iteration.getNumber() + 1)
self.debug('New iteration created: %ss' % str(clock() - t))
self.message(msg % (target.getNumber()))
## Run calculation
t = clock()
iteration = self.run_iteration(molecule, iteration, target,
is_first_it)
## cleanup
if it_number < n_iterations - 1:
it_path = infra.get_iteration_path(it_number)
engine.cleanup(it_path)
is_first_it = 0
self.message('Iteration %s done.' % target.getNumber())
self.debug('Time: %ss' % str(clock() - t))
self.debug('Total time for all iterations: %ss' % \
str(clock() - t_main_protocol))
## Write reports for final iteration
self.writeReports(iteration)
self.writeFinalReports(iteration)
## Return last iteration
return iteration
def writeFinalReports(self, iteration):
check_type(iteration, 'Iteration')
import os, MolMol
from Singleton import ProjectSingleton
## write MOLMOL file-list "file.nam"
molmol_path = self._writePDBFileList(iteration)
if molmol_path is None:
msg = 'Could not write MOLMOL file-list.'
self.warning(msg)
return
## if that worked, attempt to write restraint-lists
project = ProjectSingleton()
s = project.getReporter()['molmol']
if s['enabled'] in (YES, GZIP):
## restraint files
restraints = iteration.getPeakList()
not_merged = [r for r in restraints if not r.isMerged()]
gzip = {YES: 0, GZIP: 1}[s['enabled']]
MolMol.write_noe_restraints(not_merged, molmol_path,
gzip = gzip)
self.message('MOLMOL lower and upper bound (.lol, .upl) ' + \
'files written.')
class ProtocolXMLPickler(XMLBasePickler):
order = ['floating_assignment', 'iteration',
'water_refinement']
def _xml_state(self, x):
e = XMLElement(tag_order = self.order)
s = x.getSettings()
e.iteration = tuple(s['iteration_settings'].values())
e.water_refinement = s['water_refinement']
e.floating_assignment = s['floating_assignment']
return e
def load_from_element(self, e):
from tools import as_tuple
s = ProtocolSettings()
it_settings = as_tuple(e.iteration)
[s.addIterationSettings(i) for i in it_settings]
s['water_refinement'] = e.water_refinement
s['floating_assignment'] = str(e.floating_assignment)
d = s['iteration_settings']
it_numbers = d.keys()
if min(it_numbers) <> 0:
self.warning('Iteration numbering does not start at 0, renumbering ...')
_min = min(it_numbers)
for number, it_settings in d.items():
new_number = number - _min
it_settings['number'] = new_number
d[new_number] = it_settings
del d[number]
p = Protocol(s)
return p
IterationSettings._xml_state = IterationSettingsXMLPickler()._xml_state
Protocol._xml_state = ProtocolXMLPickler()._xml_state
from ccpnmr.analysis.ConstraintBasic import makeNmrConstraintStore, makeStructureGeneration, getFixedResonance
from ccpnmr.analysis.ExperimentBasic import getOnebondDataDims
from ccpnmr.analysis.AssignmentBasic import assignAtomsToRes, assignResToDim
from ccpnmr.analysis.PeakBasic import pickPeak, setManualPeakIntensity
from importFromCcpn import getCcpnPeakList, getKeysFromString
def getAria2ObjectsFromPickle(fileName):
"""Descrn: Given an ARIA2 pickle file (and an istallation) get the ARIA2
data model objects.
Inputs: Word (pickle file name)
Output: List of Aria2 Restraint objects
"""
import cPickle
print "Reading Pickle file... "
fileName = os.path.expanduser(fileName)
if fileName[-3:].lower() == '.gz':
import gzip
file = gzip.open(fileName)
else:
file = open(fileName)
objects = None
while 1:
try:
objects= cPickle.load(file)
except EOFError:
break
file.close()
print " ...done."
return objects
def importAria2RunData(dirName, constraintSet=None, project=None, loadStructures=False):
"""Descrn: Import ARIA2 structures plus restraints and violations from a Pickle file
given an iteration director of an ARIA2 run.
Inputs: Word (run directory name), Nmr.NmrConstraintStore, Implementation.Project
Output: NmrConstraint.ConstraintList, NmrConstraint.ViolationList
"""
assert constraintSet or project
import aria
if constraintSet:
nmrProject = constraintSet.nmrConstraintStore.nmrProject
else:
nmrProject = project.currentNmrProject
ariaPath = os.environ.get(aria.ARIA_ENV)
if not ariaPath:
print 'Warning','Could not find ARIA2 environmant variable to locate the ARIA installation'
return
ariaPath = os.path.join(ariaPath, 'src/py')
sys.path.insert(0, ariaPath)
print "Getting constraint set"
constraintSet = makeNmrConstraintStore(nmrProject) # This function is present locally
pickleFile = None
fileNames = os.listdir(dirName)
for fileName in fileNames:
if fileName[-7:] == '.pickle':
pickleFile = fileName
break
elif fileName[-10:] == '.pickle.gz':
pickleFile = fileName
break
if not pickleFile:
print 'Failure','Found no .pickle file in directory %s' % dirName
return
print "Loading PICKLE file %s" % pickleFile
pickleFile = os.path.join(dirName,pickleFile)
objects = getAria2ObjectsFromPickle(pickleFile)
if not objects:
print 'Failure','Unable to find any Aria2 restraint objects'
sys.path.remove(ariaPath)
return
print "Getting chain"
# BARDIAUX : Update return list of chains
chains = getChainsFromAria2(objects, project)
if not chain:
print 'Failure','No molecular system chain found for ARIA2 data'
return
structures = None
if loadStructures:
print "loading structures"
# BARDIAUX : give first chain of chains list
structures = getStructuresFromAria2Dir(dirName, chains[0])
if not structures:
print 'Failure','Unable to load any structures from Aria2 iteration directory'
return
structureGeneration = makeStructureGeneration(structures, constraintSet)
structureGeneration.name = 'ARIA2 import'
sys.path.remove(ariaPath)
print "Interpreting ARIA2 restraints"
# BARDIAUX : Update needs chainS
return getConstraintsFromAria2(objects, chains, constraintSet, structures=structures)
def getChainFromAria2(objects, project, aria_chain=None):
"""Descrn: Get matching (or new) Ccp.MolSystem Chain from Aria2 restraint objects
Inputs: List of Aria2 restreaint objects, Implementation.Project
Output: MolSystem.chain
"""
from ccpnmr.analysis.MoleculeBasic import findMatchingChain
try:
from memops.gui.DataEntry import askString
except:
from memops.universal.DataEntry import askString
# TBD: deal with multiple chains
# BARDIAUX: Done. See getChainsFromAria2
if aria_chain is None:
ariaChain = objects[0].getReferencePeak().getProton1Assignments()[0].getAtoms()[0].residue.chain
else:
ariaChain = aria_chain
sequence = [r.type for r in ariaChain.getResidues()]
chains = []
#by AWSS begin
# issue: ccpn project with several molecular systems, user want use 3rd molSystem.
# molSystem should be the one user defined in the aria project.
# My solution:
molSystemName = aria_chain._Chain__molecule._Molecule__name
molSystem = project.findFirstMolSystem(code=molSystemName)
chain0 = findMatchingChain(molSystem, sequence, doWarning=True)[0] #by AWSS
if chain0:
chains.append(chain0)
## it was:
#for molSystem in project.molSystems:
# chain0 = findMatchingChain(molSystem, sequence, doWarning=True)[0] #by AWSS
# if chain0:
# chains.append(chain0)
# by AWSS end
chain = None
if len(chains) == 1:
chain = chains[0]
elif len(chains) > 1:
texts = ['%s:%s' % (ch.molSystem.code, ch.code) for ch in chains]
while chain is None:
response = askString('Enter Mol System and chain codes (Available: %s)' % ' '.join(texts),texts[0])
if not response:
continue
codes = response.split(':')
if len(codes) != 2:
continue
molSystem = project.findFirstMolSystem(code=codes[0])
if not molSystem:
continue
chain = molSystem.findFirstChain(code=codes[1])
# TBD: consider making molecule, chain and molSystem
return chain
# BARDIAUX
def getChainsFromAria2(objects, project, aria_chain=None):
"""Descrn: Get matching (or new) Ccp.MolSystem Chains from Aria2 restraint objects
Inputs: List of Aria2 restreaint objects, Implementation.Project
Output: List of MolSystem.chain
"""
from ccpnmr.analysis.MoleculeBasic import findMatchingChain
if aria_chain is not None:
ariaChain = aria_chain
else:
ariaChain = objects[0].getReferencePeak().getProton1Assignments()[0].getAtoms()[0].residue.chain
sequence = [r.type for r in ariaChain.getResidues()]
chains = []
#by AWSS begin
# issue: ccpn project with several molecular systems, user want use 3rd molSystem.
# molSystem should be the one user defined in the aria project.
# My solution:
molSystemName = aria_chain._Chain__molecule._Molecule__name
molSystem = project.findFirstMolSystem(code=molSystemName)
for molChain in molSystem.chains:
chain0 = findMatchingChain(molSystem, sequence, excludeChains=chains, doWarning=True)[0] #by AWSS
if chain0:
chains.append(chain0)
## it was:
#for molSystem in project.molSystems:
# for molChain in molSystem.chains:
# chain0 = findMatchingChain(molSystem, sequence, excludeChains=chains, doWarning=True)[0] #by AWSS
# if chain0:
# chains.append(chain0)
# by AWSS end
chain = None
if len(chains) == 1:
chain = chains[0]
return chains
def getStructuresFromAria2Dir(dirName, chain):
"""Descrn: Load structures of the input chain, from the PDB style files in an Aria2 iteration directory
Inputs: Word - Aria2 run directory path. MolSystem.Chain
Output: List of MolSystem.MolStructures
"""
from ccpnmr.analysis.StructureBasic import getStructureFromFile
import os
structures = []
fileNames = os.listdir(dirName)
for fileName in fileNames:
if fileName[-4:] == '.pdb':
fileName = os.path.join(dirName,fileName)
fileName = os.path.expanduser(fileName)
structures.append( getStructureFromFile(chain.molSystem, fileName, doWarnings=False) )
return structures
def getConstraintFromAria2(ariaRestraint, constraintList):
"""Descrn: Make a distance constraint from an Aria2 restraint object
Inputs: Aria2 restraint object, NmrConstraint.DistanceConstraintList
Output: NmrRestraints.DistanceConstraint
"""
nmrProject = constraintList.topObject.nmrProject
if not hasattr(constraintList, 'peakLookup'):
constraintList.peakLookup = {}
weight = ariaRestraint.getWeight()
upperB = ariaRestraint.getUpperBound()
lowerB = ariaRestraint.getLowerBound()
rfPeak = ariaRestraint.getReferencePeak()
active = ariaRestraint.isActive()
target = ariaRestraint.getDistance()
contrb = ariaRestraint.getContributions()
# TBD deal with weight - exclude some constraints given weight threshold
peak = None
peaksDict = constraintList.peakLookup
spectrumName = rfPeak.getSpectrum().getDataSource().get('ccpn_id')
# BARDIAUX 2.2
from aria import is_type
# create the CCPN DistanceConstraint first
origData = rfPeak.getIntensity().getValue()
if not origData:
origData = rfPeak.getVolume().getValue()
# BARDIAUX 2.2 backCalcVolume
backCalcVolume = ariaRestraint.analysis.getCalculatedPeaksize()[0]
constraint = constraintList.newDistanceConstraint(weight=weight, error=abs(upperB-lowerB),
origData=origData, targetValue=target,
upperLimit=upperB, lowerLimit=lowerB,
backCalcVolume = backCalcVolume)
# BARDIAUX 2.2 now deals with ConstraintPeakContrib #typo fixed by AWSS
if spectrumName:
ccpnNames = getKeysFromString(spectrumName) # TJS
## BARDIAUX 2.2
if is_type(rfPeak.getSpectrum(), 'ConstraintList'):
# Constraint comes from a CCPN constraintList, not a spectrum
# first, retrieve the original constraintList
nmrProjectKey, storeSerial, serial = ccpnNames
nmrProject = constraintList.root.findFirstNmrProject(name=nmrProjectKey)
if not nmrProject:
raise ValueError('No NMR project with name %s' % nmrProjectKey)
constraintStore = nmrProject.findFirstNmrConstraintStore(serial=storeSerial)
if not constraintStore:
raise ValueError('No NMR constraint store with serial "%d" in CCPN project' % storeSerial)
else:
origConstraintList = constraintStore.findFirstConstraintList(serial=serial)
if origConstraintList is None:
raise ValueError('No constraint list with serial "%d" in store "%d"' % (serial,storeSerial))
# second, find the original DistanceConstraint
constraintSerial = rfPeak.getNumber()
origConstraint = origConstraintList.findFirstConstraint(serial=constraintSerial)
# track the original constraint
constraint.details = 'Orig Id %d' % constraintSerial
# now, get the original ConstraintPeakContribs and transfer them to the
# new DistanceConstraint
origConstraintPeakContribs = origConstraint.peakContribs
if origConstraintPeakContribs:
for origConstraintPeakContrib in origConstraintPeakContribs:
experimentSerial = origConstraintPeakContrib.experimentSerial
dataSourceSerial = origConstraintPeakContrib.dataSourceSerial
peakListSerial = origConstraintPeakContrib.peakListSerial
peakSerial = origConstraintPeakContrib.peakSerial
peakContrib = constraint.newConstraintPeakContrib(experimentSerial= experimentSerial,
dataSourceSerial= dataSourceSerial,
peakListSerial= peakListSerial,
peakSerial=peakSerial)
else:
# BARDIAUX 2.2
# Restraints from a Spectrum
if len(ccpnNames) == 4:
#Fix TJS: NmrProject is part of the key
#keys are serial num for experiment and spectrum not name
nName, eSerial, sSerial, plSerial = ccpnNames
experiment = nmrProject.findFirstExperiment(serial=eSerial)
if experiment:
spectrum = experiment.findFirstDataSource(serial=sSerial)
if spectrum:
peakList = spectrum.findFirstPeakList(serial=plSerial)
if peakList:
if peaksDict.get(peakList) is None:
peaksDict[peakList] = {}
for peak in peakList.peaks:
peaksDict[peakList][peak.serial] = peak
peak = peaksDict[peakList].get(rfPeak.getNumber())
if peak:
peakContrib = constraint.newConstraintPeakContrib(experimentSerial=experiment.serial,
dataSourceSerial=spectrum.serial,
peakListSerial=peakList.serial,
peakSerial=peak.serial)
return constraint
def getViolationFromAria2(ariaViolation, constraint, violationList):
"""Descrn: Make a constraint's violation object in the input violation list given an Aria2
violation analysis object.
Inputs: Aria2 violation analysis object, NmrConstraint.DistanceConstraint, NmrConstraint.ViolationList
Output: NmrConstraint.Violation
"""
violation = None
if ariaViolation.isViolated():
calcDist = ariaViolation.getAverageDistance().getValue()
calcError = ariaViolation.getAverageDistance().getError()
fracViols = max(0.0, min(ariaViolation.getDegreeOfViolation(), 1.0))
violValue = ariaViolation.getUpperBoundViolation().getValue()
violation = violationList.newViolation(violation=violValue,calcValue=calcDist,
calcValueError=calcError,constraint=constraint,
fractionViolated=fracViols)
return violation
def getAria2SpinSystems(ariaRestraint):
"""Descrn: Get the pairs of ARIA spin systems that correspond to the contributions
to an ARIA restraint
Inputs: Aria2 restraint object
Output: List of 2-List of ARIA SpinSystem objects
"""
spinSystemPairs = []
for contribution in ariaRestraint.getContributions():
ariaSpinSystems = contribution.getSpinSystems()
if len(ariaSpinSystems) == 2:
spinSystemPairs.append(list(ariaSpinSystems))
return spinSystemPairs
def getAria2AtomSetPairs(ariaRestraint):
"""Descrn: Get a list of pairs of aria atom objects from an Aria2 restraint
Which correspond to the atom sets of a DistanceConstraint.
Inputs: Aria2 restraint object
Output: List of 2-List of Aria2 atom objects
"""
ariaAtomPairs = []
for contribution in ariaRestraint.getContributions():
if contribution.getWeight() == 0:
continue
ariaAtomSetPair = []
for ariaSpinSystem in contribution.getSpinSystems():
ariaSpinSystemAtoms = ariaSpinSystem.getAtoms()
if len(ariaSpinSystemAtoms) != 2:
ariaAtomSetPair.append([ariaSpinSystemAtoms[0],])
else:
for ariaSpinSystemAtom in ariaSpinSystemAtoms:
if ariaSpinSystemAtom.getName()[-1] == '1':
ariaAtomSetPair.append([ariaSpinSystemAtom,])
break
else:
ariaAtomSetPair.append(ariaSpinSystemAtoms)
if len(ariaAtomSetPair) == 2:
ariaAtomPairs.append(ariaAtomSetPair)
return ariaAtomPairs
def getConstraintsFromAria2(objects, chains, constraintSet, structures=None, details='ARIA2 import'):
"""Descrn: Make a constraint and violation list from Aria2 restraint objects in a constraint set.
Constraint asignments will be made relative to the input chain.
The violation analysis will be linked to the input structures.
Inputs: List of Aria2 restraint objects, List of MolSystem.chain,
MolSystem.MolStructures, Nmr.NmrConstraintStore, String
Output: NmrRestraints.DistanceConstraintList, NmrConstraint.ViolationList
"""
# BARDIAUX: Update for multiple chains.
resonanceDict = {}
constraintList = constraintSet.newDistanceConstraintList(details=details)
if structures:
violationList = constraintSet.newViolationList(details=details, molStructures=structures)
else:
violationList = constraintSet.newViolationList(details=details)
# TJS collect inactive restraints for inspection in CCPN
# Only make a DistanceConstraintList for these if required
rejectConstraintList = None
for ariaRestraint in objects:
ariaAtomPairs = getAria2AtomSetPairs(ariaRestraint)
if not ariaAtomPairs:
continue
# TJS adjust to collect rejects
if ariaRestraint.isActive():
constraint = getConstraintFromAria2(ariaRestraint, constraintList)
else:
# TJS we want the restraints unmerged if one is active
# (because both have valid CCPN peaks)
equivPeak = ariaRestraint.getEquivalentPeak()
if ariaRestraint.isMerged() and equivPeak and equivPeak.isActive():
constraint = getConstraintFromAria2(ariaRestraint, constraintList)
else:
if not rejectConstraintList:
rejectConstraintList = constraintSet.newDistanceConstraintList(details=details + ' *REJECTS*')
constraint = getConstraintFromAria2(ariaRestraint, rejectConstraintList)
if not constraint:
continue
violation = getViolationFromAria2(ariaRestraint.analysis, constraint, violationList)
doneItems = {}
for ariaAtomSet1, ariaAtomSet2 in ariaAtomPairs:
for ariaAtom1 in ariaAtomSet1:
chain1 = getChainFromAria2Atom(ariaAtom1, chains)
fixedResonance1 = resonanceDict.get(ariaAtom1)
if fixedResonance1 is None:
fixedResonance1 = getFixedResonanceFromAria2Atom(ariaAtom1, chain1, constraintSet)
for ariaAtom2 in ariaAtomSet2:
chain2 = getChainFromAria2Atom(ariaAtom2, chains)
fixedResonance2 = resonanceDict.get(ariaAtom2)
if fixedResonance2 is None:
fixedResonance2 = getFixedResonanceFromAria2Atom(ariaAtom2, chain2, constraintSet)
# BARDIAUX 2.2
# check if fixedResonance1 too
if fixedResonance2 and fixedResonance1 and ( fixedResonance2 is not fixedResonance1 ):
serials = [fixedResonance1.serial,fixedResonance2.serial]
serials.sort()
serials = tuple(serials)
if doneItems.get(serials) is None:
item = constraint.newDistanceConstraintItem(resonances=[fixedResonance1,fixedResonance2])
doneItems[serials] = True
if not constraint.items:
constraint.delete()
#print "Constraint %d has no items" % constraint.serial
#print "Aria restraint atoms", ariaAtomPairs
if hasattr(constraintList, 'peakLookup'):
serials = {}
for peakList in constraintList.peakLookup.keys():
serials[ peakList.dataSource.experiment.serial ] = None
constraintList.setExperimentSerials(serials.keys())
del constraintList.peakLookup
for chain in chains:
if hasattr(chain,'residueLookup'):
del chain.residueLookup
return constraintList, rejectConstraintList, violationList
def getPeakAssignmentsFromAria2(project, ariaRestraints, namesDict=None,
aria_chain=None):
"""Descrn: Assign a CCPN peak list (new or existing) according to ARIA restraint assignments
Inputs: Implementation.Project, List of ARIA Restraints, Boolean
Output: ccp.nmr.Nmr.PeakList
"""
# BARDIAUX: Update for multiple chains
chains = getChainsFromAria2(ariaRestraints, project, aria_chain=aria_chain)
# Really need to deal with multiple chains
ariaDimDict = {}
peakListDict = {}
peaksDict = {}
bondedDimsDict = {}
for ariaRestraint in ariaRestraints:
ariaAtomPairs = getAria2AtomSetPairs(ariaRestraint)
ariaSpinSystemsPairs = getAria2SpinSystems(ariaRestraint)
resonanceDimPairs = []
refPeak = ariaRestraint.getReferencePeak()
for i in range(len(ariaAtomPairs)):
atomPair = ariaAtomPairs[i]
spinSystemPair = ariaSpinSystemsPairs[i]
atomPair = [a[0] for a in atomPair]
# BARDIAUX: Update for multiple chains
chain0, chain1 = getChainFromAria2Atom(atomPair[0], chains), getChainFromAria2Atom(atomPair[1], chains)
resonance1 = getResonanceFromAria2Atom(chain0, atomPair[0])
resonance2 = getResonanceFromAria2Atom(chain1, atomPair[1])
dim1 = (refPeak.getDimension(spinSystemPair[0]) -1) *2
dim2 = (refPeak.getDimension(spinSystemPair[1]) -1) *2
if dim1 == dim2:
print "WARNING: Problem resolving dimensions of peak assignments for CCPN"
if dim2 == 2:
dim2 = 0
else:
dim2 = 2
resonanceDimPairs.append((dim1,resonance1,dim2,resonance2))
peakListKeys = tuple(getKeysFromString(refPeak.getSpectrum().getDataSource().get('ccpn_id')))
if not peakListKeys:
print "Aria Restraint ", ariaRestraint, " has no associated CCPN spectrum/peak list"
continue
peakList = peakListDict.get(peakListKeys)
if peakList is None:
peakList = getCcpnPeakList(project, peakListKeys)
spectrum = peakList.dataSource
if namesDict:
new_name = namesDict[refPeak.getSpectrum().getDataSource().get('ccpn_id')]
peakList = spectrum.newPeakList(name=new_name)
peakList.details = new_name
peakListDict[peakListKeys] = peakList
else:
spectrum = peakList.dataSource
if peaksDict.get(peakList) is None:
peaksDict[peakList] = {}
for peak0 in peakList.peaks:
peaksDict[peakList][peak.serial] = peak0
onebondDims = bondedDimsDict.get(spectrum)
if onebondDims is None:
onebondDims = {}
for dataDim1, dataDim2 in getOnebondDataDims(spectrum):
onebondDims[dataDim1.dim] = dataDim2.dim
onebondDims[dataDim2.dim] = dataDim1.dim
bondedDimsDict[spectrum] = onebondDims
ppmH1 = refPeak.getProton1ChemicalShift().getValue()
ppmX1 = refPeak.getHetero1ChemicalShift().getValue()
ppmH2 = refPeak.getProton2ChemicalShift().getValue()
ppmX2 = refPeak.getHetero2ChemicalShift().getValue()
ppms = [ppmH1, ppmX1, ppmH2, ppmX2]
ariaDims = ariaDimDict.get(spectrum)
if not ariaDims:
ariaDims = [] #[0,1,2]
dataDims = spectrum.sortedDataDims()
if len(dataDims) == 3:
for dataDim in dataDims:
expDimRef = dataDim.findFirstDataDimRef().expDimRef
if '1H' in expDimRef.isotopeCodes: # 0 or 2
if onebondDims.get(dataDim.dim):
if ppmX1 is None:
ariaDims.append(2)
else:
ariaDims.append(0)
else:
if ppmX1 is None:
ariaDims.append(0)
else:
ariaDims.append(2)
else: # 1 or 3
if ppmX1 is None:
ariaDims.append(3)
else:
ariaDims.append(1)
else:
transfer = spectrum.experiment.findFirstExpTransfer(transferType='NOESY')
for dataDim in dataDims:
expDimRefs = [dataDimRef.expDimRef for dataDimRef in dataDim.dataDimRefs]
i = 0
for expDimRef in transfer.sortedExpDimRefs():
if expDimRef in expDimRefs:
ariaDims.append(i)
boundDim = onebondDims.get(dataDim.dim)
if boundDim:
ariaDims.append(i+1)
i += 2
ariaDimDict[spectrum] = ariaDims
if namesDict:
# Fix TJS # # # # # # # # # # # # # # # #
position = [ppms[dim] for dim in ariaDims]
peak = pickPeak(peakList, position, unit='ppm', doFit=False)
setManualPeakIntensity(peak, refPeak.getVolume().getValue(), intensityType='volume')
else:
peak = peaksDict[peakList].get(refPeak.getNumber())
if not peak:
data = [spectrum.experiment.name,
spectrum.name,
peakList.serial,
peak.serial]
print 'Cannot find CCPN Peak %s %s %d number %d' % data
continue # TJS should not try to assign
# TJS mark rejected peaks - would not be used next time by default
if not ariaRestraint.isActive():
peak.figOfMerit = 0.0
peak.details = 'ARIA2 REJECT'
# TJS rearranges this to put in checks for isotope mismatch
for dim1, resonance1, dim2, resonance2 in resonanceDimPairs:
peakAssignments = []
ccpnDim = 0 # Tim's
for dim in ariaDims:
ccpnDim += 1 # Fix TJS # # # # # # # # # # # # # # # #
peakDim = peak.findFirstPeakDim(dim=ccpnDim)
if dim1 == dim:
peakAssignments.append((peakDim, resonance1))
dimX = onebondDims.get(ccpnDim)
if dimX:
peakDimX = peak.findFirstPeakDim(dim=dimX)
bound = resonance1.findFirstCovalentlyBound()
if bound:
peakAssignments.append((peakDimX, bound))
elif dim2 == dim:
peakAssignments.append((peakDim, resonance2))
dimX = onebondDims.get(ccpnDim)
if dimX:
peakDimX = peak.findFirstPeakDim(dim=dimX)
bound = resonance2.findFirstCovalentlyBound()
if bound:
peakAssignments.append((peakDimX, bound))
# TJS the actual isotope mismatch check
for peakDim, resonance in peakAssignments:
if resonance.isotopeCode not in peakDim.dataDimRef.expDimRef.isotopeCodes:
print 'WARNING: Attempt to assign %s resonance to peak dim %s' % (resonance.isotopeCode,peakDim)
break
else:
# TJS do the peak assignment if all OK
for peakDim, resonance in peakAssignments:
if resonance.isotopeCode == '1H':
tolerance=0.5
else:
tolerance=2.0
assignResToDim(peakDim, resonance, tolerance=tolerance, doWarning=False)
return peakListDict.keys()
def getResonanceFromAria2Atom(chain, ariaAtom):
"""Descrn: Get the corresponding (or new) resonance for the input Aria2 atom
given the input chain
Inputs: MolSystem.Chain, Aria2 atom object
Output: Nmr.Resonance
"""
project = chain.root
if hasattr(project, 'ariaAtomCache'):
cache = project.ariaAtomCache
else:
project.ariaAtomCache = cache = {}
resonance = cache.get(ariaAtom)
if resonance:
return resonance
if hasattr(chain,'residueLookup'):
residueDict = chain.residueLookup
else:
# Dictionary for quick residue lookup
residueDict = {}
for residue in chain.residues:
residueDict[residue.seqCode] = residue
chain.residueLookup = residueDict
tlc = ariaAtom.getResidue().getType()
seqId = ariaAtom.getResidue().getNumber()
residue = residueDict.get(seqId)
if (not residue) or (residue.ccpCode != tlc):
print 'Warning','Cannot find residue %s %s' % (seqId,tlc)
return
atom = residue.findFirstAtom(name=ariaAtom.name)
if not atom:
print 'Warning','Cannot find atom %s %s %s' % (seqId,tlc,ariaAtom.name)
return
resonance = None
atomSet = atom.atomSet
if not atomSet:
print 'Missing atom set %s %s %s' % (seqId,tlc,atom.name)
# BARDIAUX if no atom set, return FixedAtomSet
atomSet = atom.findFirstFixedAtomSet()
if not atomSet:
return resonance
if not list(atomSet.resonanceSets):
return resonance
resonanceSets = list(atomSet.resonanceSets)
if resonanceSets:
for resonanceSet in resonanceSets:
numSets = len(resonanceSet.atomSets)
if numSets == 1:
resonance = resonanceSet.findFirstResonance()
if not resonance:
resonanceSet = resonanceSets[0]
atomSets = resonanceSet.sortedAtomSets()
#index = atomSets.index(atom.atomSet)
index = atomSets.index(atomSet)
resonances = resonanceSet.sortedResonances()
resonance = resonances[min(index,len(resonances)-1)]
else:
# Make new resonance if atom not assigned
resonance = chain.root.newResonance(isotopeCode='1H')
assignAtomsToRes([atomSet,],resonance)
cache[ariaAtom] = resonance
return resonance
def getFixedResonanceFromAria2Atom(ariaAtom, chain, constraintSet):
"""Descrn: Get the corresponding (or new) FixedResonance for the input Aria2 atom
given the input chain
Inputs: Aria2 atom object, MolSystem.Chain, Nmr.NmrConstraintHead
Output: NmrRestraints.FixedResonance
"""
resonance = getResonanceFromAria2Atom(chain, ariaAtom)
# BARDIAUX 2.2
# in case ccpnAtom has no atomSet
# and FixedResonance has no resonance
if not resonance:
return resonance
if resonance.className == 'FixedResonance':
if not resonance.resonance:
return copyFixedResonance(resonance, constraintSet)
return getFixedResonance(constraintSet,resonance)
# BARDIAUX
def getChainFromAria2Atom(ariaAtom, chains):
"""Descrn: Get the corresponding MolSystem.chain for the input Aria2 atom
Inputs: Aria2 atom object, List of MolSystem.Chain
Output: MolSystem.Chain
"""
from tools import string_to_segid
project = chains[0].root
if hasattr(project, 'ariaAtomCache'):
cache = project.ariaAtomCache
else:
project.ariaAtomCache = cache = {}
for chain in chains:
if ariaAtom.getSegid() <> string_to_segid(chain.code):
continue
if hasattr(chain,'residueLookup'):
residueDict = chain.residueLookup
else:
# Dictionary for quick residue lookup
residueDict = {}
for residue in chain.residues:
residueDict[residue.seqCode] = residue
chain.residueLookup = residueDict
tlc = ariaAtom.getResidue().getType()
seqId = ariaAtom.getResidue().getNumber()
residue = residueDict.get(seqId)
if (not residue) or (residue.ccpCode != tlc):
#print 'Warning','Cannot find residue %s %s' % (seqId,tlc)
continue
atom = residue.findFirstAtom(name=ariaAtom.name)
if not atom:
#print 'Warning','Cannot find atom %s %s %s' % (seqId,tlc,ariaAtom.name)
continue
else:
return chain
## Tim Stevens code
def copyFixedResonance(fixedResonance, nmrConstraintStore):
"""Descrn: Make an equivalent fixed resonance in the input constrsint store
Inputs: NmrConstraint.FixedResonance, Nmr.NmrConstraintStore
Output: NmrConstraint.FixedResonance
"""
# BARDIAUX 2.2
# this function is used in case of fixedResonance have no Nmr.Resonance
fixedResonance2 = None
if fixedResonance.resonanceSerial:
fixedResonance2 = nmrConstraintStore.findFirstFixedResonance(resonanceSerial=fixedResonance.resonanceSerial)
if fixedResonance2:
return fixedResonance2
serial = fixedResonance.resonanceSerial
name = fixedResonance.name
details = fixedResonance.details
isotope = fixedResonance.isotopeCode
resonanceSet = fixedResonance.resonanceSet
if resonanceSet:
atomSets2 = [getFixedAtomSet(nmrConstraintStore, as.atoms) for as in resonanceSet.atomSets]
resonanceSet2 = nmrConstraintStore.findFirstFixedResonanceSet(atomSets=atomSets2)
if resonanceSet2:
index = list(resonanceSet.resonances).index(fixedResonance)
fixedResonance2 = list(resonanceSet2.resonances)[index]
else:
fixedResonance2 = nmrConstraintStore.newFixedResonance(name=name,
details=details,
isotopeCode=isotope,
resonanceSerial=serial)
fixedResonances2 = [fixedResonance2]
for fixedResonance1 in resonanceSet.resonances:
if fixedResonance1 is not fixedResonance:
fixedResonance3 = nmrConstraintStore.newFixedResonance(name=fixedResonance1.name,
details=fixedResonance1.details,
isotopeCode=fixedResonance1.isotopeCode,
resonanceSerial=fixedResonance1.resonanceSerial)
fixedResonances2.append(fixedResonance3)
resonanceSet2 = nmrConstraintStore.newFixedResonanceSet(atomSets=atomSets2,resonances=fixedResonances2)
else:
fixedResonance2 = nmrConstraintStore.newFixedResonance(name=name,
details=details,
isotopeCode=isotopeCode,
resonanceSerial=resonanceSerial)
return fixedResonance2
def getFixedAtomSet(nmrConstraintStore, atoms):
"""Descrn: Finds or creates a fixed set of atoms that is used in an NMR constraint head object (equivalent to one
NmrConstraint file). Creating fixed atom sets allows assignments to change but old constraints to be preserved.
Inputs: Nmr.NmrConstraintStore, List of MolSystem.Atoms
Output: NmrConstraint.FixedAtomSet
"""
if not hasattr(nmrConstraintStore, 'quickAtomSets'):
nmrConstraintStore.quickAtomSets = {}
fixedAtomSet = nmrConstraintStore.quickAtomSets.get(atoms)
if not fixedAtomSet:
fixedAtomSet = nmrConstraintStore.findFirstFixedAtomSet(atoms=atoms)
if not fixedAtomSet:
if atoms[0].atomSet:
atomSet = atoms[0].atomSet
fixedAtomSet = nmrConstraintStore.newFixedAtomSet(atoms=atomSet.atoms,name=atomSet.name)
# BARDIAUX 2.2
# if atom have no atomSet
if not fixedAtomSet:
fixedAtomSet2 = atoms[0].findFirstFixedAtomSet(atoms=atoms)
if fixedAtomSet2:
fixedAtomSet = nmrConstraintStore.newFixedAtomSet(atoms=fixedAtomSet2.atoms,name=fixedAtomSet2.name)
nmrConstraintStore.quickAtomSets[atoms] = fixedAtomSet
return fixedAtomSet
from Molecule import Molecule
from Chain import TYPE_PROTEIN, TYPE_DNA, TYPE_RNA, TYPE_NONPOLYMER, Chain, ChainSettings
from Residue import Residue
from Atom import Atom, AtomSettings
from ccpn2top import setupEquivalentGroup
from Assignment import Assignment, ASSIGNMENT_TYPE_MANUAL
from NOESYSpectrum import NOESYSpectrum
from CrossPeak import CrossPeak
from Datum import ChemicalShift, Datum
from ChemicalShiftList import ChemicalShiftList
from ShiftAssignment import AVERAGING_METHOD_FAST, AVERAGING_METHOD_SLOW, AVERAGING_METHOD_NONE,\
ASSIGNMENT_METHOD_STEREO_SPECIFIC, ASSIGNMENT_METHOD_EQUIVALENT, \
ASSIGNMENT_METHOD_FLOATING, SpinSystem, ShiftAssignment
from tools import string_to_segid
from ccpnmr.format.general.Constants import ccpNmr_kw, isAriaInput_kw
from memops.api import Implementation
from memops.universal.Constants import True
from Settings import Settings
from xmlutils import XMLBasePickler
chainTypeMapping = {'protein' : TYPE_PROTEIN,
'DNA' : TYPE_DNA,
'RNA' : TYPE_RNA,
'nonpolymer': TYPE_NONPOLYMER}
def getKeysFromString(word, delimiter='|'):
"""Descrn: Get a list of CCPN data model keys given an object identifier string.
Inputs: String
Output: List of Keys (Words or Ints)
"""
items = word.split(delimiter)
keys = []
for item in items:
try:
key = int(item)
except:
key = item
keys.append(key)
return keys
def getObjectKeys(object):
"""Descrn: Get a list of CCPN data model keys for a CCPN data model object.
Inputs: CCPN data model object
Output: List of Keys (Words or Ints or a further list of keys)
"""
return object.getExpandedKey()
def getObjectKeyString(object, delimiter='|'):
"""Descrn: Make an object identifier string for a CCPN data model object
Inputs: CCPN data model object
Output: String
"""
keys = object.getExpandedKey()
for i in range(len(keys)):
key = keys[i]
keyType = type(key)
if keyType is type([]):
keys[i] = delimiter.join([str(k) for k in key])
elif keyType is not type(''):
keys[i] = str(key)
return delimiter.join(keys)
# BARDIAUX
def getObjectsKeyString(o):
s = []
for a in o:
s.append(getKeysFromString(getObjectKeyString(a)))
base = s[0]
for i,j in zip(*s):
if i <> j:
base += j
return "|".join(base)
def getShiftLists(project):
"""Descrn: Get a list of shift lists from a CCPN project
Inputs: Implementation.Project
Output: List of ccp.nmr.Nmr.ShiftLists
"""
nmrProject = project.currentNmrProject
shiftLists = []
for measurementList in nmrProject.sortedMeasurementLists():
if measurementList.className == 'ShiftList':
shiftLists.append(measurementList)
return shiftLists
def getNoesyPeakLists(project, molSystem=None):
"""Descrn: Get the NOE peak lists from a CCPN project. Can filter if appropriate to a given
molecular system if passed in.
Inputs: Implementation.Project, ccp.molecule.MolSystem.MolSystem
Output: List of ccp.nmr.Nmr.Peaks
"""
peakLists = []
for experiment in project.currentNmrProject.experiments:
if experiment.molSystems:
if molSystem and (molSystem not in experiment.molSystems):
continue
if experiment.refExperiment:
if experiment.refExperiment.nmrExpPrototype.name.find('NOESY') < 0:
continue
else:
transfer = experiment.findFirstExpTransfer(transferType='NOESY')
if not transfer:
continue
for spectrum in experiment.dataSources:
if (spectrum.dataType == 'processed') and (spectrum.numDim > 1):
isotopes = []
for dataDim in spectrum.dataDims:
for expDimRef in dataDim.expDim.expDimRefs:
if expDimRef.measurementType in ('shift','Shift'):
isotope = ','.join(expDimRef.isotopeCodes)
isotopes.append(isotope)
break
if isotopes.count('1H') > 1:
for peakList in spectrum.peakLists:
if peakList.findFirstPeak():
peakLists.append(peakList)
return peakLists
def getCcpnPeakList(ccpnProject, peakListKey):
"""Descrn: Fetch a CCPN peak list from a CCPN project using a list of object keys
Inputs: Implementation.Project, List of keys (Words or Ints)
Output: ccp.nmr.Nmr.PeakList
"""
nmrProjectName, exptSerial, specSerial, plSerial = peakListKey
nmrProject = ccpnProject.findFirstNmrProject(name=nmrProjectName)
if not nmrProject:
ValueError('No NMR project with name %s' % nmrProjectName)
experiment = nmrProject.findFirstExperiment(serial=exptSerial)
if not experiment:
ValueError('Could not find experiment %d in CCPN project.' % (exptSerial))
return
spectrum = experiment.findFirstDataSource(serial=specSerial)
if not spectrum:
ValueError('Could not find spectrum %d in experiment %d.' % (specSerial, exptSerial))
return
peakList = spectrum.findFirstPeakList(serial=plSerial)
if not peakList:
ValueError('Could not find peak list serial %d in spectrum %d:%d' % (plSerial, exptSerial, specSerial))
return peakList
def getCcpnShiftList(ccpnProject, shiftListKey):
"""Descrn: Fetch a CCPN shift list from a CCPN project using a list of object keys
Inputs: Implementation.Project, List of keys (Words or Ints)
Output: ccp.nmr.Nmr.ShiftList
"""
nmrProjectName, shiftListSerial = shiftListKey
nmrProject = ccpnProject.findFirstNmrProject(name=nmrProjectName)
if not nmrProject:
ValueError('No NMR project with name %s' % nmrProjectName)
shiftList = nmrProject.findFirstMeasurementList(serial=shiftListSerial)
if not shiftList:
ValueError('Could not find shift list %d in CCPN project' % (int(shiftListSerial)))
return shiftList
def getCcpnExperimentShiftList(ccpnExperiment):
"""Descrn: Fetch the CCPN shift list that corresponds to a CCPN NMR experiment
Inputs: ccp.nmr.Nmr.Experiment
Output: ccp.nmr.Nmr.ShiftList
"""
shiftList = ccpnExperiment.shiftList
if not shiftList:
for spectrum in ccpnExperiment.dataSources:
for peakList in spectrum.peakLists:
for peak in peakList.peaks:
for peakDim in peak.peakDims:
for contrib in peakDim.peakDimContribs:
shift = contrib.resonance.findFirstShift()
if shift:
return shift.parentList
return shiftList
def getCcpnChain(ccpProject, chainKey):
"""Descrn: Fetch a CCPN chain from a CCPN project using a list of object keys
Inputs: Implementation.Project, List of keys (Words or Ints)
Output: ccp.molecule.MolSystem.Chan
"""
molSystemCode, chainCode = chainKey
molSystem = ccpProject.findFirstMolSystem(code=molSystemCode)
if not molSystem:
ValueError('No molecular system with code "%s" in CCPN project' % molSystemCode)
chain = molSystem.findFirstChain(code=chainCode)
if not chain:
ValueError('No chain found with code "%s" in molecular system "%s"' % (chainCode, molSystemCode))
return chain
# BARDIAUX
def getCcpnChains(ccpProject, chainKey):
"""Descrn: Fetch a CCPN chain from a CCPN project using a list of object keys
Inputs: Implementation.Project, List of keys (Words or Ints)
Output: list of ccp.molecule.MolSystem.Chains
"""
molSystemCode, chainCodes = chainKey[0], chainKey[1:]
molSystem = ccpProject.findFirstMolSystem(code=molSystemCode)
if not molSystem:
ValueError('No molecular system with code "%s" in CCPN project' % molSystemCode)
chains = []
for chainCode in chainCodes:
chain = molSystem.findFirstChain(code=chainCode)
if not chain:
ValueError('No chain found with code "%s" in molecular system "%s"' % (chainCode, molSystemCode))
chains.append(chain)
return chains
def getCcpnPeakAndShiftLists(ccpProject, molSystem, peakShiftListKeys):
"""Descrn: Fetch a list of peak list and corresponding shift list pairs for a CCPN project
given pairs of corresponding object key lists
Inputs: Implementation.Project, List of 2-List of keys (Words or Ints)
Output: List of (ccp.nmr.Nmr.PeakList, ccp.nmr.Nmr.ShiftList)
"""
peakAndShiftLists = []
for peakListKey, shiftListKey in peakShiftListKeys:
shiftlist = None
peakList = getCcpnPeakList(ccpProject, peakListKey)
if peakList:
experiment = peakList.dataSource.experiment
if (experiment.molSystems) and (molSystem not in experiment.molSystems):
data = (molSystem.code, experiment.serial, experiment.name)
print 'WARNING: Selected molecular system "%s" does not match experiment %d (%s) molecular systems' % data
shiftList = getCcpnExperimentShiftList(experiment)
shiftList0 = getCcpnShiftList(ccpProject, shiftListKey)
if shiftList0:
shiftList = shiftList0
if shiftList and (shiftList0 is not shiftList):
print 'WARNING: Selected CCPN shift list does not match experiment %d (%s) shift list' % (experiment.serial, experiment.name)
if peakList and shiftList:
peakAndShiftLists.append([peakList, shiftList])
if not peakAndShiftLists:
print 'ERROR: No selectable CCPN peak lists.'
return []
return peakAndShiftLists
def getCcpnConstraintLists(ccpnProject, restraintsNames):
"""Descrn: Fetch a CCPN constraint list from a CCPN project using list
of restraint type and corresponding object keys
Inputs: Implementation.Project, ARIA restraint names (what object is this?)
Output: List if ccp.nmr.NmrConstraint.AbstractConstraintLists
"""
constraintLists = {}
for restraintType, constraintListKeys in restraintsNames.items():
constraintLists[restraintType] = []
for constraintListKey in constraintListKeys:
keys = getKeysFromString(constraintListKey)
constraintList = getCcpnConstraintList(ccpnProject, keys)
if constraintList:
constraintLists[restraintType].append(constraintList)
return constraintLists
def getCcpnConstraintList(ccpnProject, keys):
"""Descrn: Fetch a CCPN constraint list from a CCPN project using list of object keys
Inputs: Implementation.Project, List of keys (Words or Ints)
Output: ccp.nmr.NmrConstraint.AbstractConstraintList
"""
constraintList = None
nmrProjectKey, storeSerial, serial = keys
nmrProject = ccpnProject.findFirstNmrProject(name=nmrProjectKey)
if not nmrProject:
raise ValueError('No NMR project with name %s' % nmrProjectKey)
constraintStore = nmrProject.findFirstNmrConstraintStore(serial=storeSerial)
if not constraintStore:
raise ValueError('No NMR constraint store with serial "%d" in CCPN project' % storeSerial)
else:
constraintList = constraintStore.findFirstConstraintList(serial=serial)
if constraintList is None:
raise ValueError('No constraint list with serial "%d" in store "%d"' % (serial,storeSerial))
return constraintList
def makeAriaMolecule(ccpMolSystem, ccpChain=None):
"""Descrn: Make an ARIA Molecule given a CCPN MolSystem object
Inputs: ccp.molecule.MolSystem.MolSystem
Output: ARIA Molecule
"""
aria_molecule = Molecule(name=ccpMolSystem.code)
for chain in ccpMolSystem.chains:
if not ccpChain or chain in ccpChain:
aria_chain = makeAriaChain(chain)
aria_molecule.add_chain(aria_chain)
return aria_molecule
def makeAriaChain(ccpChain):
"""Descrn: Make an ARIA Chain given a CCPN Chain object
Inputs: ccp.molecule.MolSystem.Chain
Output: ARIA Chain
"""
# Does below work for DNA/RNA?
aria_settings = ChainSettings()
aria_settings['type'] = chainTypeMapping[ccpChain.molecule.molType]
aria_chain = Chain(settings=aria_settings, segid=string_to_segid(ccpChain.code))
for residue in ccpChain.residues:
aria_residue = makeAriaResidue(residue)
aria_chain.addResidue(aria_residue)
return aria_chain
def makeAriaResidue(ccpResidue):
"""Descrn: Make an ARIA Residue given a CCPN Residue object
Inputs: ccp.molecule.MolSystem.Residue
Output: ARIA Residue
"""
## TODO: have to define chemComp ARIA sysName for residues
## that match, throw error if not recognized (have to define
## own topology file for CNS for now)
try:
aria_residue = Residue(number=ccpResidue.seqCode, residue_type=ccpResidue.ccpCode)
except:
message = 'Residue %s not recognized in ARIA - define topology!'
raise Exception(message % ccpResidue.ccpCode)
ariaAtomDict = {}
idNum = 0
for atom in ccpResidue.atoms:
aria_atom = makeAriaAtom(atom, idNum)
aria_residue.addAtom(aria_atom)
ariaAtomDict[atom.name] = aria_atom
idNum += 1
setupEquivalentGroup(ccpResidue.chemCompVar, aria_residue, ariaAtomDict)
return aria_residue
def makeAriaAtom(ccpAtom, idNum, heteroElements=('N', 'C')):
"""Descrn: Make an ARIA Atom given a CCPN Atom object
Also sets up the hetero atom name.
Inputs: ccp.molecule.MolSystem.Atom
Output: ARIA Atom
"""
chemAtom = ccpAtom.chemAtom
elementSymbol = chemAtom.elementSymbol
heteroName = None
if elementSymbol == 'H':
for bound in chemAtom.findFirstChemBond().chemAtoms:
if bound.elementSymbol in heteroElements:
heteroName = bound.name
break
aria_settings = AtomSettings()
aria_settings['type'] = elementSymbol
aria_settings['hetero_atom_name' ] = heteroName
aria_atom = Atom(settings=aria_settings, name=ccpAtom.name, id=idNum)
aria_atom._setSegid(string_to_segid(ccpAtom.residue.chain.code))
return aria_atom
def makeAriaChemicalShift(ccpShift):
"""Descrn: Make an ARIA ChemicalShift given a CCPN Shift object
Inputs: ccp.nmr.Nmr.Shift
Output: ARIA ChemicalShift
"""
if ccpShift is None:
val = None
err = None
else:
val = ccpShift.value
err = ccpShift.error
return ChemicalShift(val, err)
def makeAriaShiftList(ccpShiftList, ccpMolSystem, ariaMolecule):
"""Descrn: Make a populated ARIA ChemicalShiftList for a given molecule given a CCPN ShiftList object
Inputs: ccp.nmr.Nmr.Shift, ccp.molecule.MolSystem.MolSystem, ARIA Molecule
Output: ARIA ChemicalShiftList
"""
aria_shiftList = ChemicalShiftList()
ss = []
shiftLookup = {}
for shift in ccpShiftList.measurements:
shiftLookup[shift.resonance] = shift
for chain in ccpMolSystem.chains:
for residue in chain.residues:
atomSetDict = {}
for atom in residue.atoms:
atomSet = atom.atomSet
if atomSet and atomSet.resonanceSets:
atomSetDict[atomSet] = True
for atomSet in atomSetDict.keys():
if atomSetDict[atomSet] is False:
continue
resonanceSets = list(atomSet.resonanceSets)
if len(resonanceSets) > 1:
print 'MESSAGE: CCPN atom set %d%s %s has multiple resonance assignments' % (residue.seqCode, residue.ccpCode, atomSet.name)
nonstereo = None
stereo = None
for resonanceSet0 in resonanceSets:
if len(resonanceSet0.atomSets) > 1:
nonstereo = resonanceSet0
else:
stereo = resonanceSet
if stereo:
if nonstereo:
atomSets2 = list(nonstereo.atomSets)
atomSets2.remove(atomSet)
resonanceSet = nonstereo
for resonanceSet0 in atomSets2[0].resonanceSets: # resonance set of other atomSet in pair
if len(resonanceSet0.atomSets) == 1: # if other atom set stereospecifically assigned
resonanceSet = stereo
else:
resonanceSet = stereo
else:
resonanceSet = nonstereo
else:
resonanceSet = resonanceSets[0]
nA = len(resonanceSet.atomSets)
nR = len(resonanceSet.resonances)
if (nA==1) and (nR==1):
ariaAtoms = getAriaAtomsFromResonance(resonanceSet.findFirstResonance(), ariaMolecule)[0]
if len(atomSet.atoms) == 1:
method = AVERAGING_METHOD_NONE
assignment_type = ASSIGNMENT_METHOD_STEREO_SPECIFIC
else:
method = AVERAGING_METHOD_FAST
assignment_type = ASSIGNMENT_METHOD_EQUIVALENT
shift = shiftLookup.get(resonanceSet.findFirstResonance())
if shift:
aria_shift = makeAriaChemicalShift(shift)
aria_spinsystem = SpinSystem(method)
aria_spinsystem.setAtoms(tuple(ariaAtoms))
aria_spinsystem.setChemicalShifts((aria_shift,))
if (aria_spinsystem,) in ss:
continue
ss.append((aria_spinsystem,))
aria_shiftassignment = ShiftAssignment(assignment_type)
aria_shiftassignment.setSpinSystems((aria_spinsystem,))
aria_shiftList.addShiftAssignment(aria_shiftassignment)
atomSetDict[atomSet] = False
elif (nA==2) and (nR<3):
ariaShifts = []
for resonance in resonanceSet.resonances:
shift = shiftLookup.get(resonance)
if shift:
ariaShifts.append(makeAriaChemicalShift(shift))
while len(ariaShifts) < nA:
ariaShifts.append(ChemicalShift(None))
aria_spinsystems = []
for ariaAtoms in getAriaAtomsFromResonance(resonanceSet.findFirstResonance(), ariaMolecule):
if len(ariaAtoms) == 1:
method = AVERAGING_METHOD_NONE
else:
method = AVERAGING_METHOD_FAST
aria_spinsystem = SpinSystem(method)
aria_spinsystem.setAtoms(tuple(ariaAtoms))
aria_spinsystem.setChemicalShifts(tuple(ariaShifts))
aria_spinsystems.append(aria_spinsystem)
atomSetDict[atomSet] = False
if aria_spinsystems in ss:
continue
ss.append(aria_spinsystems)
aria_shiftassignment = ShiftAssignment(ASSIGNMENT_METHOD_FLOATING)
aria_shiftassignment.setSpinSystems(tuple(aria_spinsystems))
aria_shiftList.addShiftAssignment(aria_shiftassignment)
return aria_shiftList
# TJS: Addef option to filter out peaks rejected on CCPN
def makeAriaSpectrum(peakList, ariaMolecule, filterRejected=True):
"""Descrn: Make an ARIA NOESY Spectrum (and assign it)
given a CCPN peak list and ARIA Molecule
Inputs: ccp.nmr.Nmr.PeakList, ARIA Molecule
Output: ARIA NOESYSpectrum
"""
ariaDimNames = ('Proton1','Proton2','Hetero1','Hetero2')
spectrum = peakList.dataSource
experiment = spectrum.experiment
shiftList = experiment.shiftList or experiment.nmrProject.findFirstMeasurementList(className='ShiftList')
transfer = experiment.findFirstExpTransfer(transferType='NOESY')
if not transfer:
raise Exception('Not a NOESY.')
light_dim1 = None
light_dim2 = None
heavy_dim1 = None
heavy_dim2 = None
expDimRefDict = {}
for expDimRef in transfer.sortedExpDimRefs():
if expDimRef.isotopeCodes != ('1H',):
raise Exception('Not an H-H NOESY')
onebondTransfer = expDimRef.findFirstExpTransfer(transferType='onebond')
if onebondTransfer:
expDimRefs = list(onebondTransfer.expDimRefs)
expDimRefs.remove(expDimRef)
expDimRefX = expDimRefs[0]
if light_dim1 is None:
heavy_dim1 = findCcpnDataDim(spectrum, expDimRefX.expDim)
light_dim1 = findCcpnDataDim(spectrum, expDimRef.expDim)
else:
heavy_dim2 = findCcpnDataDim(spectrum, expDimRefX.expDim)
light_dim2 = findCcpnDataDim(spectrum, expDimRef.expDim)
else:
if light_dim1 is None:
heavy_dim1 = None
light_dim1 = findCcpnDataDim(spectrum, expDimRef.expDim)
else:
heavy_dim2 = None
light_dim2 = findCcpnDataDim(spectrum, expDimRef.expDim)
cross_peaks = []
for peak in peakList.peaks:
# TJS added to remove rejected peaks
if filterRejected and peak.figOfMerit == 0.0:
continue
volume = peak.findFirstPeakIntensity(intensityType='volume')
if volume:
value = volume.value
err = volume.error
else:
value = None
err = None
volume = Datum(value, err)
height = peak.findFirstPeakIntensity(intensityType='height')
if height:
value = height.value
err = height.error
else:
value = None
err = None
intensity = Datum(value, err)
cross_peak = CrossPeak(number=peak.serial, volume=volume, intensity=intensity)
peakDims = []
for dim in (light_dim1, light_dim2, heavy_dim1 ,heavy_dim2):
peakDims.append(peak.findFirstPeakDim(dim=dim))
for i in range(len(ariaDimNames)):
peakDim = peakDims[i]
if peakDim:
setShiftFunc = getattr(cross_peak, 'set%sChemicalShift' % ariaDimNames[i])
setShiftFunc(ChemicalShift(peakDim.value, peakDim.valueError))
assignments = []
for peakDim in peakDims:
assignments0 = []
if peakDim:
for contrib in peakDim.peakDimContribs:
resonance = contrib.resonance
for ariaAtoms in getAriaAtomsFromResonance(resonance, ariaMolecule):
assi = Assignment(ariaAtoms, assignment_type=ASSIGNMENT_TYPE_MANUAL)
if assi not in assignments0:
assignments0.append(Assignment(ariaAtoms, assignment_type=ASSIGNMENT_TYPE_MANUAL))
assignments.append(assignments0)
for i in range(len(ariaDimNames)):
addAssignFunc = getattr(cross_peak, 'add%sAssignment' % ariaDimNames[i])
for assignment in assignments[i]:
addAssignFunc(assignment)
cross_peaks.append(cross_peak)
spectrum = NOESYSpectrum(name=experiment.name, noes=cross_peaks)
return spectrum
def getAriaAtomsFromResonance(resonance, ariaMolecule, cache={}):
"""Descrn: Get the corresponding Aria2 atom in a molecule
given the input CCPN Resonance object
Inputs: Nmr.Resonance, Aria2 Molecule object
Output: List of List of Aria2 Atom objects
"""
from tools import string_to_segid
ariaAtoms = cache.get(resonance)
if ariaAtoms:
return ariaAtoms
else:
ariaAtoms = []
resonanceSet = resonance.resonanceSet
if not resonanceSet:
#print 'Attempt to get atoms from unassigned CCPN resonance (%d)' % resonance.serial
return []
atomSets = resonanceSet.sortedAtomSets()
residue = atomSets[0].findFirstAtom().residue
ariaChain = ariaMolecule.getChain(string_to_segid(residue.chain.code)) # Func warns if failure
ariaResidue = ariaChain.residues.get(residue.seqCode)
if not ariaResidue:
print 'Could not find ARIA Residue for CCPN residue %d%s' % (residue.seqCode,residue.ccpCode)
return []
#for atomSet in atomSets:
# - Removed this loop and added below.
# - Loop was giving two ARIA atoms for prochirals: now using only one
# because the restraint re-export was giving two resonances where
# it couldn't tell which was which: both linked to the same two atoms
# TJS Added
index = resonanceSet.sortedResonances().index(resonance)
atomSet = atomSets[min(index,len(atomSets)-1)]
atoms = []
for atom in atomSet.atoms:
ariaAtom = ariaResidue.atoms.get(atom.name)
if not ariaAtom:
print 'Could not find ARIA Atom for CCPN atom %d%s %s' % (residue.seqCode,residue.ccpCode, atom.name)
return []
atoms.append(ariaAtom)
ariaAtoms.append(atoms)
cache[resonance] = ariaAtoms
return ariaAtoms
def findCcpnDataDim(spectrum, expDim):
"""Descrn: Get the data dimension number that corresponds
to a given experimental dimension of a spectrum
Inputs: ccp.nmr.Nmr.DataSource, ccp.nmr.Nmr.ExpDim
Output: Int
"""
dataDim = None
if expDim:
dataDim = spectrum.findFirstDataDim(expDim=expDim)
if dataDim:
return dataDim.dim
# BARDIAUX
def dumpRestraintsList(constraintList, data_dir, type):
"""
Dump CCPN RestraintList to CNS Format TBL file.
"""
import DataContainer as DC
import os
from ccpnmr.format.converters.CnsFormat import CnsFormat
ccpProject = constraintList.root#nmrConstraintHead.project
cns_format = CnsFormat(ccpProject)
f = {DC.DATA_UNAMBIGUOUS : cns_format.writeDistanceConstraints,
DC.DATA_AMBIGUOUS : cns_format.writeDistanceConstraints,
DC.DATA_HBONDS : cns_format.writeHBondConstraints,
DC.DATA_DIHEDRALS : cns_format.writeDihedralConstraints,
DC.DATA_RDCS : cns_format.writeRdcConstraints,
DC.DATA_KARPLUS : cns_format.writeJCouplingConstraints,
DC.DATA_SSBONDS : cns_format.writeDistanceConstraints}
kw = {'compressResonances' : 0,
'minimalPrompts' : 1,
'forceNamingSystem' : 'XPLOR'}
fileName = getObjectKeyString(constraintList, delimiter='_').replace(" ", "")
fileName = "%s.tbl" % fileName
dst = os.path.join(data_dir, fileName)
# assumes that Resonnances are already linked
# overwrite existing file without warning
# attribute 'usePeakInfo' has default True now, but not all ccpn projects has
# this attribute and it's not necessary to ARIA 2 any longer
if type in (DC.DATA_AMBIGUOUS, DC.DATA_UNAMBIGUOUS): #by AWSS
kw['usePeakInfo'] = False #by AWSS
f[type](dst, constraintList = constraintList, **kw)
return dst
# BARDIAUX
def getAriaDistanceRestraintsList(constraint_list, constraint_type, aria_mol):
"""
Convert CCPN DistanceContraintList
to list of ARIA DistanceRestraint
"""
from ShiftAssignment import SpinSystem, AVERAGING_METHOD_NONE
import Contribution as C
from Singleton import SpinPairFactory
from AriaPeak import DistanceRestraint
# BARDIAUX 2.2 NEW
from CrossPeak import CrossPeak
from NOESYSpectrum import ConstraintList
from N import power
from DataContainer import PeakData#, LowerBoundCorrection, UpperBoundCorrection
from Datum import Datum
ccpn_project = constraint_list.root #nmrConstraintHead.project
SpinPair = SpinPairFactory()
restraints = []
peaks = []
#source = getObjectKeyString(constraint_list)
source = "%s %s" % (constraint_list.name, getObjectKeyString(constraint_list))
msg = "WARNING : %s ignored."
for distConstr in constraint_list.constraints:
spin_systems_pairs = []
# BARDIAUX 2.2 NEW
if distConstr.origData:
vol = distConstr.origData
else:
vol = power(distConstr.targetValue, -6)
vol = Datum(vol, 0)
xpk = CrossPeak(distConstr.serial, vol, vol)
contributions = []
restraint = DistanceRestraint()
target, upper, lower = distConstr.targetValue, distConstr.upperLimit, distConstr.lowerLimit
weight = distConstr.weight
restraint.setDistance(target)
restraint.setLowerBound(lower)
restraint.setUpperBound(upper)
restraint.setWeight(weight)
for constrItem in distConstr.items:
reso1, reso2 = constrItem.resonances
ss1 = SpinSystem(AVERAGING_METHOD_NONE)
atoms = getAriaAtomsFromResonance(reso1, aria_mol)
if not atoms:
print msg % (constrItem)
continue
elif len(atoms)> 1:
atoms = map(lambda x: x[0], atoms)
else:
atoms = atoms[0]
ss1.setAtoms(tuple(atoms))
ss2 = SpinSystem(AVERAGING_METHOD_NONE)
atoms = getAriaAtomsFromResonance(reso2, aria_mol)
if not atoms:
print msg % (constrItem)
continue
elif len(atoms)> 1:
atoms = map(lambda x: x[0], atoms)
else:
atoms = atoms[0]
ss2.setAtoms(tuple(atoms))
spin_pairs = []
for a1 in ss1.getAtoms():
for a2 in ss2.getAtoms():
if a1 == a2:
continue
sp = SpinPair(a1, a2)
spin_pairs.append(sp)
if not len(spin_pairs):
continue
if (ss1,ss2) in spin_systems_pairs or \
(ss2,ss1) in spin_systems_pairs :
continue
spin_systems_pairs.append((ss1,ss2))
c = C.Contribution(0,
C.CONTRIBUTION_TYPE_FAST_EXCHANGE,
spin_pairs, spin_systems = (ss1, ss2))
## TODO: hack
c.setSpinSystems((ss1, ss2))
contributions.append(c)
## Need a weight if no further Contribution evaluator
c.setWeight(1.)
if len(contributions) == 0:
print msg % (distConstr)
continue
restraint.setContributions(contributions)
peaks.append(xpk)
restraint.setReferencePeak(xpk)
restraints.append(restraint)
# assign a ARIA ConstraintList as Spectrum
s = ConstraintList(source, peaks)
pd = PeakData()
pd.reset()
pd['ccpn_id'] = getObjectKeyString(constraint_list)
s.setDataSource(pd)
return restraints, s
# return restraints
# BARDIAUX 2.2
def getCcpnExperimentData(peakList):
experiment = peakList.dataSource.experiment
frequency = mixing = None
# Spec freq
spec = experiment.spectrometer
if spec:
frequency = spec.protonFreq
# Mixing time
# CCPN return seconds, ARIA needs ms
transfer = experiment.findFirstExpTransfer(transferType='NOESY')
if transfer:
mixing = transfer.mixingTime
if mixing:
mixing = mixing *1e3
# Correlation Time ?
return frequency, mixing
## class CCPNData(Settings):
## def create(self):
## from Settings import MultiTypeEntity, AbsolutePath, TypeEntity
## from TypeChecking import TUPLE
## d = {}
## d['filename'] = AbsolutePath(exists=0)
## return d
## def create_default_values(self):
## d = {'filename': ''}
## return d
## class CCPNDataXMLPickler(XMLBasePickler):
## order = 'filename',
## def create(self):
## return CCPNData()
## def _xml_state(self, x):
## from xmlutils import XMLElement
## e = XMLElement()
## e.filename = x['filename']
## order = list(self.order)
## e.set_tag_order(order)
## return e
## def load_from_element(self, e):
## from tools import as_tuple
## s = self.create()
## filename = str(e.filename).strip()
## E = s.getEntity('filename')
## E.reset()
## E.mandatory(filename <> '')
## s['filename'] = filename
## return s
## CCPNData._xml_state = CCPNDataXMLPickler()._xml_state
|
