Source code for ccpnmodel.ccpncore.lib.V2Upgrade

"""Mapping of Resonances and ResonancGroups in version 2 to new 4-string assignment style

"""
#=========================================================================================
# Licence, Reference and Credits
#=========================================================================================
__copyright__ = "Copyright (C) CCPN project (http://www.ccpn.ac.uk) 2014 - 2017"
__credits__ = ("Wayne Boucher, Ed Brooksbank, Rasmus H Fogh, Luca Mureddu, Timothy J Ragan & Geerten W Vuister")
__licence__ = ("CCPN licence. See http://www.ccpn.ac.uk/v3-software/downloads/license",
               "or ccpnmodel.ccpncore.memops.Credits.CcpnLicense for licence text")
__reference__ = ("For publications, please use reference from http://www.ccpn.ac.uk/v3-software/downloads/license",
               "or ccpnmodel.ccpncore.memops.Credits.CcpNmrReference")

#=========================================================================================
# Last code modification
#=========================================================================================
__modifiedBy__ = "$modifiedBy: CCPN $"
__dateModified__ = "$dateModified: 2017-07-07 16:33:09 +0100 (Fri, July 07, 2017) $"
__version__ = "$Revision: 3.0.0 $"
#=========================================================================================
# Created
#=========================================================================================

__author__ = "$Author: CCPN $"
__date__ = "$Date: 2017-04-07 10:28:48 +0000 (Fri, April 07, 2017) $"
#=========================================================================================
# Start of code
#=========================================================================================

import os
import operator

from ccpnmodel.ccpncore.lib.chemComp.Util import chemAtomSetFromAtoms


