Here it is (and a good point!).
Wayne
On Thu, 24 Nov 2005, Magnus Helgstrand wrote:
> Hi!
>
> Could you please post it on the list. I can't run Analysis, so I cant
> use the upgrade functionality...
>
> Magnus
>
> 2005/11/24, Wayne Boucher <[log in to unmask]>:
> > I've added the offending file, MoleculeBasic, to the server.
> > Unfortunately I don't know whether that will be the only extra file
> > required, but I hope it is (let us know if not!). (Tim is in Sheffield
> > today so will not be able to look at this properly until he gets back.)
> >
> > Wayne
> >
> > On Thu, 24 Nov 2005, Dr Andy Herbert wrote:
> >
> > > Hi,
> > >
> > > I've just installed the StructureBasic.py upgrade using the upgrade
> > > feature of analysis. The result is that now analysis doesn't work with
> > > the following traceback.
> > >
> > > Traceback (most recent call last):
> > > File
> > > "/usr/progs/ccpnmr_floyd/ccpnmr1.0/python/ccpnmr/analysis/AnalysisGui.py", line 62, in ?
> > > from ccpnmr.analysis.AnalysisPopup import AnalysisPopup
> > > File
> > > "/usr/progs/ccpnmr_floyd/ccpnmr1.0/python/ccpnmr/analysis/AnalysisPopup.py", line 92, in ?
> > > from ccpnmr.analysis.BrowseConstraintsPopup import
> > > BrowseConstraintsPopup
> > > File
> > > "/usr/local/ccpnmr.1.0.7/ccpnmr/ccpnmr1.0/python/ccpnmr/analysis/BrowseConstraintsPopup.py", line 68, in ?
> > > from ccpnmr.analysis.StructureBasic import getAtomSetsDistance
> > > File
> > > "/usr/progs/ccpnmr_floyd/ccpnmr1.0/python/ccpnmr/analysis/StructureBasic.py", line 61, in ?
> > > from ccpnmr.analysis.MoleculeBasic import getSequenceResidueMapping
> > > ImportError: cannot import name getSequenceResidueMapping
> > > >>>
> > >
> > >
> > > Replacing with the old version of StructureBasic.py fixes the problem.
> > >
> > > Cheers
> > >
> > > Andy
> > >
> > > --
> > > Dr Andy Herbert
> > > Department of Chemistry
> > > University of Edinburgh
> > > West Mains Road
> > > Edinburgh
> > > UK
> > > EH9 3JJ
> > > Tel: +44 (0)131 650 4704 or 650 7372
> > > Email: [log in to unmask]
> > >
> >
>
"""
======================COPYRIGHT/LICENSE START==========================
MoleculeBasic.py: Part of the CcpNmr Analysis program
Copyright (C) 2005 Wayne Boucher and Tim Stevens (University of Cambridge)
=======================================================================
This file contains reserved and/or proprietary information
belonging to the author and/or organisation holding the copyright.
It may not be used, distributed, modified, transmitted, stored,
or in any way accessed, except by members or employees of the CCPN,
and by these people only until 31 December 2005 and in accordance with
the guidelines of the CCPN.
A copy of this license can be found in ../../../license/CCPN.license.
======================COPYRIGHT/LICENSE END============================
for further information, please contact :
- CCPN website (http://www.ccpn.ac.uk/)
- email: [log in to unmask]
- contact the authors: [log in to unmask], [log in to unmask]
=======================================================================
If you are using this software for academic purposes, we suggest
quoting the following references:
===========================REFERENCE START=============================
R. Fogh, J. Ionides, E. Ulrich, W. Boucher, W. Vranken, J.P. Linge, M.
Habeck, W. Rieping, T.N. Bhat, J. Westbrook, K. Henrick, G. Gilliland,
H. Berman, J. Thornton, M. Nilges, J. Markley and E. Laue (2002). The
CCPN project: An interim report on a data model for the NMR community
(Progress report). Nature Struct. Biol. 9, 416-418.
Wim F. Vranken, Wayne Boucher, Tim J. Stevens, Rasmus
H. Fogh, Anne Pajon, Miguel Llinas, Eldon L. Ulrich, John L. Markley, John
Ionides and Ernest D. Laue (2005). The CCPN Data Model for NMR Spectroscopy:
Development of a Software Pipeline. Proteins 59, 687 - 696.
===========================REFERENCE END===============================
"""
import re
import os
from memops.general.Util import copySubTree
from ccp.api import MolSystem, Nmr, Molecule
from ccpnmr.api import Analysis
from ccp.general.Io import getStdChemCompHeads, getChemCompHead, getChemCompStoragePath
from ccp.general.ChemCompOverview import chemCompStandardOverview
from ccp.util.Molecule import findAtomSetResonances, nextChainCode, makeChain
from ccpnmr.analysis.AssignmentBasic import assignAtomsToRes, assignResonanceResidue, getResidueResonances
#from ccpnmr.analysis.MoleculeBasic import findAtomSetResonances
from memops.gui.MessageReporter import showWarning
BLOSUM62={
'A':{'A': 4,'R':-1,'N':-2,'D':-2,'C': 0,'Q':-1,'E':-1,'G': 0,'H':-2,'I':-1,'L':-1,'K':-1,'M':-1,'F':-2,'P':-1,'S': 1,'T': 0,'W':-3,'Y':-2,'V': 0,'X':0},
'R':{'A':-1,'R': 5,'N': 0,'D':-2,'C':-3,'Q': 1,'E': 0,'G':-2,'H': 0,'I':-3,'L':-2,'K': 2,'M':-1,'F':-3,'P':-2,'S':-1,'T':-1,'W':-3,'Y':-2,'V':-3,'X':0},
'N':{'A':-2,'R': 0,'N': 6,'D': 1,'C':-3,'Q': 0,'E': 0,'G': 0,'H': 1,'I':-3,'L':-3,'K': 0,'M':-2,'F':-3,'P':-2,'S': 1,'T': 0,'W':-4,'Y':-2,'V':-3,'X':0},
'D':{'A':-2,'R':-2,'N': 1,'D': 6,'C':-3,'Q': 0,'E': 2,'G':-1,'H':-1,'I':-3,'L':-4,'K':-1,'M':-3,'F':-3,'P':-1,'S': 0,'T':-1,'W':-4,'Y':-3,'V':-3,'X':0},
'C':{'A': 0,'R':-3,'N':-3,'D':-3,'C': 9,'Q':-3,'E':-4,'G':-3,'H':-3,'I':-1,'L':-1,'K':-3,'M':-1,'F':-2,'P':-3,'S':-1,'T':-1,'W':-2,'Y':-2,'V':-1,'X':0},
'Q':{'A':-1,'R': 1,'N': 0,'D': 0,'C':-3,'Q': 5,'E': 2,'G':-2,'H': 0,'I':-3,'L':-2,'K': 1,'M': 0,'F':-3,'P':-1,'S': 0,'T':-1,'W':-2,'Y':-1,'V':-2,'X':0},
'E':{'A':-1,'R': 0,'N': 0,'D': 2,'C':-4,'Q': 2,'E': 5,'G':-2,'H': 0,'I':-3,'L':-3,'K': 1,'M':-2,'F':-3,'P':-1,'S': 0,'T':-1,'W':-3,'Y':-2,'V':-2,'X':0},
'G':{'A': 0,'R':-2,'N': 0,'D':-1,'C':-3,'Q':-2,'E':-2,'G': 6,'H':-2,'I':-4,'L':-4,'K':-2,'M':-3,'F':-3,'P':-2,'S': 0,'T':-2,'W':-2,'Y':-3,'V':-3,'X':0},
'H':{'A':-2,'R': 0,'N': 1,'D':-1,'C':-3,'Q': 0,'E': 0,'G':-2,'H': 8,'I':-3,'L':-3,'K':-1,'M':-2,'F':-1,'P':-2,'S':-1,'T':-2,'W':-2,'Y': 2,'V':-3,'X':0},
'I':{'A':-1,'R':-3,'N':-3,'D':-3,'C':-1,'Q':-3,'E':-3,'G':-4,'H':-3,'I': 4,'L': 2,'K':-3,'M': 1,'F': 0,'P':-3,'S':-2,'T':-1,'W':-3,'Y':-1,'V': 3,'X':0},
'L':{'A':-1,'R':-2,'N':-3,'D':-4,'C':-1,'Q':-2,'E':-3,'G':-4,'H':-3,'I': 2,'L': 4,'K':-2,'M': 2,'F': 0,'P':-3,'S':-2,'T':-1,'W':-2,'Y':-1,'V': 1,'X':0},
'K':{'A':-1,'R': 2,'N': 0,'D':-1,'C':-3,'Q': 1,'E': 1,'G':-2,'H':-1,'I':-3,'L':-2,'K': 5,'M':-1,'F':-3,'P':-1,'S': 0,'T':-1,'W':-3,'Y':-2,'V':-2,'X':0},
'M':{'A':-1,'R':-1,'N':-2,'D':-3,'C':-1,'Q': 0,'E':-2,'G':-3,'H':-2,'I': 1,'L': 2,'K':-1,'M': 5,'F': 0,'P':-2,'S':-1,'T':-1,'W':-1,'Y':-1,'V': 1,'X':0},
'F':{'A':-2,'R':-3,'N':-3,'D':-3,'C':-2,'Q':-3,'E':-3,'G':-3,'H':-1,'I': 0,'L': 0,'K':-3,'M': 0,'F': 6,'P':-4,'S':-2,'T':-2,'W': 1,'Y': 3,'V':-1,'X':0},
'P':{'A':-1,'R':-2,'N':-2,'D':-1,'C':-3,'Q':-1,'E':-1,'G':-2,'H':-2,'I':-3,'L':-3,'K':-1,'M':-2,'F':-4,'P': 7,'S':-1,'T':-1,'W':-4,'Y':-3,'V':-2,'X':0},
'S':{'A': 1,'R':-1,'N': 1,'D': 0,'C':-1,'Q': 0,'E': 0,'G': 0,'H':-1,'I':-2,'L':-2,'K': 0,'M':-1,'F':-2,'P':-1,'S': 4,'T': 1,'W':-3,'Y':-2,'V':-2,'X':0},
'T':{'A': 0,'R':-1,'N': 0,'D':-1,'C':-1,'Q':-1,'E':-1,'G':-2,'H':-2,'I':-1,'L':-1,'K':-1,'M':-1,'F':-2,'P':-1,'S': 1,'T': 5,'W':-2,'Y':-2,'V': 0,'X':0},
'W':{'A':-3,'R':-3,'N':-4,'D':-4,'C':-2,'Q':-2,'E':-3,'G':-2,'H':-2,'I':-3,'L':-2,'K':-3,'M':-1,'F': 1,'P':-4,'S':-3,'T':-2,'W':11,'Y': 2,'V':-3,'X':0},
'Y':{'A':-2,'R':-2,'N':-2,'D':-3,'C':-2,'Q':-1,'E':-2,'G':-3,'H': 2,'I':-1,'L':-1,'K':-2,'M':-1,'F': 3,'P':-3,'S':-2,'T':-2,'W': 2,'Y': 7,'V':-1,'X':0},
'V':{'A': 0,'R':-3,'N':-3,'D':-3,'C':-1,'Q':-2,'E':-2,'G':-3,'H':-3,'I': 3,'L': 1,'K':-2,'M': 1,'F':-1,'P':-2,'S':-2,'T': 0,'W':-3,'Y':-1,'V': 4,'X':0},
'X':{'A': 0,'R': 0,'N': 0,'D': 0,'C': 0,'Q': 0,'E': 0,'G': 0,'H': 0,'I': 0,'L': 0,'K': 0,'M': 0,'F': 0,'P': 0,'S': 0,'T': 0,'W': 0,'Y': 0,'V': 0,'X':0}
}
def newMolSystem(project):
"""Descrn: Get a new molSystem for a project with a unique code.
Inputs: Implementation.Project
Output: MolSystem.MolSystem
"""
i = 1
while project.findFirstMolSystem(code='MS%d' % (i)):
i += 1
molSystem = project.newMolSystem(code='MS%d' % (i))
molSystem.name = molSystem.code
return molSystem
def transferChainAssignments(chainA, chainB):
"""Descrn: Transfer any atomic assignments from one chain to another where possible.
Inputs: MolSystem.Chain, MolSystem.Chain
Output: None
"""
mapping = getChainResidueMapping(chainA, chainB)
for residueA, residueB in mapping:
resonancesB = getResidueResonances(residueB)
if resonancesB:
showOkCancel('Warning','Destination residue %s%d has assignments. Continue?.' % (residueB.seqCode,residueB.ccpCode))
for residueA, residueB in mapping:
if residueA:
if residueB is None:
showWarning('Warning','Residue %s%d has no equivalent in destination chain' % (residueA.seqCode,residueA.ccpCode))
else:
transferResidueAssignments(residueA,residueB)
def copyMolecule(molecule, newName=None):
"""Descrn: Make a new molecule based upon the sequence of an existing ome
Inputs: Molecule.Molecule
Output: Molecule.Molecule
"""
project = molecule.root
i = len(project.molecules) + 1
newName = newName or 'Molecule %d' % (i)
newMolecule = copySubTree(molecule, project, topObjectParameters={'name':newName,}, maySkipCrosslinks=1 )
return newMolecule
def transferResidueAssignments(residueA,residueB):
"""Descrn: Move any atomic assignments from one residue to another where possible.
Inputs: MolSystem.Residue, MolSystem.Residue
Output: None
"""
resonancesA = getResidueResonances(residueA)
for resonance in resonancesA:
assignResonanceResidue(resonance, residueB)
def getChainResidueMapping(chainA, chainB):
"""Descrn: Find the corresponding pairs of residues in two sequence similar chains.
Inputs: MolSystem.Chain, MolSystem.Chain
Output: List of List of [MolSystem.Residue or None]
"""
seq1 = chainA.molecule.stdSeqString
seq2 = chainB.molecule.stdSeqString
seqA, seqB = sequenceAlign(seq1,seq2,BLOSUM62)
x = 0
y = 0
mapping = []
for i in range(len(seqA)):
mapping.append([None,None])
if seqA[i] != '-':
mapping[i][0] = chainA.residues[x]
x += 1
if seqB[i] != '-':
mapping[i][1] = chainB.residues[y]
y += 1
return mapping
def areAtomsBound(atom1, atom2):
"""Descrn: Dertemine whether two atoms are bonded together
Inputs: MolSystem.Atom, MolSystem.Atom
Output: Boolean
"""
for chemBond in atom1.chemAtom.chemBonds:
if atom2.chemAtom in chemBond.chemAtoms:
return 1
return 0
def areResonancesBound(resonance1,resonance2):
"""Descrn: Determine whether two resonances are assigned to
directly bonded atoms
Inputs: Nmr.Resonance, Nmr.Resonance
Output: Boolean
"""
if resonance1 is resonance2:
return 0
# check both resonances are assigned to atoms
# (else the return is None)
if resonance1.resonanceSet and resonance2.resonanceSet:
residue1 = resonance1.resonanceSet.atomSets[0].atoms[0].residue
residue2 = resonance2.resonanceSet.atomSets[0].atoms[0].residue
# quickly discard on the basis of residue
if residue1 is not residue2:
return 0
for atomSet1 in resonance1.resonanceSet.atomSets:
for atom1 in atomSet1.atoms:
for atomSet2 in resonance2.resonanceSet.atomSets:
for atom2 in atomSet2.atoms:
if areAtomsBound(atom1,atom2):
return 1
return 0
# The default No
def findBoundResonances(resonance):
"""Descrn: Find any resonances which are assigned to atoms covalently bound
to the assigned atoms of the input resonance
Inputs: Nmr.Resonance
Output: List of Nmr.Resonances
"""
boundResonances = ()
if resonance.resonanceSet:
residue = resonance.resonanceSet.atomSets[0].atoms[0].residue
for atomSet in resonance.resonanceSet.atomSets:
for atom in atomSet.atoms:
for atom2 in residue.atoms:
if atom2.atomSet and atom2.atomSet.resonanceSet and (atom2.atomSet is not atomSet):
if areAtomsBound(atom, atom2):
resonances = atom2.atomSet.resonanceSet.resonances
boundResonances.extend(resonances)
return boundResonances
def greekSortAtomNames(dataList, molType='protein'):
"""Descrn: Sorts a list of atom names according to greek/sidechain order when letters are latinised
Inputs: List of Strings
Output: List of Strings (sorted)
"""
sortList = []
for x in dataList:
if type(x) in (tuple,list):
sortName = x[0]
else:
sortName = x
sortName = re.sub('(.+\')', 'zzz@\\1', sortName)
sortName = re.sub('^(\d)','zz@\\1', sortName)
sortName = re.sub('N(\S*)','\'\'@N\\1',sortName)
sortName = re.sub('CO','\'\'@CO',sortName)
sortName = re.sub('g', 'c@g', sortName)
sortName = re.sub('G', 'C@G', sortName)
sortName = re.sub('z', 'f@z', sortName)
sortName = re.sub('Z', 'f@Z', sortName)
if molType == 'protein':
sortName = re.sub('Hn', 'H1n', sortName)
sortList.append( (sortName,x) )
sortList.sort()
dataList = [e[1] for e in sortList]
return dataList
def getMolTypeCcpCodes(molType='all'):
"""Descrn: Gives ccpCodes for chemComps according to molecule type: e.g. DNA
Inputs: String (ChemCompHead.molType or 'all')
Output: List of Words (ChemCompHead.CcpCodes)
"""
ccpCodes = []
if molType == 'all':
molTypes = ['protein','DNA','RNA','other']
else:
molTypes = [molType,]
for molType in molTypes:
if molType == 'other':
from ccp.general.ChemCompOverview import chemCompOtherOverview
ccpCodes.extend( chemCompOtherOverview.keys() )
else:
ccpCodes.extend( chemCompStandardOverview[molType].keys() )
if ccpCodes:
ccpCodes.sort()
return ccpCodes
def makeMolSystemLink(residueA,residueB,linkEndA,linkEndB):
"""Descrn: Make a molSystemLink given two residues and the linkEnds to be joined
Inputs: MolSystem.Residue, MolSystem.Residue, ChemComp.LinkEnd, ChemComp.LinkEnd
Output: MolSystem.MolSystemLink
"""
molSystem = linkAtoms[0].residue.chain.molSystem
project = molSystem.project
molSysLinkEndA = MolSystem.MolSystemLinkEnd(residue,linkCode=linkEndA.linkCode)
molSysLinkEndB = MolSystem.MolSystemLinkEnd(residue,linkCode=linkEndB.linkCode)
molSysLinkEnds = (linkEndA,linkEndB)
try:
molSystemBond = MolSystem.MolSystemLink(molSystem, molSystemLinkEnds=molSysLinkEnds)
except:
return
removeAtoms = []
for atom in removeAtoms:
atomSet = atom.atomSet
if atomSet:
for resonanceSet in atomSet.resonanceSets:
deleteResonanceSet(resonanceset)
atomSet.delete()
atom.delete()
for residue in (residueA,residueB):
residueMapping = getResidueMapping(residue)
for atomSetMapping in residueMapping.atomSetMappings:
atomSetMapping.delete()
makeResidueAtomSets(residue)
return molSystemLink
def makeAtomSet(guiName,atoms,chemAtomSet,mappingType):
"""Descrn: Make atomSet and atomSetMapping for equivalent atoms
Inputs: Word (AtomSet.name), List of Nmr.Atoms, ChemComp.ChemAtomSet, Word (AtomSetMapping.mappingType)
Output: Nmr.AtomSet
"""
project = atoms[0].root
for atom in atoms:
if hasattr(atom, 'atomSet'):
if atom.atomSet != atoms[0].atomSet:
# was made incorrectly previously
for atom2 in atoms:
if atom2.atomSet and (not atom2.atomSet.resonanceSets):
atom2.atomSet.delete()
break
#print guiName,atoms,chemAtomSet,mappingType
#raise 'Tried to make atomSet for non-equivalent atoms'
atom = atoms[0]
if hasattr(atom, 'atomSet'):
if atom.atomSet is None:
atomSet = Nmr.AtomSet(project,atoms=atoms)
else:
atomSet = atom.atomSet
else:
atomSet = Nmr.AtomSet(project,atoms=atoms)
elementSymbol = atom.chemAtom.elementSymbol
residue = atom.residue
residueMapping = getResidueMapping(residue)
if not residueMapping.findFirstAtomSetMapping(name=guiName):
makeAtomSetMapping(residueMapping,guiName,(atomSet,),chemAtomSet,mappingType)
atomSet.name = guiName
return atomSet
def getResidueMapping(residue, aromaticsEquivalent=1):
"""Descrn: Gives the Analysis.ResidueMapping for a residue
Makes a new one with new AtomSetsMappings if not exists
Makes a ChainMapping too if needed.
Inputs: Nmr.Residue
Output: Analysis.ResidueMapping
"""
if hasattr(residue, 'residueMapping'):
return residue.residueMapping
residueMapping = None
chainMapping = residue.root.findFirstChainMapping(chain=residue.chain)
if not chainMapping:
chainMapping = Analysis.ChainMapping(residue.root, molSystemCode=residue.chain.molSystem.code,chainCode=residue.chain.code)
chainMapping.residueMappingDict = {}
else:
if not hasattr(chainMapping, 'residueMappingDict'):
chainMapping.residueMappingDict = {}
residueMapping = chainMapping.residueMappingDict.get(residue.seqId)
if not residueMapping:
residueMapping = chainMapping.findFirstResidueMapping(seqId=residue.seqId)
chainMapping.residueMappingDict[residue.seqId] = residueMapping
if not residueMapping:
residueMapping = Analysis.ResidueMapping(chainMapping, seqId=residue.seqId)
#makeResidueAtomSets(residue)
residue.residueMapping = residueMapping
if not residueMapping.atomSetMappings:
makeResidueAtomSets(residue)
return residueMapping
def makeAtomSetMapping(residueMapping,name,atomSets,chemAtomSet,mappingType,resonances=None):
"""Descrn: Make atomSetMapping given atomSets and mapping type
Inputs: Analysis.ResidueMapping, Word, Nmr.Residue, List of Nmr.AtomSets, ChemComp.ChemAtomSet, Word, Word, List of Nmr.Resonances
Output: Analysis.AtomSetMapping
"""
atom = atomSets[0].atoms[0]
elementSymbol = atom.chemAtom.elementSymbol
serials = []
for atomSet in atomSets:
serials.append(atomSet.serial)
molType = residueMapping.residue.molResidue.molType
guiName = makeGuiName(name, elementSymbol, molType)
atomSetMapping = residueMapping.newAtomSetMapping(name=guiName,mappingType=mappingType,atomSetSerials=serials,
chemAtomSet=chemAtomSet,elementSymbol=elementSymbol)
if resonances is not None:
resSerials = []
for resonance in resonances:
resSerials.append(resonance.serial)
atomSetMapping.setResonanceSerials(resSerials)
return atomSetMapping
def makeGuiName(name, elementSymbol, molType):
"""Descrn: Convert atom or atomSet name into name for gui: e.g H becomes Hn
Inputs: Word (Nmr.AtomSet.name), Word, Word
Output: Word
"""
if molType == 'protein':
letters = re.sub('\d|\*','',name)
if letters == 'H':
name = name[0:1] + 'n' + name[1:]
guiName = elementSymbol + name[len(elementSymbol):].lower()
return guiName
def makeResidueLocTag(residue, chemAtoms):
"""Descrn: Make unique identifier for a given residue type in a given chain location
Inputs: Nmr.Residue, List of ChemCmp.chemAtoms
Output: Word
"""
# check if any atoms have been deleted (to make a molSystemLink)
tag = residue.ccpCode + residue.linking
if len(chemAtoms) != len(residue.atoms):
names = []
for chemAtom in chemAtoms:
names.append(chemAtom.name)
for atom in residue.atoms:
if atom.name in names:
names.remove(atom.name)
for name in names:
tag = tag + '-' + name
return tag
def makeResidueAtomSetsNonEquivalent(residue):
"""Descrn: Remake a residue's atom sets if they are found to be non-equivalent: e.g.
if an aromatic ring does not rotate quickly on the NMR timescale
Inputs: Nmr.Residue
Output: None
"""
residueMapping = getResidueMapping(residue)
project = residue.root
chain = residue.chain
molType = residue.molResidue.molType
nonequivalent = {}
elementSymbolDict = {}
chemAtomSetDict = {}
for atom in residue.atoms:
chemAtom = atom.chemAtom
chemAtomSetDict[atom] = chemAtom
elementSymbol = chemAtom.elementSymbol
if chemAtom.chemAtomSet:
chemAtomSet = chemAtom.chemAtomSet
name = chemAtomSet.name
if chemAtomSet.isEquivalent is None: # i.e. not False, aromatic rotation
if nonequivalent.get(name) is None:
nonequivalent[name] = []
nonequivalent[name].append(atom)
elementSymbolDict[name] = chemAtom.elementSymbol
chemAtomSetDict[name] = chemAtomSet
for groupName in nonequivalent.keys():
atoms = nonequivalent[groupName]
atomSet = atoms[0].atomSet
if atomSet:
for atom in atoms[1:]:
if atom.atomSet is not atoms[0].atomSet:
return # already seperate
elementSymbol = elementSymbolDict[groupName]
chemAtomSet = chemAtomSetDict[groupName]
resonances = []
if atomSet:
for resonanceSet in atomSet.resonanceSets:
for resonance in resonanceSet.resonances:
resonances.append(resonance)
atomSet.delete()
name = makeGuiName(groupName, elementSymbol, molType)
atomSetMapping = residueMapping.findFirstAtomSetMapping(name=name)
if atomSetMapping:
atomSetMapping.delete()
atomSets = []
atomSetNames = []
for atom in atoms:
name = chemAtomSetDict[atom].name
atomSet = makeAtomSet(name,(atom,),chemAtomSet,'stereo')
atomSets.append(atomSet)
atomSetNames.append(name)
resonanceSet = None
for resonance in resonances:
resonanceSet = assignAtomsToRes(atomSets,resonance,resonanceSet)
n = 0
for atom in atoms:
name = chemAtomSetDict[atom].name
name2 = makeNonStereoName(molType, name, n)
n += 1
makeGuiMultiAtomSet(residue, name2, atomSetNames,elementSymbol,'nonstereo',chemAtomSet)
makeGuiMultiAtomSet(residue, groupName, atomSetNames,elementSymbol,'ambiguous',chemAtomSet)
def makeResidueAtomSetsEquivalent(residue):
"""Descrn: Remake a residue's atom sets if they are found to be equivalent: e.g.
if an aromatic ring rotates quickly on the NMR timescale
Inputs: Nmr.Residue
Output: None
"""
project = residue.root
chain = residue.chain
molType = residue.molResidue.molType
residueMapping = getResidueMapping(residue)
equivalent = {}
elementSymbolDict = {}
chemAtomSetDict = {}
for atom in residue.atoms:
chemAtom = atom.chemAtom
chemAtomSetDict[atom] = chemAtom
elementSymbol = chemAtom.elementSymbol
if chemAtom.chemAtomSet:
chemAtomSet = chemAtom.chemAtomSet
name = chemAtomSet.name
if chemAtomSet.isEquivalent is None: # i.e. not False, aromatic rotation
if equivalent.get(name) is None:
equivalent[name] = []
equivalent[name].append(atom)
elementSymbolDict[name] = chemAtom.elementSymbol
chemAtomSetDict[name] = chemAtomSet
for groupName in equivalent.keys():
atoms = equivalent[groupName]
if atoms[0].atomSet:
for atom in atoms[1:]:
if atom.atomSet is atoms[0].atomSet:
return
elementSymbol = elementSymbolDict[groupName]
chemAtomSet = chemAtomSetDict[groupName]
resonances = []
for atom in atoms:
# TBD more nested layers?
# delete ambiguous
name = makeGuiName(chemAtomSet.name, elementSymbol, molType)
atomSetMapping = residueMapping.findFirstAtomSetMapping(name=name)
if atomSetMapping:
atomSetMapping.delete()
# delete stereospecific
name = makeGuiName(atom.chemAtom.name, elementSymbol, molType)
atomSetMapping = residueMapping.findFirstAtomSetMapping(name=name)
if atomSetMapping:
atomSetMapping.delete()
# delete non-stereospecific
name = makeGuiName(makeNonStereoName(molType, atom.chemAtom.name), elementSymbol, molType)
atomSetMapping = residueMapping.findFirstAtomSetMapping(name=name)
if atomSetMapping:
atomSetMapping.delete()
atomSet = atom.atomSet
if atomSet:
for resonanceSet in atomSet.resonanceSets:
resonances.append(resonanceSet.resonances)
atomSet.delete()
elementSymbol = elementSymbolDict[groupName]
chemAtomSet = chemAtomSetDict[groupName]
# make single equivalent group
#makeAtomSet(,groupName,atoms,chemAtomSet,'simple')
atomSet = makeAtomSet(groupName,atoms,chemAtomSet,'simple')
for resonance in resonances:
assignAtomsToRes([atomSet,],resonance)
def makeResidueAtomSets(residue, aromaticsEquivalent=1):
"""Descrn: Make all atomSets and atomSetMappings for a given residue
Aromatic Phe, Tyr (Hd1,Hd2), (He1,He2) can be made into
single equivalent atom sets due to rotation.
Inputs: Nmr.Residue, Boolean
Output: None
"""
residueMapping = getResidueMapping(residue)
equivalent = {}
elementSymbolDict = {}
nonequivalent = {}
multiSet = {}
chemAtomSetDict = {}
inMultiSet = {}
molType = residue.molResidue.molType
for atom in residue.atoms:
chemAtom = atom.chemAtom
chemAtomSetDict[atom] = chemAtom
elementSymbol = chemAtom.elementSymbol
if chemAtom.chemAtomSet is None:
name = chemAtom.name
makeAtomSet(name,(atom,),None,'simple')
else:
chemAtomSet = chemAtom.chemAtomSet
name = chemAtomSet.name
elementSymbolDict[name] = elementSymbol
chemAtomSetDict[name] = chemAtomSet
if chemAtomSet.isEquivalent:
if equivalent.get(name) is None:
equivalent[name] = []
equivalent[name].append(atom)
elif (chemAtomSet.isEquivalent is None) and atom.atomSet and (len(atom.atomSet.atoms) > 1):# aromatic rotation prev set
if equivalent.get(name) is None:
equivalent[name] = []
equivalent[name].append(atom)
elif (chemAtomSet.isEquivalent is None) and (not atom.atomSet) and aromaticsEquivalent:# aromatic rotation to be set
if equivalent.get(name) is None:
equivalent[name] = []
equivalent[name].append(atom)
else:
if nonequivalent.get(name) is None:
nonequivalent[name] = []
nonequivalent[name].append(atom)
if chemAtomSet.chemAtomSet is not None:
multiName = chemAtomSet.chemAtomSet.name
chemAtomSetDict[multiName] = chemAtomSet.chemAtomSet
elementSymbolDict[multiName] = elementSymbol
if multiSet.get(multiName) is None:
multiSet[multiName] = {}
multiSet[multiName][name] = 1
inMultiSet[name] = multiName
for groupName in equivalent.keys():
atoms = equivalent[groupName]
elementSymbol = elementSymbolDict[groupName]
chemAtomSet = chemAtomSetDict[groupName]
if len(atoms)==2:
# not enough atoms for multi sets!
makeAtomSet(groupName,atoms,chemAtomSet,'simple')
else:
if inMultiSet.get(groupName):
# e.g. for Val Hg1*
makeAtomSet(groupName,atoms,chemAtomSet,'stereo')
else:
makeAtomSet(groupName,atoms,chemAtomSet,'simple')
for groupName in nonequivalent.keys():
atoms = nonequivalent[groupName]
elementSymbol = elementSymbolDict[groupName]
chemAtomSet = chemAtomSetDict[groupName]
atomSetNames = []
for atom in atoms:
name = chemAtomSetDict[atom].name
makeAtomSet(name,(atom,),chemAtomSet,'stereo')
atomSetNames.append(name)
n = 0
for atom in atoms:
name = chemAtomSetDict[atom].name
name2 = makeNonStereoName(molType, name, n)
n += 1
makeGuiMultiAtomSet(residue, name2, atomSetNames,elementSymbol,'nonstereo',chemAtomSet)
makeGuiMultiAtomSet(residue, groupName, atomSetNames,elementSymbol,'ambiguous',chemAtomSet)
for groupName in multiSet.keys():
atomSetNames = multiSet[groupName].keys()
elementSymbol = elementSymbolDict[groupName]
chemAtomSet = chemAtomSetDict[groupName]
if "|" in groupName:
# we don't do these pseudoatoms in Analysis
continue
# e.g. for Val Hga*
n = 0
for atomSetName in atomSetNames:
name2 = makeNonStereoName(molType, atomSetName, n)
n += 1
makeGuiMultiAtomSet(residue, name2, atomSetNames, elementSymbol,'nonstereo',chemAtomSet)
makeGuiMultiAtomSet(residue, groupName, atomSetNames,elementSymbol,'ambiguous',chemAtomSet)
def makeNonStereoName(molType, name, n=None):
"""Descrn: Convert a sterospecific atom name into a non-stereospecific one for a GUI
Inputs: Word, Int (naming offset from start of alphabet)
Output: Word
"""
match = re.match('(\w+)(\d|\'+)(\D*)', name)
if not match:
#print molType, name, n
return name
letters = match.group(1)
number = match.group(2)
prime = ''
if number == '\'':
number = 1
prime = '\''
elif number == '\'\'':
number = 2
prime = '\''
if n is None:
n = int(number) - 1
if molType == 'protein':
if (letters == 'H'):
letters = 'Hn'
name = letters + prime + chr(ord('a')+n)+ match.group(3)
return name
def makeGuiMultiAtomSet(residue,multiGuiName,guiSetsNames,elementSymbol,mappingType,chemAtomSet):
"""Descrn: Make atom set mappings for multiple atom set selections
Inputs: Nmr.Residue, Word (Analysis.AtomSetMapping.name), List of Words (Analysis.AtomSetMapping.names), Word, Word, ChemComp.ChemAtomSet
Output: Analysis.AtomSetMapping
"""
if "|" in multiGuiName:
return
residueMapping = getResidueMapping(residue)
molType = residue.molResidue.molType
for guiName in guiSetsNames:
atomSetMapping = residueMapping.findFirstAtomSetMapping(name=makeGuiName(guiName, elementSymbol, molType))
if atomSetMapping is None:
print "Nonexistant group error in makeGuiMultiAtomSet for", residue.molResidue.ccpCode, residue.seqCode, guiName
return
#atomSet = atomSetMapping.atomSets[0]
chemAtomSet1 = atomSetMapping.chemAtomSet
for guiName2 in guiSetsNames:
atomSetMapping2 = residueMapping.findFirstAtomSetMapping(name=makeGuiName(guiName2, elementSymbol, molType))
if atomSetMapping2 is None:
print "Nonexistant group error in makeGuiMultiAtomSet for", residue.molResidue.ccpCode
return
#atomSet = atomSetMapping2.atomSets[0]
chemAtomSet2 = atomSetMapping2.chemAtomSet
if chemAtomSet2 and chemAtomSet1:
if chemAtomSet1.isProchiral != chemAtomSet2.isProchiral:
print "Prochiratity error in makeGuiMultiAtomSet for", residue.molResidue.ccpCode
return
if chemAtomSet1.isEquivalent != chemAtomSet2.isEquivalent:
print "Equivalenct arror error makeGuiMultiAtomSet for ", residue.molResidue.ccpCode
return
atomSets = []
for guiName in guiSetsNames:
atomSetSerials = residueMapping.findFirstAtomSetMapping(name=makeGuiName(guiName, elementSymbol, molType)).atomSetSerials
for atom in residue.atoms:
atomSet = atom.atomSet
if atomSet:
if atomSet.serial in atomSetSerials:
atomSets.append(atomSet)
break
if not residueMapping.findFirstAtomSetMapping(name=multiGuiName):
atomSetMapping = makeAtomSetMapping(residueMapping, multiGuiName,atomSets,chemAtomSet,mappingType)
return atomSetMapping
def moveMolSystemChain(chain, molSystem):
"""Descrn: Moves a chain from one molSystem to another
Inputs: MolSystem.Chain, MolSystem.MolSystem
Output: None
"""
from MergeObjects import mergeObjects
if chain.molSystem is molSystem:
return
molecule = chain.molecule
project = molSystem.root
startNum = chain.residues[0].seqCode
seq = []
code = nextChainCode(molSystem)
molType = chain.molecule.molType
for residue in chain.residues:
seq.append( residue.ccpCode )
newChain = makeChain(molSystem,molecule)
for chainMapping in project.chainMappings:
if chainMapping.chain is chain:
chainMapping.delete()
residues1 = list(chain.residues)
for i in range(len(seq)):
residue1 = residues1[i]
residue2 = newChain.residues[i]
for atom in residue1.atoms:
atom2 = residue2.findFirstAtom(name=atom.name)
if atom2:
mergeObjects(atom, atom2)
else:
print "missing %d %s %s" % (residue.seqCode,residue.ccpCode,atom.name)
mergeObjects(residue1, residue2)
mergeObjects(chain, newChain)
def sequenceAlign(seq1,seq2,matrix,treshold=50,inspen=3,extpen=1,minScore=-10):
"""Descrn: Aligns two sequence strings (of one letter codes) - results are gapped with '-'
Inputs: String, String, Dict of Dicts (homology score matrix), Int, Int, Int, Int
Output: String, String (aligned & gapped sequence strings)
"""
seq1 = re.sub('-|\s','',seq1)
seq2 = re.sub('-|\s','',seq2)
seq1 = re.sub('\*','X',seq1)
seq2 = re.sub('\*','X',seq2)
nX = len(seq1)
nY = len(seq2)
M = [[0] * nY]
R = [[3] * nY]
#maximum = 0
#bestPoint = [0,0]
route = 0
for x in range(1,nX):
M.append([0,])
R.append([2,])
for y in range(1,nY):
p1 = inspen
p2 = inspen
if route == 2:
p1 = extpen
elif route == 3:
p2 = extpen
score = matrix[seq1[x]][seq2[y]]
try:
paths = [minScore, M[x-1][y-1]+score, M[x-1][y]-p1, M[x][y-1]-p2]
except:
raise '%d %d %d %d %d %d' % (len(M),len(M[0]),x-1,y-1,x,y)
best = max(paths)
route = paths.index(best)
M[x].append( best )
R[x].append( route )
#if best >= maximum:
# maximum = best
# bestPoint = [x,y]
x = nX-1
y = nY -1
score = M[x][y]
route = R[x][y]
seqA = ''
seqB = ''
while (score > minScore) and ( x > 0 or y > 0 ):
if route == 1:
seqA = seq1[x] + seqA
seqB = seq2[y] + seqB
x -=1
y -=1
elif route == 2:
seqB = '-' + seqB
seqA = seq1[x] + seqA
x -=1
elif route == 3:
seqA = '-' + seqA
seqB = seq2[y] + seqB
y -=1
score = M[x][y]
route = R[x][y]
del M
del R
return seqA, seqB
|