[docs]def mapResonanceGroupResidues(apiNmrProject, molSystem=None, chainMap=None) -> dict: """Map resonanceGroup:assignmentTuple for fully assigned ResonanceGroups """ result = {} # We need either chainMap or molSystem, and chainMap takes precedence if chainMap: molSystem = None # Handle properly assigned ResonanceGroups for resonanceGroup in apiNmrProject.sortedResonanceGroups(): # Find residue residue = resonanceGroup.residue if residue is None: ll = [] for residueProb in resonanceGroup.residueProbs: if residueProb.weight: ll.append(residueProb) if len(ll) == 1: # In principle this should never happen, but it does not hurt to be careful residue = ll[0].possibility if residue: # Remap chain and residue if necessary, chain = residue.chain if chainMap: # If there is a chainMap the chain MUST be in it. chain = chainMap[chain] # get residue in new chain residue = chain.findFirstResidue(seqCode=residue.seqCode, seqInsertCode=residue.seqInsertCode) elif residue.topObject is not molSystem: raise ValueError("Cannot generate consistent assignment names from mixed MolSystems - 1") # # # set residue assignment strings # chainCode = chain.code # sequenceCode = str(residue.seqCode)+ (residue.seqInsertCode or '').strip() # residueType = residue.code3Letter if residue: # Should nevere happen, bunless there is an error in teh cahinMapping result[resonanceGroup] = residue # return result
[docs]def mapUnAssignedFixedResonances(nmrConstraintStore): """Map unassigned resonances for NmrConstraintStore NBNB This must be done AFTER assignments and resonances are done, as it transfers teh assignment of the attached NmrProject resonance (if any) to the fixedResonance""" nmrProject = nmrConstraintStore.nmrProject result = {} separator1 = '@' separator2 = '@@' # To distinguish fixedResonance serial from resonance serial for resonance in nmrConstraintStore.sortedFixedResonances(): if not resonance.resonanceSet: # unassigned - treat it nmrResonance = nmrProject.findFirstResonance(serial=resonance.resonanceSerial) if nmrResonance: resonanceGroup = nmrResonance.resonanceGroup result[resonance] = (resonanceGroup.nmrChain.code, resonanceGroup.sequenceCode, resonanceGroup.residueType, nmrResonance.name) else: # No resonance found - make name from fixedResonance info name = regularisedResonanceName(resonance) # Add resonance serial to name, if it is not in the name already. separator = separator1 # use resonanceSerial if available serial = resonance.resonanceSerial if not serial: # use FixedResonance.serial instead, and use '@@' to distinguish serial = resonance.serial separator = separator2 ss = '%s%s' % (separator1, serial) if ss not in name: name = '%s%s%s' % (name, separator, serial) # result[resonance] = (None, None, None, name) # return result
[docs]def addOffsetResonanceGroup(addToGroup:'ResonanceGroup', addGroup:'ResonanceGroup', offset:int) -> bool: """Add addGroup as satellite to addToGroup with offset offset. Return True if successful""" previous = addToGroup.nmrProject.findFirstResonanceGroup(mainGroupSerial=addToGroup.serial, relativeOffset=offset) if previous is None: # New. Set as offset. Multisteps to avoid name clashes and ensure undo. # We do not need to guard against connected stretches etc. in V2 conversion addGroup.sequenceCode = None addGroup.directNmrChain = addToGroup.nmrChain addGroup.sequenceCode = '%s%+d' % (addToGroup.sequenceCode, offset) # set offset residueType residueType = addGroup.residueType if residueType is None: chemComp = addToGroup.root.findFirstChemComp(molType=addGroup.molType, ccpCode=addGroup.ccpCode) if chemComp: addGroup.residueType = chemComp.code3Letter return True else: # Duplicate. Merge. # NBNB check that resonance names are dealt with properly later for resonance in addGroup.resonances: resonance.resonanceGroup = previous addGroup.delete() return False
[docs]def mapAssignedResonances(topObject, molSystem=None, chainMap=None): """Make/extend {resonance:assignmentTuple} map in V2 for either Resonances or fixedResonances chainMap remaps chains to new ones with different chainCodes (for V2-V3 upgrade). NB, does NOT use ResonanceGroup information""" result = {} if topObject.className == 'NmrProject': resonanceSets = topObject.sortedResonanceSets() else: resonanceSets = topObject.sortedFixedResonanceSets() for resonanceSet in resonanceSets: # set up loop-level parameters resonances = list(resonanceSet.resonances) if resonanceSet.findFirstResonance(name=None): resonances.sort(key=operator.attrgetter('serial')) else: resonances.sort(key=operator.attrgetter('name')) atomSets = list(resonanceSet.atomSets) if resonanceSet.findFirstAtomSet(name=None): atomSets.sort(key=operator.attrgetter('serial')) else: atomSets.sort(key=operator.attrgetter('name')) allAtoms = [x for y in atomSets for x in y.atoms] chemAtomSet = chemAtomSetFromAtoms(allAtoms) if len(set(x.residue for x in allAtoms)) == 1: residue = allAtoms[0].residue else: residue = None if residue: # Remap chain and residue if necessary chain = residue.chain if chainMap: # If there is a chainMap the chain MUST be in it. chain = chainMap[chain] # get residue in new chain residue = chain.findFirstResidue(seqCode=residue.seqCode, seqInsertCode=residue.seqInsertCode) elif chain.molSystem is not molSystem: raise ValueError("Cannot generate consistent assignment names from mixed MolSystems - 2") if residue: # Should always be true - but in case something went wrong with the chain mapping # We have the residue - and we need the ChemCompVar below chemComp = residue.molResidue.chemComp chemCompVar = residue.chemCompVar or chemComp.findFirstChemCompVar(isDefaultVar=True) # Now for the atom name if chemAtomSet and len(atomSets) == 2 and len(resonances) <= 2: # prochiral pair. Use _getAmbigProchiralLabel for priority, and set 'x' for 'a', 'y' for 'b' # get non-stereo names atomSetNames = [x.name for x in atomSets] starpos = chemAtomSet.name.find('*') newNames = [] if [x[starpos:] for x in sorted(atomSetNames)] == ["", "'"]: # special case - names are "abc'", "abc''" ss = atomSetNames[0][:starpos-1] newNames = [ss + 'x', ss + 'y'] else: for ii, newChar in enumerate('xy'): chars = list(atomSetNames[ii]) try: chars[starpos] = newChar except Exception as err: from ccpn.util.Logging import getLogger getLogger().warn(f'Impossible to convert to new nomenclature. AtomSet:{atomSetNames[ii]}, Error:{err}') newNames.append(''.join(chars)) # select new name to use if topObject.className == 'NmrProject': # Real resonance - use normal procedure # First one resonance. NB we test prochiral label on one resonance only # - if data are inconsistent both might give 'a' and we must avoid a name clash. if _getAmbigProchiralLabel(resonances[0]) != 'a': # First resonance does not match first name - reverse name order # NB we test on one resonance only. If data are inconsistent and both give 'a' # we still want the names to be different newNames.reverse() for ii,resonance in enumerate(resonances): result[resonance] = (residue, newNames[ii].replace('*','%')) else: # NmrConstraintStore - these are fixed resonances if len(resonances) == 2: # resonances are sorted by name, as are newNames. Match in order. # NB this being in an NmrConstraintStore, we have to assume that names are # consistent, so e.g. HG1* is bound to CG1 and not CG2 for ii in range(2): result[resonances[ii]] = (residue, newNames[ii].replace('*','%')) else: # Only one resonance. Must choose which. # Use various heuristics, as we can not assume assignments in NmrProject still match. # assert len(resonances) == 1 resonance = resonances[0] if resonance.name is None: # Aribtratily pck the first one. Wnat else? indx = 0 else: for ii in range(2): if (resonance.name in(atomSetNames[ii], newNames[ii]) or resonance.name.upper() in(atomSetNames[ii], newNames[ii])): # name matches one of the atomSets - use matching name indx = ii break else: # Name did not match either possibility, so these are non-standard. Try anyway if 'b' in resonance.name or '3' in resonance.name: # Heuristic - these names are likely to be second in sorting # e.g. 'HBb' or 'HB3' indx = 1 elif 'a' in resonance.name or '1' in resonance.name: # Heuristic = this name is likely to be first in sorting (e.g. 'HBa' # NB both XY1, XY1*, XY11, XY21 (Asn, Gln,Arg) XY12 (Ile) should sort first # Arg HH12 should sort second, but this cannot be helped indx = 0 elif '2' in resonance.name: # This will not sort second for e,g, HB2/HB3, but that cannot be helped. # The only one that matters is # where e.g. HD1 must bind to CD1 and HD2 to CD2, and it gets those right # Anyway the real atom names are caught before this indx = 1 else: if not hasattr(resonance, 'resonance'): from ccpn.util.Logging import getLogger getLogger().warn( f'Impossible to convert to new nomenclature. resonance:{resonance} has not attribute resonance.') continue realResonance = resonance.resonance if realResonance: # We can not be sure that assignments have not changed, but # better have it match the resonance than not if _getAmbigProchiralLabel(realResonance) == 'a': indx = 0 else: indx = 1 else: # Stuff it. We just do not know and pick the first one. # Anyway the cases where is makes a difference are cared for above. indx = 0 resonanceName = newNames[indx] result[resonance] = (residue, resonanceName.replace('*','%')) elif len(resonances) == 1: # Single resonance, not assigned to prochiral resonance = resonances[0] if len(atomSets) == 1: # simple one-to-one stereospecific assignment if chemAtomSet: # assignment to atomSet resonanceName = chemAtomSet.name.replace('*', '%') else: #asssignment to single atom resonanceName = allAtoms[0].name else: # multiple atomSets # NB we do it this way because it must work in pure V2 as well as the intermediate model nuc = allAtoms[0].elementSymbol residueChemAtoms = [chemCompVar.findFirstChemAtom(name=x.name) for x in residue.atoms] residueChemAtoms = [x for x in residueChemAtoms if x is not None] if len(allAtoms) == len(residueChemAtoms): # All single atoms in residue resonanceName = '*' elif (all((x.elementSymbol == nuc) for x in allAtoms) and len([x for x in residueChemAtoms if x.elementSymbol == nuc]) == len(allAtoms)): # All atoms of a given nucleus resonanceName = nuc + '*' else: # random multiple atom selection resonanceName = '/'.join(sorted((str(x.name) for x in atomSets))) result[resonance] = (residue, resonanceName) else: # multiple resonances not matching chemAtomSet atomsName = '/'.join(sorted((str(x.name) for x in atomSets))) for resonance in resonances: # NB this name can not be in use already, so we do not need to check resonanceName = '%s@%s' % (atomsName, resonance.serial) result[resonance] = (residue, resonanceName.replace('*','%')) else: # assigned to multiple residues - cannot be helped # Same naming style for single and multiple resonances partNames = [] for atomSet in atomSets: rr = atomSet.findFirstAtom().residue partNames.append('%s%s-%s' % (rr.seqCode, (rr.seqInsertCode or '').strip(), atomSet.name)) ss = '/'.join(sorted(partNames)) for resonance in resonances: resonanceName = '%s@%s' % (ss, resonance.serial) result[resonance] = (None, resonanceName.replace('*','%')) # return result
################################################################################### # # Functions for V2/upgrade only: # ###################################################################################
[docs]def regularisedResonanceName(resonance): """V2: Get resonance name, starting with element type and adding @serial to impossible names""" # NB names like '*', 'C*' etc. will not make it through regularisation # but it is too dangerous to arrive at those from name strings only # For that you need to have them properly assigned. resonanceName = resonance.name or '' # get elementCode isotope = resonance.isotope if isotope is None: elementCode = '?' else: elementCode = isotope.chemElement.symbol.upper() result = None if resonance.className == 'Resonance': # Exclude fixedResonances assignNames = list(set(resonance.assignNames)) if len(assignNames) == 1: # One assignName - use for assignment name = assignNames[0] result = name.replace('*', '%') elif len(assignNames) == 2: # 2 assignNames - if they match prochiral return nonstereo variant assignNames.sort() prefix = os.path.commonprefix(assignNames) lenPrefix = len(prefix) if assignNames[0][lenPrefix].isdigit() and assignNames[1][lenPrefix].isdigit(): if assignNames[0][lenPrefix+1:] == assignNames[1][lenPrefix+1:]: # Heuristics, try to get x or y right # NB 'x' will not always match the lowest sorting - it fails for # ordinary HB2/HB3 methylenes and must be caught later # But the crucial isopropyl and aromatic cases should work newChar = None for char in reversed(resonanceName): if char.isdigit(): break else: char = None if assignNames[0][lenPrefix] == char: newChar = 'x' elif assignNames[1][lenPrefix] == char: newChar = 'y' elif 'a' in resonanceName and not 'b' in resonanceName: newChar = 'x' elif 'b' in resonanceName and not 'a' in resonanceName: newChar = 'y' if newChar is not None: ll = list(assignNames[0]) ll[lenPrefix] = newChar result = ''.join(ll).replace('*', '%') if result is None: # If we are still here, assignNames did not help. Use resonanceName only if resonanceName: if elementCode and resonanceName.upper().startswith(elementCode): # name is OK except possibly for casing - fix the casing ss = resonanceName[len(elementCode):] if 'x' in ss or 'y' in ss: # Necessary to avoid potential clashes with XY names set above # No actual assigned names contain 'x' or 'y' anyway result = '%s@%s-%s' % (elementCode, resonance.serial, resonanceName) else: # Name might be proper assignment name. Change to new wildcard convention result = (elementCode + ss).replace('*', '%') else: # Set unique default name result = '%s@%s-%s' % (elementCode, resonance.serial, resonanceName) elif hasattr(resonance, 'resonanceSerial') and resonance.resonanceSerial: result = '%s@%s' % (elementCode, resonance.resonanceSerial) else: result = '%s@%s' % (elementCode, resonance.serial) # return result
[docs]def upgradeConstraintList(constraintList): """Upgrade ConstraintList from early V3 to newer V3 - this avoids redoing earlier function and anyway data must be copied to a new set of objects Will also work if called on old-type V3 ConstraintLists (use only internally)""" constraintStore = constraintList.nmrConstraintStore # Get defining parameters className = constraintList.className restraintType = className[:-14] if restraintType in ('Distance', 'HBond', 'JCoupling', 'Rdc',): itemLength = 2 elif restraintType in ('Csa', 'ChemicalShift'): itemLength = 1 elif restraintType == 'Dihedral': itemLength = 4 else: raise ValueError("Restraint list named %s not recognized by code (BUG2?)" % className) # Make new ConstraintList params = {'constraintType':restraintType, 'itemLength':itemLength} for tag in ('name', 'details', 'potentialType', 'unit', 'usedForCalculation', 'experimentSerial', 'measureListSerial', 'origin', 'tensorIsotropicValue', 'tensorMagnitude', 'tensorRhombicity', 'tensorChainCode', 'tensorSequenceCode', 'tensorResidueType'): if hasattr(constraintList, tag): val = getattr(constraintList, tag) if val is not None: params[tag] = getattr(constraintList, tag) # RESETS: # reset name to unique string, to free the name to give to new ConstraintList: constraintList.name = '@@@%s' % constraintList.serial if restraintType == 'HBond': # RESET: HBond lists are type 'Distance' origin 'hbond' restraintType = params['constraintType'] = 'Distance' params['origin'] = 'hbond' newList = constraintStore.newGenericConstraintList(**params) constraintStore.root.override = True constraintTags = ('serial', 'origData', 'details') contributionTags = ('serial', 'weight', 'targetValue', 'error', 'upperLimit', 'lowerLimit', 'additionalUpperLimit', 'additionalLowerLimit', 'combinationId', ) contributionTags2 = ('scale', 'isDistanceDependent') try: for constraint in constraintList.sortedConstraints(): # Make new Constraint params = {} for tag in constraintTags: if hasattr(constraint, tag): val = getattr(constraint, tag) if val is not None: params[tag] = val # params = {tag:getattr(constraint,tag) for tag in constraintTags} if hasattr(constraint, 'vectorLength'): params['vectorLength'] = constraint.vectorLength newConstraint = newList.newGenericConstraint(**params) for peakContrib in constraint.sortedPeakContribs(): dd = {} for tag in ('experimentSerial', 'dataSourceSerial', 'peakListSerial', 'peakSerial'): dd[tag] = getattr(peakContrib, tag) newConstraint.newConstraintPeakContrib(**dd) for contribution in constraint.sortedContributions(): # Make new Contribution params = {} for tag in contributionTags + contributionTags2: if hasattr(contribution, tag): val = getattr(contribution, tag) if val is not None: params[tag] = val # params = {tag:getattr(contribution,tag) for tag in contributionTags} # for tag in contributionTags2: # if hasattr(contribution, tag): # params[tag] = getattr(contribution, tag) newContribution = newConstraint.newGenericContribution(**params) for item in contribution.items: # Make new Item if itemLength == 1: newContribution.newSingleAtomItem(resonance=item.resonance) elif itemLength == 4: newContribution.newFourAtomItem(resonances=item.resonances) else: # assert itemLength == 2 newContribution.newAtomPairItem(resonances=item.resonances) # delete old list and copy across serial serial = constraintList.serial tempSerial = newList.serial constraintList.delete() newList.__dict__['serial'] = serial constraintStore.__dict__['constraintLists'][serial] = newList del constraintStore.__dict__['constraintLists'][tempSerial] constraintStore.__dict__['_serialDict']['constraintLists'] -= 1 finally: constraintStore.root.override = False
[docs]def findSpinSystemStretch(resonanceGroup, excludedSpinSystems=set()): """Find (one of the) longest sequential spin system stretch(es) containing resonanceGroup """ # To avoid modifying external data excludedSpinSystems = excludedSpinSystems.copy() excludedSpinSystems.add(resonanceGroup) stretches = [] for delta in (-1,1): queue = [resonanceGroup] dd = {resonanceGroup:[]} for rg in queue: ll = dd[rg] for rg2 in findConnectedSpinSystems(rg, delta): if rg2 not in excludedSpinSystems: excludedSpinSystems.add(rg2) queue.append(rg2) ll2 = dd[rg2] = ll.copy() ll2.append(rg2) stretch = [] for ll3 in sorted(dd.values()): if len(ll3) > len(stretch): stretch = ll3 # stretches.append(stretch) # result = list(reversed(stretches[0])) + [resonanceGroup] + stretches[1] return result
[docs]def findConnectedSpinSystems(spinSystem, delta=None): """ Find spin systems sequentially connected to the input one with given sequence offset. .. describe:: Input Nmr.ResonanceGroup, Int .. describe:: Output Nmr.ResonanceGroup """ result = [] for link in getSeqSpinSystemLinks(spinSystem, delta=delta): if spinSystem is link.parent: result.append(link.possibility) else: result.append(link.parent) # return result
[docs]def findIdentityResonanceGroup(resonanceGroup): """V2: find unique resonanceGroup linked as identical to teh input""" ll = resonanceGroup.findAllResonanceGroupProbs(linkType='identity', isSelected=True) identicals = list(x.possibility for x in ll) ll2 = resonanceGroup.findAllFromResonanceGroups(linkType='identity', isSelected=True) identicals.extend(x.parent for x in ll2) if len(identicals) == 1: return identicals[0] else: return None
def _getAmbigProchiralLabel(resonance): """ V2: Deterimine if an ambigous prochiral resonance (non-stereospecifically assigned) Has an "a" label or a "b" label. "a" is reserved for the upfield proton and any other nulceus bound to it. .. describe:: Input Nmr.Resonance .. describe:: Output Character """ letter = '' if hasattr(resonance, 'onebond'): del resonance.onebond resonanceSet = resonance.resonanceSet if resonanceSet: if resonance.isotopeCode == '1H': data = [] for resonance2 in resonanceSet.sortedResonances(): if resonance2.shifts: data.append( (resonance2.findFirstShift().value, resonance2.serial, resonance2) ) else: data.append( (999999.999, resonance2.serial, resonance2) ) data.sort() resonances = [x[2] for x in data] i = resonances.index(resonance) letter = chr(ord('a')+i) else: resonance2 = _getOnebondResonance(resonance, isotopeCode='1H') if resonance2 and resonance2.resonanceSet and (len(resonance2.resonanceSet.atomSets) > 1): letter = _getAmbigProchiralLabel(resonance2) resonance2.onebond = resonance elif (len(resonanceSet.resonances) > 1) and (len(resonanceSet.atomSets) > 1): for resonance2 in resonanceSet.resonances: if resonance2 is not resonance: resonance3 = _getOnebondResonance(resonance2) if (resonance3 and resonance3.resonanceSet and (len(resonance3.resonanceSet.atomSets) > 1)): letter = 'b' break if not letter: data = [] for resonance2 in resonanceSet.resonances: if resonance2.shifts: data.append( (resonance2.findFirstShift().value,resonance2) ) else: data.append( (resonance2.serial,resonance2) ) data.sort() resonances = [x[1] for x in data] i = resonances.index(resonance) letter = chr(ord('a')+i) #keyword = 'ambigProchiralLabel' #appData = resonance.findFirstApplicationData(application=app, keyword=keyword) # #if appData and (appData.value != letter): # appData.delete() # appData = None # #if not appData: # AppDataString(resonance,application=app,keyword=keyword, value=letter) return letter def _getOnebondResonance(resonance, isotopeCode=None): """ V2: Find any resonance that may have a single bond connection to the input resonance Option to specify the isotope type .. describe:: Input Nmr.Resonance, Nmr.Resonance.isotopeCode .. describe:: Output Nmr.Resonance """ resonances = _getBoundResonances(resonance) if resonances: if isotopeCode: for resonance1 in resonances: if resonance1.isotopeCode == isotopeCode: return resonance1 else: return resonances[0] resonance2 = None for contrib in resonance.peakDimContribs: peakDim = contrib.peakDim expDimRef1 = peakDim.dataDimRef.expDimRef expTransfers = expDimRef1.expTransfers for expTransfer in expTransfers: if expTransfer.transferType in ('onebond','CP'): expDimRef2 = None for expDimRef in expTransfer.expDimRefs: if expDimRef is not expDimRef1: expDimRef2 = expDimRef break if expDimRef2: if (not isotopeCode) or (isotopeCode in expDimRef2.isotopeCodes): for peakDim2 in peakDim.peak.peakDims: if peakDim2.dataDimRef and (peakDim2.dataDimRef.expDimRef is expDimRef2): for contrib2 in peakDim2.peakDimContribs: if (not isotopeCode) or (contrib2.resonance.isotopeCode == isotopeCode): resonance2 = contrib2.resonance break if resonance2: break return resonance2 def _getBoundResonances(resonance, recalculate=False, contribs=None, recursiveCall=False): """ V2: Find all resonances that have a single bond connection to the input resonance Option to recalculate given assignment status (e.g. if something changes) Option to specify peakDimContribs to search .. describe:: Input Nmr.Resonance, Boolean, List of Nmr.PeakDimContribs .. describe:: Output List of Nmr.Resonances """ if (not recalculate) and resonance.covalentlyBound: return list(resonance.covalentlyBound) resonances = set() # Linked by bound atoms irrespective of spectra pairResonances = set() # prochiral or other pairs that can not be determined imemdiately resonanceSet = resonance.resonanceSet funnyResonances = set() if resonanceSet: #residue = resonanceSet.findFirstAtomSet().findFirstAtom().residue atomSets = resonanceSet.atomSets for atomSet in atomSets: #for atom in atomSet.atoms: atom = atomSet.findFirstAtom() for atom2 in _getBoundAtoms(atom): atomSet2 = atom2.atomSet if atomSet2 and atomSet2.resonanceSets: usePaired = False if len(atomSets) > 1: chemAtomSet = atom2.chemAtom.chemAtomSet if chemAtomSet: usePaired = (chemAtomSet.isProchiral or (chemAtomSet.chemAtomSet and chemAtomSet.chemAtomSet.isProchiral)) for resonanceSet2 in atomSet2.resonanceSets: for resonance2 in resonanceSet2.resonances: if resonance2 is resonance: # should not happen if resonance not in funnyResonances: resonance.root._logger.warning( 'in _getBoundResonances():' 'resonance %d tried to be linked to itself' % resonance.serial) funnyResonances.add(resonance) elif usePaired: pairResonances.add(resonance2) else: resonances.add(resonance2) if not contribs: contribs = resonance.peakDimContribs expResonances = set() foundBothPaired = False for contrib in contribs: peakDim = contrib.peakDim expDimRef1 = peakDim.dataDimRef.expDimRef expTransfers = expDimRef1.expTransfers for expTransfer in expTransfers: if expTransfer.transferType == 'onebond': expDimRef2 = None for expDimRef in expTransfer.expDimRefs: if expDimRef is not expDimRef1: expDimRef2 = expDimRef break if expDimRef2: for peakDim2 in peakDim.peak.peakDims: if peakDim2.dataDimRef and (peakDim2.dataDimRef.expDimRef is expDimRef2): expBound = set() for contrib2 in peakDim2.peakDimContribs: if (not contrib.peakContribs) and (not contrib2.peakContribs): resonance2 = contrib2.resonance if resonance is not resonance2: expBound.add(resonance2) else: for peakContrib in contrib.peakContribs: if peakContrib in contrib2.peakContribs: resonance2 = contrib2.resonance if resonance is not resonance2: expBound.add(resonance2) break if len(expBound) > 1: # Ambiguity for bound in expBound: # Leave the covalently bound one if bound in resonances: break else: aSet = set(x for x in expBound if x in resonance.covalentlyBound) if aSet and aSet != pairResonances: # Resonances found. Previously linked. # Not the pairResonances. Use them expResonances.update(aSet) else: # check presence of prochiral pairs ll = [x for x in pairResonances if x in expBound] if len(pairResonances) == 2 and len(ll) == 2: foundBothPaired= True elif ll: # found some prochiral pair resonances - use them expResonances.update(ll) else: expResonances.update(expBound) if foundBothPaired and not [x for x in expResonances if x in pairResonances]: # particular special case. # Resonance is bound to both prochiral alternatives but always as a pair. if recursiveCall: # This was called from elsewhere. We could resolve nothing, so send back to caller pass else: # call for sister resonances and see resons = resonanceSet.sortedResonances() newResonances = set() if len(resons)> 1: # there are sister resonances resons.remove(resonance) for reson in resons: boundResons = _getBoundResonances(reson, recalculate=True, contribs=contribs, recursiveCall=True) ll = [x for x in pairResonances if x not in boundResons] if not ll: # One sister was bound to both. Incorrect data. Bind to both here too newResonances.update(pairResonances) break elif len(ll) < len(pairResonances): # Some resonances were taken. Use the free ones. newResonances.update(ll) if newResonances: expResonances.update(newResonances) else: # No data anywhere to resolve which is which. Match on serials pairResonList = list(sorted(pairResonances, key=operator.attrgetter('serial'))) rr = pairResonList[resonanceSet.sortedResonances().index(resonance)] expResonances.add(rr) resonances.update(expResonances) #if doWarning and (resonance.isotopeCode == '1H') and (len(resonances) > 1): # pass if resonances: resonance.setCovalentlyBound(resonances) else: resonance.setCovalentlyBound([]) return list(resonances) def _getBoundAtoms(atom): """Get a list of atoms bound to a given atom. NB ONLY for use in upgrade .. describe:: Input MolSystem.Atom .. describe:: Output List of MolSystem.Atoms """ if hasattr(atom, 'boundAtoms'): return atom.boundAtoms atoms = [] chemAtom = atom.chemAtom residue = atom.residue chemAtomDict = {} for atom2 in residue.atoms: # Only atoms specific to ChemCompVar :-) chemAtomDict[atom2.chemAtom] = atom2 for chemBond in chemAtom.chemBonds: for chemAtom2 in chemBond.chemAtoms: if chemAtom2 is not chemAtom: atom2 = chemAtomDict.get(chemAtom2) if atom2: atoms.append(atom2) linkEnd = residue.chemCompVar.findFirstLinkEnd(boundChemAtom=chemAtom) if linkEnd: molResLinkEnd = residue.molResidue.findFirstMolResLinkEnd(linkEnd=linkEnd) if molResLinkEnd: molResLink = molResLinkEnd.molResLink if molResLink: for molResLinkEnd2 in molResLink.molResLinkEnds: if molResLinkEnd2 is not molResLinkEnd: residue2 = residue.chain.findFirstResidue(molResidue=molResLinkEnd2.molResidue) if residue2: chemAtom2 = molResLinkEnd2.linkEnd.boundChemAtom atom2 = residue2.findFirstAtom(chemAtom=chemAtom2) if atom2: atoms.append(atom2) break atom.boundAtoms = atoms return atoms