# Copyright 2007-2009 Nanorex, Inc. See LICENSE file for details. """ PAM_Atom_methods.py - methods for class Atom, which are only meant to be called on PAM Atoms. @author: Bruce, Mark, Ninad @version: $Id$ @copyright: 2007-2009 Nanorex, Inc. See LICENSE file for details. History: bruce 080327 split out of chem.py class Atom, in which these methods had been written TODO: REVIEW whether any of these methods are sometimes called on non-PAM atoms, always doing nothing or returning a null value, e.g. getDnaBaseName. also refactor some code that remains in chem.py but has PAM-element-specific sections, like the transmute context menu entries for PAM elements, and some bond_geometry_error_string things (which needs a separate refactoring first, for other reasons) """ # WARNING: this module is imported by chem.py and therefore indirectly by # chunk.py. Therefore its import dependence on other things needs to be # minimized. However, it seems ok at this point to import quite a bit of dna # updater code. Today I moved all imports to toplevel (including those with # old warnings that that didn't work when they were written) and experienced # no problems. # Even so, further refactoring (which among other things would mitigate this) # is desirable. (The ultimate issue is having PAM-specific methods on class # Atom, as opposed to on a PAM-specific subclass of class Atom; and having # that class Atom known directly to class Chunk.) [bruce 090119] from Numeric import dot from geometry.VQT import norm, V ##from foundation.state_constants import S_CHILDREN from utilities import debug_flags from utilities.debug import print_compact_stack from utilities.constants import Pl_STICKY_BOND_DIRECTION, MODEL_PAM5 from utilities.constants import average_value from model.bond_constants import DIRBOND_CHAIN_MIDDLE from model.bond_constants import DIRBOND_CHAIN_END from model.bond_constants import DIRBOND_NONE from model.bond_constants import DIRBOND_ERROR from model.bond_constants import find_bond from model.elements import Pl5 from model.global_model_changedicts import _changed_otherwise_Atoms from dna.model.pam3plus5_math import baseframe_abs_to_rel from dna.model.pam3plus5_math import baseframe_rel_to_abs from dna.model.pam3plus5_math import default_Pl_relative_position from dna.model.pam3plus5_math import default_Gv_relative_position from dna.model.pam3plus5_ops import Pl_pos_from_neighbor_PAM3plus5_data from dna.model.pam3plus5_ops import _f_baseframe_data_at_baseatom from dna.model.pam3plus5_ops import _f_find_new_ladder_location_of_baseatom from dna.model.pam3plus5_ops import kill_Pl_and_rebond_neighbors from dna.model.pam_conversion_mmp import Fake_Pl from dna.updater.dna_updater_prefs import legal_numbers_of_strand_neighbors_on_axis from dna.updater.fix_bond_directions import PROPOGATED_DNA_UPDATER_ERROR from dna.updater.fix_bond_directions import _f_detailed_dna_updater_error_string VALID_ELEMENTS_FOR_DNABASENAME = ('Ss5', 'Ss3', 'Sh5', 'Se3', 'Sj5', 'Sj3',) # TODO: use a boolean element attribute instead. # review: not the same as role == 'strand'... but is it if we exclude Pl? def make_vec_perp_to_vec(vec1, vec2): #bruce 080405 # todo: refile, or find an existing func the same """ return a modified copy of vec1 with its component parallel to vec2 subtracted out. """ remove_this_direction = norm(vec2) return vec1 - dot(vec1, remove_this_direction) * remove_this_direction #k _GV5_DATA_INDEX = 2 # == class PAM_Atom_methods: """ Intended only as a pre-mixin for class Atom. """ # Someday these methods will be refactored into subclasses of Atom. # Then those subclasses can be moved inside the dna package. # (Maybe this module could be moved there right now? Not easily, # regarding package import conventions/worries.) # === # # The Atom methods below this point might be moved into subclasses for # specific types of atoms. # # (Some of these methods might need trivial default defs on class Atom # until old code is fully revised to only call them on the subclasses.) # # Several general methods which remain in class Atom # also have special cases that might be # revised to be in subclass methods that extend them. These include: # # drawing_color # _draw_atom_style # draw_atom_sphere # draw_wirespheres # max_pixel_radius # writemmp # readmmp_info_atom_setitem # getInformationString # # (But some of these might be purely graphical, perhaps usable by more than # one subclass, and might thus remain in the superclass, or in a separately # refactored Atom drawing controller class.) # # [bruce 071113 comment, and reordered existing methods to move these here] # # === # default values of instance variables (some not needed): _dna_updater__error = "" # intentionally not undoable/copyable # this is now higher up, with an undo _s_attr decl: ## _dnaBaseName = "" #bruce 080319 revised # note: _dnaStrandId_for_generators is set when first demanded, or can be # explicitly set using setDnaStrandId_for_generators(). _fake_Pls = None # see comments where used _DnaLadder__ladder = None # for private use by class DnaLadder # note: not named with _f_, so I can use the pseudo-name-mangling # convention to indicate ownership by objects of a different class # [bruce 080413 added class definition, revised using code] # not: _s_attr__nonlive_Pls = S_CHILDREN # for default value assignments of _PAM3plus5_Pl_Gv_data and # _f_Pl_posn_is_definitive, and comments about them, see class Atom # (which inherits this class); they are only used by methods in this # class, except for a few foundational methods in class Atom # (for copy and soon for mmp save/load), but since it uses them # at all, we define them there. [bruce 080523] # # _PAM3plus5_Pl_Gv_data = ... # # _f_Pl_posn_is_definitive = ... # == methods for either strand or axis PAM atoms def dna_updater_error_string(self, include_propogated_error_details = True, newline = '\n' ): #bruce 080206; 080214 also covers check_bond_geometry ### REVIEW: integrate with bad_valence etc """ Return "" if self has no dna updater error (recorded by the dna updater in the private attribute self._dna_updater__error), or an error string if it does. By default, the error string is expanded to show the source of propogated errors from elsewhere in self's basepair (assuming the updater has gotten through the step of propogating them, which it does immediately after assigning/updating all the direct per-atom error strings). Any newlines in the error string (which only occur if it was expanded) are replaced with the optional newline argument (by default, left unchanged). @param include_propogated_error_details: see main docstring @type include_propogated_error_details: boolean @param newline: see main docstring @type newline: string @see: helper functions like _atom_set_dna_updater_error which should eventually become friend methods in a subclass of Atom. """ if self.check_bond_geometry(): # bruce 080214; initial kluge -- # if this happens, ignore the rest (non-ideal) return self.check_bond_geometry() # (redoing this is fast) if not self._dna_updater__error: # report errors from self's DnaLadder, if it has one ladder = getattr(self.molecule, 'ladder', None) # kluge for .ladder if ladder and ladder.error: return ladder.error return "" # otherwise report error from self res = self._dna_updater__error if include_propogated_error_details and \ res == PROPOGATED_DNA_UPDATER_ERROR: res = _f_detailed_dna_updater_error_string(self) res = res.replace('\n', newline) # probably only needed in then clause return res def _pam_atom_cmenu_entries(self): #bruce 080401 """ Return a menu_spec list of context menu entries specific to self being a PAM atom, or None. """ assert self.element.pam ladder = self.molecule.ladder res = [] if ladder: dnaladder_menu_spec = ladder.dnaladder_menu_spec(self) if dnaladder_menu_spec: if res: # never happens yet res.append(None) res.extend(dnaladder_menu_spec) # more? pass return res def reposition_baggage_using_DnaLadder(self, dont_use_ladder = False, only_bondpoints = False ): #bruce 080404; added options, bruce 080501 """ Reposition self's bondpoints (or other monovalent "baggage" atoms if any), making use of structural info from self.molecule.ladder. @note: as of late 080405, this doesn't yet actually make use of self.molecule.ladder. (Perhaps it could be optimized by making it do so?) @note: since we don't expect any baggage except bondpoints on PAM atoms, we don't bother checking for that, even though, if this did ever happen, it might be arbitrary which baggage element ended up at what position. This can easily be fixed if necessary. @warning: by default, it's only safe to call this when self.molecule.ladder is correct (including all its chunks and rails -- no need for its parent Group structure or markers), which if the dna updater is running means after self's new ladder's chunks were remade. @param dont_use_ladder: if true, don't assume self.molecule.ladder is correct. (As of 080501, this option has no effect, since the current implem doesn't use self.molecule.ladder or anything else fixed by the dna updater, except possibly bond directions on open bonds.) @type dont_use_ladder: boolean @param only_bondpoints: if true, only reposition bondpoints, not other baggage. @type only_bondpoints: boolean """ # And, if it turns out this ever needs to look inside # neighbor atom dnaladders (not in the same base pair), # we'll need to revise the dna updater's indirect caller # to call it in a separate loop over ladders than the one # that calls remake_chunks on all of the new ladders. # Note: I don't think there's any issue of the changes we do here # causing an infinite repeat (or even a single repeat) of the # dna updater -- we reset the flag that makes it call us, # and moving bondpoints (or real atoms) doesn't count as # changing the atoms for purposes of dna updater, and even # if it did, the updater ignores new changes from this step. # If we ever kill and remake bondpoints here, then we'll # need to also reset the flag at the end, but even if we forget, # it won't cause an immediate rerun of the dna updater # (due to its ignoring the changes from this step). self._f_dna_updater_should_reposition_baggage = False # even in case of exceptions (also simplifies early return) del dont_use_ladder # correct, since present implem never uses it baggage, others = self.baggage_and_other_neighbors() if only_bondpoints: # bruce 080501 for n in baggage[:]: if not n.is_singlet(): baggage.remove(n) others.append(n) continue pass if not baggage: return # optimization and for safety; might simplify following code # Notes: # # - If baggage might not be all bondpoints, we might want to filter # it here and only look at the bondpoints. (See docstring.) # [Now this is done by the only_bondpoints option. [bruce 080501]] # # - I'm assuming that Ax, Pl, Ss neighbors can never be baggage. # This depends on their having correct valence, since the test # used by baggage_and_other_neighbors is len(bonds), not valence. # how to do it depends on type of PAM element: if self.element.role == 'axis': # Ax3, Ax5, or Gv5 # which real neighbors do we have? # (depends on ladder len1, single strand or duplex) # ... # short term important case: we have one bondpoint # along axis, make it opposite the other axis bond, # and maybe we also have another "strand bondpoint". strand_neighbors = [] axis_neighbors = [] for other in others: # would also work if baggage was in this loop role = other.element.role if role == 'strand': strand_neighbors += [other] elif role == 'axis': axis_neighbors += [other] # other roles are possible too, but we can ignore them for now # (unpaired-base(?), None) # (in future we'll need to count unpaired-bases to understand # how many axis bondpoints are missing) ### TODO continue if len(axis_neighbors) == 1: # move a bondpoint to the opposite position from this one # (for now, just move the already-closest one) selfpos = self.posn() axisvec = axis_neighbors[0].posn() - selfpos newpos = selfpos - axisvec # correct only in direction, not distance moved_this_one = self.move_closest_baggage_to( newpos, baggage, remove = True, min_bond_length = 3.180 / 2 ) # this first corrects newpos for distance (since no option # prevents that) (an option does set the minimum distance # -- after 080405 fix in _compute_bond_params this is no # longer needed, but remains as a precaution) # (doing it first is best for finding closest one) # (someday we might pass bond order or direction to help # choose? then bond order would affect distance, perhaps) # passing baggage means only consider moving something # in that list; # remove = True means remove the moved one from that list; # return value is the one moved, or None if none was found # to move. if not moved_this_one: return # optim? pass ### todo: also fix bondpoints meant to go to strand atoms pass elif self.element.role == 'strand': # might be Pl5 or Ss5 or Ss3 if self.element.symbol == 'Ss3': ### TODO: split this section into its own method # if we have one real strand bond, make the open directional # bond opposite it (assume we have correct structure of # directional bonds with direction set -- one in, one out) strand_neighbors = [] axis_neighbors = [] for other in others: role = other.element.role if role == 'strand': strand_neighbors += [other] elif role == 'axis': axis_neighbors += [other] continue honorary_strand_neighbor = None # might be set below reference_strand_neighbor = None # might be set below if len(strand_neighbors) == 1: reference_strand_neighbor = strand_neighbors[0] elif not strand_neighbors and not axis_neighbors: # totally bare (e.g. just deposited) reference_strand_neighbor = self.next_atom_in_bond_direction(1) # should exist, but no need to check here if reference_strand_neighbor is not None: direction_to_it = find_bond(self, reference_strand_neighbor).bond_direction_from(self) moveme = self.next_atom_in_bond_direction( - direction_to_it ) if moveme is not None and moveme in baggage: honorary_strand_neighbor = moveme selfpos = self.posn() strandvec = reference_strand_neighbor.posn() - selfpos newpos = selfpos - strandvec # correct only in direction, not distance moved_this_one = self.move_closest_baggage_to( newpos, [moveme] ) if moved_this_one is not moveme: # should never happen print "error: moved_this_one %r is not moveme %r" % \ (moved_this_one, moveme) pass pass if len(axis_neighbors) == 0: # see if we can figure out which baggage should point towards axis candidates = list(baggage) if honorary_strand_neighbor is not None: # it's in baggage, by construction above candidates.remove(honorary_strand_neighbor) if reference_strand_neighbor is not None \ and reference_strand_neighbor in candidates: candidates.remove(reference_strand_neighbor) if len(candidates) != 1: # I don't know if this ever happens (didn't see it in test, # can't recall if I thought it should be possible) ## print " for %r: candidates empty or ambiguous:" % self, candidates pass else: # this is the honorary axis neighbor; we'll repoint it. honorary_axis_neighbor = candidates[0] # but which way should it point? just "not along the strand"? # parallel to the next strand atom's? that makes sense... # we'll start with that if we have it (otherwise with # whatever it already is), # then correct it to be perp to our strand. axisvecs_from_strand_neighbors = [] for sn in strand_neighbors: # real ones only its_axisbond = None # might be set below its_rung_bonds = [b for b in sn.bonds if b.is_rung_bond() ] if len(its_rung_bonds) == 1: its_axisbond = its_rung_bonds[0] elif not its_rung_bonds: its_runglike_bonds = \ [b for b in sn.bonds if not b.bond_direction_from(sn) ] # might be rung bonds, might be open bonds with # no direction set if len(its_runglike_bonds) == 1: its_axisbond = its_runglike_bonds[0] if its_axisbond is not None: # that gives us the direction to start with axisvec_sn = its_axisbond.spatial_direction_from(sn) # but correct it to be perpendicular to our bond to sn axisvec_sn = make_vec_perp_to_vec( axisvec_sn, sn.posn() - self.posn() ) axisvecs_from_strand_neighbors += [axisvec_sn] continue axisvec = None # might be set below if axisvecs_from_strand_neighbors: axisvec = average_value(axisvecs_from_strand_neighbors) if axisvec is None: # i.e. no strand_neighbors axisvec = honorary_axis_neighbor.posn() - self.posn() # this value (which it has now) would not change it # print " for %r: axisvec %r starts as" % (self, axisvec,) # now remove the component parallel to the presumed # strand bonds' average direction presumed_strand_baggage = filter( None, [reference_strand_neighbor, honorary_strand_neighbor] ) vecs = [norm(find_bond(self, sn).bond_direction_vector()) for sn in presumed_strand_baggage] if vecs: strand_direction = average_value(vecs) axisvec = make_vec_perp_to_vec( axisvec, strand_direction ) pass assert axisvec is not None self.move_closest_baggage_to( self.posn() + axisvec, [honorary_axis_neighbor] ) # assume it worked (since all candidates were in baggage) pass pass pass pass # (no other cases are possible yet, at least when called by dna updater) # (flag was reset above, so we're done) return # == end of methods for either strand or axis PAM atoms # == PAM Pl atom methods def Pl_preferred_Ss_neighbor(self): # bruce 080118, revised 080401 """ For self a Pl atom (PAM5), return the Ss neighbor it prefers to be grouped with (e.g. in the same chunk, or when one of its bonds is broken) if it has a choice. (If it has no Ss atom, print bug warning and return None.) @warning: the bond direction constant hardcoded into this method is an ARBITRARY GUESS as of 080118. Also it ought to be defined in some dna-related constants file (once this method is moved to a dna-related subclass of Atom). """ assert self.element is Pl5 candidates = [ # note: these might be None, or conceivably a non-Ss atom self.next_atom_in_bond_direction( Pl_STICKY_BOND_DIRECTION), self.next_atom_in_bond_direction( - Pl_STICKY_BOND_DIRECTION) ] candidates = [c for c in candidates if c is not None and \ c.element.symbol.startswith("Ss") # KLUGE, matches Ss3 or Ss5 ] # note: this cond excludes X (good), Pl (bug if happens, but good). # It also excludes Sj and Hp (bad), but is only used from dna updater # so that won't be an issue. Non-kluge variant would test for # "a strand base atom". ## if debug_pref("DNA: Pl5 stick with PAM5 over PAM3 chunk, " \ ## "regardless of bond direction?", ...): if 0: # see if this caused my last bridging Pl bug, bruce 080413 330pm PT: candidates_PAM5 = [c for c in candidates if c.element.pam == MODEL_PAM5] # Try these first, so we prefer Ss5 over Ss3 if both are present, # regardless of bond direction. [bruce 080401] candidates = candidates_PAM5 + candidates ## for candidate in candidates: ## if ...: ## return candidate # all necessary tests were done above -- just return first one, # if any are there if candidates: return candidates[0] print "bug: Pl with no Ss: %r" % self # only a true statement (that this is a bug) when dna updater # is enabled, but that's ok since we're only used then return None def _f_Pl_finish_converting_if_needed(self): """ [friend method for dna updater] Assume self is a Pl5 atom between two Ss3 or Ss5 atoms (or between one such atom and a bondpoint), in the same or different DnaLadders, WHICH MAY NOT HAVE YET REMADE THEIR CHUNKS (i.e. which may differ from atom.molecule.ladder for their atoms, if that even exists), which have finished doing a pam conversion or decided not to or didn't need to (i.e. which are in their proper post-conversion choice of model), and that self has proper directional bonds. Figure out whether self should continue to exist. (This is true iff either neighbor is a PAM5 atom, presumably Ss5.) If not, kill self and directly bond its neighbors, also recording its position on the neighbors which are Ss3 or Ss5 if its position is definitive (if it's not, it means this already happened -- maybe not possible, not sure). (This won't be called twice on self; assert that by asserting self is not killed.) If so, ensure self's position is definitive, which means, if it's not, make it so by setting it from the "+5 data" on self's Ss neighbors intended for self, and removing that data. (This method can be called twice on self in this case; self's position would typically be non-definitive the first time and (due to our side effect that time) definitive the second time.) """ assert not self.killed() sn = self.strand_neighbors() # doesn't include bondpoints pam5_neighbors = [n for n in sn if n.element.pam == MODEL_PAM5] should_exist = not not pam5_neighbors if not should_exist: # make sure it's not due to being called when we have no strand # neighbors (bare Pl) (this method is not legal to call then) assert sn, "error: %r._f_Pl_finish_converting_if_needed() " \ "illegal since no strand neighbors" % self if self._f_Pl_posn_is_definitive: self._f_Pl_store_position_into_Ss3plus5_data() # sets flag to false if debug_flags.DEBUG_DNA_UPDATER_VERBOSE: # 080413 print "fyi: stored pos of %r, "\ "will kill it and rebond neighbors" % self kill_Pl_and_rebond_neighbors(self) ###REVIEW: does killing self mess up its chain or its DnaLadderRailChunk? # (when called from dna updater, those are already invalid, they're # not the new ones we worry about -- UNLESS self was old and untouched # except by this corner Pl. That issue is a predicted unresolved bug # as of 080412 11:54pm PT. #### @@@@) else: if not self._f_Pl_posn_is_definitive: self._f_Pl_set_position_from_Ss3plus5_data() # sets flag to true if debug_flags.DEBUG_DNA_UPDATER_VERBOSE: # 080413 print "fyi: fixed pos of %r, keeping it" % self return def _f_Pl_set_position_from_Ss3plus5_data(self): #bruce 080402 """ [friend method for dna updater] Assume self is a Pl5 atom between two Ss3 or Ss5 atoms, in the same or different DnaLadders, WHICH MAY NOT HAVE YET REMADE THEIR CHUNKS (i.e. which may differ from atom.molecule.ladder for their atoms, if that even exists), and that self has proper directional bonds, and that self's current position is meaningless or out of date (and should be ignored). Set self's position using the "PAM3plus5 Pl-position data" (if any) stored on its neighbors, and remove that data from its neighbors. (If that data is not present, use reasonable defaults.) (Note that the analogous method on class Fake_Pl would *not* remove that data from its neighbors.) Assume that self's neighbors' DnaLadders (which, like self's, may not yet have remade their chunks) have up-to-date stored PAM basepair baseframe data to help do the necessary coordinate conversions. (This might not be checked. Using out of date data would cause hard-to-notice bugs.) No effect on other "PAM3plus5 data" (if any) on those neighbors (e.g. Gv-position data). """ assert not self._f_Pl_posn_is_definitive # note: the following function is also called by class Fake_Pl abspos = Pl_pos_from_neighbor_PAM3plus5_data( self.bond_directions_to_neighbors(), remove_data_from_neighbors = True ) # print "_f_Pl_set_position_from_Ss3plus5_data will set %r " \ # "on %r, now at %r" % (abspos, self, self.posn()) if abspos is None: print "bug: _f_Pl_set_position_from_Ss3plus5_data " \ "can't set %r on %r, now at %r" % \ (abspos, self, self.posn()) else: self.setposn(abspos) del self._f_Pl_posn_is_definitive # like doing self._f_Pl_posn_is_definitive = True, # but expose class default value to save RAM return def _f_Pl_store_position_into_Ss3plus5_data(self): #bruce 080402 """ [friend method for dna updater] Assume self is a Pl5 atom between two Ss3 or Ss5 atoms, in the same or different newly made DnaLadders -- WHICH MAY NOT HAVE YET REMADE THEIR CHUNKS (i.e. which may differ from atom. molecule.ladder for their atoms, if that even exists) -- and that self has proper directional bonds, and that self's current position is definitive (i.e. that any "PAM3plus5 Pl-position data" on self's neighbors should be ignored, and can be entirely replaced). Set the "PAM3plus5 Pl-position data" on self's Ss3 or Ss5 neighbors to correspond to self's current position. (This means that self will soon be killed, with its neighbors rebonded, and transmuted to Ss3 if necessary, but all that is up to the caller.) Assume that self's neighbors' DnaLadders have up-to-date stored PAM basepair baseframe data to help do the necessary coordinate conversions. (This might not be checked. Using out of date data would cause hard-to-notice bugs. Note that the global formerly called ladders_dict only assures that the baseframes are set at the ends.) No effect on other "PAM3plus5 data" (if any) on those neighbors (e.g. Gv-position data). @see: _f_Gv_store_position_into_Ss3plus5_data """ assert self.element is Pl5 assert self._f_Pl_posn_is_definitive # note: unlike with _f_Pl_set_position_from_Ss3plus5_data, # this method has no analogous method used during writemmp-only # PAM conversion -- either we're saving as PAM5 (moving data # onto Pl, not off of it), or as pure PAM3 (no Pl data at all). # PAM3+5 itself has no mmp format. pos = self.posn() for direction_to, ss in self.bond_directions_to_neighbors(): if ss.element.role == 'strand': # (avoid bondpoints or (erroneous) non-PAM or axis atoms) ss._f_store_PAM3plus5_Pl_abs_position( - direction_to, pos) continue self._f_Pl_posn_is_definitive = False return def _f_Gv_store_position_into_Ss3plus5_data(self): # todo: refile within class #bruce 080409, modified from _f_Pl_store_position_into_Ss3plus5_data """ [friend method for dna updater] Assume self is (or very recently was transmuted from??) a Gv5 atom axially between two Ss3 or Ss5 atoms (connected by rung bonds, in the same base pair), in the same DnaLadder (valid and error-free), and that self's current position is definitive (i.e. that any "PAM3plus5 Gv-position data" on self's neighbors should be ignored, and can be entirely replaced), as is true for all live Gv's except *very* transiently when they are created or removed during PAM conversion. Set the "PAM3plus5 Gv-position data" on self's Ss3 or Ss5 neighbors to correspond to self's current position. (This means that self will soon be moved or transmuted (or perhaps just has been??), with its neighbors also transmuted to Ss3 if necessary, but all that is up to the caller.) Assume that self's DnaLadder -- WHICH MAY NOT HAVE YET REMADE ITS CHUNKS (i.e. which may differ from atom.molecule.ladder for its atoms, if that even exists) -- has up-to-date stored PAM basepair baseframe data to help do the necessary coordinate conversions. (This might not be checked. Using out of date data would cause hard-to-notice bugs. Note that the global formerly called ladders_dict only assures that the baseframes are set for the basepairs at the ends of a ladder.) No effect on other "PAM3plus5 data" (if any) on self's neighbors (e.g. Pl-position data). @see: _f_Pl_store_position_into_Ss3plus5_data """ assert self.element.role == 'axis' pos = self.posn() sn = self.strand_neighbors() if len(sn) != 2: print "bug? %r has unexpected number of strand_neighbors %r" % (self, sn) # since ghost bases should be added to bring this up to 2 for ss in sn: assert ss.element.role == 'strand' # (avoid bondpoints or (erroneous) non-PAM or axis atoms) ss._f_store_PAM3plus5_Gv_abs_position( pos) continue return # == end of PAM Pl atom methods # == PAM strand atom methods (some are more specific than that, # e.g. not on Pl or only on Pl, but these are not yet divided up # except that some methods or attrs have Pl or Ss in their names, # and some of those have been moved into the Pl atom methods section above) def _writemmp_PAM3plus5_Pl_Gv_data( self, mapping): #bruce 080523 """ Write the mmp info record (or similar extra data) which represents a non-default value of self._PAM3plus5_Pl_Gv_data. """ vecs = self._PAM3plus5_Pl_Gv_data assert vecs is not None # should be a list of 3 standard "atom position vectors" # (they are relative rather than absolute, but can still # be written in the same manner as atom positions) record = "info atom +5data =" # will be extended below for vec in vecs: if vec is None: vecstring = " ()" # (guessing this is easier to read than None) else: xs, ys, zs = mapping.encode_atom_coordinates( vec ) vecstring = " (%s, %s, %s)" % (xs, ys, zs) # note: 4 of these chars could be left out if we wanted to # optimize the format record += vecstring record += "\n" mapping.write( record) return def _readmmp_3plus5_data(self, key, val, interp): #bruce 080523 """ Read the value part of the mmp info record which represents a non-default value of self._PAM3plus5_Pl_Gv_data. Raise exceptions on errors. @param val: @type val: string """ del key val = val.replace("()", "(None, None, None)") # kluge, for convenience for ignore in "(),": val = val.replace(ignore, " ") # ditto val = val.strip() coord_strings = val.split() assert len(coord_strings) == 9, "wrong length: %r" % (coord_strings,) def decode(coord_string): if coord_string == "None": return None return interp.decode_atom_coordinate(coord_string) coords = map( decode, coord_strings) # each element is a float or None def decode_coords( x, y, z ): if x is None or y is None or z is None: return None return V(x, y, z) x1, y1, z1, x2, y2, z2, x3, y3, z3 = coords d1 = decode_coords(x1, y1, z1) d2 = decode_coords(x2, y2, z2) d3 = decode_coords(x3, y3, z3) self._PAM3plus5_Pl_Gv_data = [d1, d2, d3] return def setDnaBaseName(self, dnaBaseName): # Mark 2007-08-16 #bruce 080319 revised, mainly to support undo/copy """ Set the Dna base letter. This is only valid for PAM atoms in the list VALID_ELEMENTS_FOR_DNABASENAME, i.e. strand sugar atoms, presently defined as ('Ss5', 'Ss3', 'Sh5', 'Se3', 'Sj5', 'Sj3'). @param dnaBaseName: The DNA base letter. This must be "" or one letter (this requirement is only partially enforced here). @type dnaBaseName: str @raise: AssertionError, if self is not a strand sugar atom or if we notice the value is not permitted. """ assert type(dnaBaseName) == type("") #bruce 080326 assert self.element.symbol in VALID_ELEMENTS_FOR_DNABASENAME, \ "Can only assign dnaBaseNames to PAM strand sugar atoms. " \ "Attempting to assign dnaBaseName %r to %r of element %r." \ % (dnaBaseName, self, self.element.name) if dnaBaseName == 'X': dnaBaseName = "" assert len(dnaBaseName) <= 1, \ "dnaBaseName must be empty or a single letter, not %r" % \ (dnaBaseName,) #bruce 080326 # todo: check that it's a letter # maybe: canonicalize the case; if not, make sure it doesn't matter if self._dnaBaseName != dnaBaseName: self._dnaBaseName = dnaBaseName _changed_otherwise_Atoms[self.key] = self #bruce 080319 # todo: in certain display styles (which display dna base letters), # call self.molecule.changeapp(0) self.changed() # bruce 080319 return def getDnaBaseName(self): """ Returns the value of attr I{_dnaBaseName}. @return: The DNA base name, or None if the attr I{_dnaBaseName} does not exist. @rtype: str """ # Note: bruce 080311 revised all direct accesses of atom._dnaBaseName # to go through this method, and renamed it to be private. # (I also stopped storing it in mmp file when value is X, unassigned. # This is desirable to save space and speed up save and load. # If some users of this method want the value on certain atoms # to always exist, this method should be revised to look at # the element type and return 'X' instead of "" for appropriate # elements.) #UPDATE: The following is now revised per above comment. i.e. if it #can't find a baseName for a valid element symbol (see list below) #it makes the dnaBaseName as 'X' (unassigned base). This is useful #while reading in the strand sequence. # See DnaStrand.getStrandSequence() for an example. --Ninad 2008-03-12 valid_element_symbols = VALID_ELEMENTS_FOR_DNABASENAME allowed_on_this_element = (self.element.symbol in valid_element_symbols) baseNameString = self.__dict__.get('_dnaBaseName', "") # modified below if not allowed_on_this_element: #bruce 080319 change: enforce this never existing on disallowed # element (probably just a clarification, since setting code was # probably already enforcing this) baseNameString = "" else: if not baseNameString: baseNameString = 'X' # unassigned base pass if len(baseNameString) > 1: #bruce 080326 precaution, should never happen baseNameString = "" return baseNameString def get_strand_atom_mate(self): """ Returns the 'mate' of this dna pseudo atom (the atom on another strand to which this atom is "base-paired"), or None if it has no mate. @return: B{Atom} (PAM atom) """ #Note: This method was created to support assignment of strand sequence #to strand chunks. This should be moved to dna_model and #can be revised further. -- Ninad 2008-01-14 # (I revised it slightly, to support all kinds of single stranded # regions. -- Bruce 080117) if self.element.role != 'strand': # REVIEW: return None, or raise exception? [bruce 080117 Q] return None #First find the connected axis neighbor axisAtom = self.axis_neighbor() if axisAtom is None: # single stranded region without Ax; no mate return None #Now find the strand atoms connected to this axis atom strandAtoms = axisAtom.strand_neighbors() #... and we want the mate atom of self for atm in strandAtoms: if atm is not self: return atm # if we didn't return above, there is no mate # (single stranded region with Ax) return None def setDnaStrandId_for_generators(self, dnaStrandId_for_generators): # Mark 070904 """ Set the Dna strand name. This is only valid for PAM atoms in the list 'Se3', 'Pe5', 'Pl5' (all deprecated when the dna updater is active). @param dnaStrandId_for_generators: The DNA strand id used by the dna generator @type dnaStrandId_for_generators: str @raise: AssertionError if self is not a Se3 or Pe5 or Pl5 atom. @see: self.getDnaStrandId_for_generators() for more comments. """ # Note: this (and probably its calls) needs revision # for the dna data model. Ultimately its only use will be # to help when reading pre-data-model mmp files. Presently # it's only called when reading "info atom dnaStrandName" # records. Those might be saved on the wrong atomtypes # by current dna updater code, but the caller tolerates # exceptions here (but complains using a history summary). # [bruce 080225/080311 comment] assert self.element.symbol in ('Se3', 'Pe5', 'Pl5'), \ "Can only assign dnaStrandNames to Se3, Pe5, or Pl5 (PAM) atoms. " \ "Attempting to assign dnaStrandName %r to %r of element %r." \ % (dnaStrandId_for_generators, self, self.element.name) # Make sure dnaStrandId_for_generators has all valid characters. #@ Need to allow digits and letters. Mark 2007-09-04 # ## for c in dnaStrandId_for_generators: ## if not c in string.letters: ## assert 0, "Strand id for generatos %r has an invalid " \ ## "character (%r)." % \ ## (dnaStrandId_for_generators, c) self._dnaStrandId_for_generators = dnaStrandId_for_generators def getDnaStrand(self): """ Returns the DnaStrand(group) node to which this atom belongs to. Returns None if there isn't a parent DnaStrand group. @see: Chunk.getDnaStrand() which is used here. """ chunk = self.molecule if chunk and not chunk.isNullChunk(): return chunk.getDnaStrand() return None def getDnaSegment(self): """ Returns the DnaSegment(group) node to which this atom belongs to. Returns None if there isn't a parent DnaSegment group. @see: Chunk.getDnaSegment() which is used here. """ chunk = self.molecule if chunk and not chunk.isNullChunk(): return chunk.getDnaSegment() return None def getDnaStrandId_for_generators(self): """ Returns the value of attr I{_dnaStrandId_for_generators}, or "" if it doesn't exist. @return: The DNA strand id used by dna generator, or "" if the attr I{_dnaStrandId_for_generators} does not exist. @rtype: str """ # Note: this (and probably its calls) need revision for the # dna data model. [bruce 080225/080311 comment] # Note: bruce 080311 revised all direct accesses of # atom._dnaStrandId_for_generators to go through this method, and # renamed it to make it private. #UPDATE renamed previous attr dnastrandName to # dnaStrandId_for_generators based on a discussion with Bruce. #It was renamed to this new name #in order to avoid confusion with the dna strand name which can #be acceesed as node.name. The new name 'dnaStrandId_for_generators' #and this comment makes it clear enough that this will only be used #by generators ... i.e. while creating a duplex from scratch by reading #in the standard mmp files in cad/plugins/PAM*/*.mmp. See #DnaDuplex._regroup to see how this is used. -- Ninad 2008-03-12 return self.__dict__.get('_dnaStrandId_for_generators', "") def directional_bond_chain_status(self): # bruce 071016, revised 080212 """ Return a tuple (statuscode, bond1, bond2) indicating the status of self's bonds with respect to chains of directional bonds. The possible return values are: DIRBOND_CHAIN_MIDDLE, bond1, bond2 -- inside a chain involving these two bonds (note: there might be other directional bonds (open bonds) which should be ignored) DIRBOND_CHAIN_END, bond1, None -- at the end of a chain, which ends with this bond DIRBOND_NONE, None, None -- not in a chain DIRBOND_ERROR, None, None -- local error in directional bond structure, so caller should treat this as not being in a chain Note that all we consider is whether a bond is directional, not whether a direction is actually set (except for atoms with more than one open bond, to disambiguate bare chain ends). Similarly, when two bonds have directions set, we don't consider whether their directions are consistent. (One reason is that we need to grow a chain that covers both bonds so the user can set the entire chain's direction. REVIEW: better to stop at such errors, so only a consistent part of the chain would be changed at once??) But if self is monovalent (e.g. a bondpoint) and its neighbor is not, we consider its neighbor's status in determining its own. Note that when drawing a bond, each of its atoms can have an almost-independent directional_bond_chain_status (due to the possibility of erroneous structures), so both of its atoms need their directional_bond_chain_status checked for errors. """ # note: I think this implem is correct with or without open bonds # being directional [bruce 071016] [revised 080212 to make more true] if not self.element.bonds_can_be_directional: # optimization return DIRBOND_NONE, None, None if len(self.bonds) == 1: # Special cases. This then-clause covers all situations for # self being monovalent, except a few that I think never happen. # (But if they do, fall through to general case below.) bond = self.bonds[0] neighbor = bond.other(self) if len(neighbor.bonds) > 1: # Monovalents defer to non-monovalent neighbors # (note: this applies to bondpoints (after mark 071014 changes) # or to "strand termination atoms".) # Those neighbors may decide bond is not in their chain due # to their other bonds. statuscode, bond1, bond2 = neighbor.directional_bond_chain_status() if statuscode == DIRBOND_NONE or statuscode == DIRBOND_ERROR: return statuscode, None, None elif statuscode == DIRBOND_CHAIN_MIDDLE: # it matters whether we're in the neighbor's chain if bond is bond1 or bond is bond2: return DIRBOND_CHAIN_END, bond, None else: # we're attached to the chain but not in it. # REVIEW: return DIRBOND_ERROR in some cases?? # (For example, when an atom has ._dna_updater__error # set on it?) Note that for open bonds on bare # strands, this happens routinely. return DIRBOND_NONE, None, None # DIRBOND_ERROR? pass elif statuscode == DIRBOND_CHAIN_END: # it matters whether the neighbor's chain includes us. # (for bare strand ends with two open bonds, this is up # to that neighbor even if arbitrary, as of 080212) if bond is bond1: return DIRBOND_CHAIN_END, bond, None else: return DIRBOND_NONE, None, None # DIRBOND_ERROR? pass else: assert 0, "%r got unrecognized statuscode %r from " \ "%r.directional_bond_chain_status" % \ (self, statuscode, neighbor) return DIRBOND_ERROR, None, None pass else: # two connected monovalent atoms, one maybe-directional # bond... for now, proceed with no special case. If this ever # happens, review it. (e.g. we might consider it an error.) pass pass dirbonds = self.directional_bonds() num = len(dirbonds) if num == 2: # it doesn't matter how many of them are open bonds, in this case return DIRBOND_CHAIN_MIDDLE, dirbonds[0], dirbonds[1] elif num == 1: # whether or not it's an open bond return DIRBOND_CHAIN_END, dirbonds[0], None elif num == 0: return DIRBOND_NONE, None, None else: # more than 2 -- see if some of them can be ignored # [behavior at ends of bare strands was revised 080212] real_dirbonds = filter( lambda bond: not bond.is_open_bond(), dirbonds ) num_real = len(real_dirbonds) if num_real == 2: # This works around the situation in which a single strand # (not at the end) has open bonds where axis atoms ought to # be, by ignoring those open bonds. (Note that they count as # directional, even though if they became real they would not # be directional since one atom would be Ax.) # POSSIBLE BUG: the propogate caller can reach this, if it can # start on an ignored open bond. Maybe we should require that # it is not offered in the UI in this case, by having it check # this method before deciding. ### REVIEW return DIRBOND_CHAIN_MIDDLE, real_dirbonds[0], real_dirbonds[1] else: # we need to look at bond directions actually set # (on open bonds anyway), to decide what to do. # # WARNING: this happens routinely at the end of a "bare strand" # (no axis atoms), since it has one real and two open bonds, # all directional. # # We might fix this by: # - making that situation never occur [unlikely] # - making bonds know whether they're directional even if # they're open (bond subclass) # - atom subclass for bondpoints # - notice whether a direction is set on just one open bond # [done below, 080212]; # - construct open bonds on directional elements so the right # number are set [dna updater does that now, but a user bond # dir change can make it false before calling us] # - (or preferably, marked as directional bonds without a # direction being set) # REVIEW: return an error message string? # [bruce 071112, 080212 updated comment] if num_real < 2: # new code (bugfix), bruce 080212 -- look at bonds with # direction set (assume dna updater has made the local # structure make sense) (only works if cmenu won't set dir # on open bond to make 3 dirs set ### FIX) # kluge (bug): assume all real bonds have dir set. # Easily fixable in the lambda if necessary. dirbonds_set = filter( lambda bond: bond._direction, dirbonds ) #k non-private method? if len(dirbonds_set) == 2: return DIRBOND_CHAIN_MIDDLE, dirbonds_set[0], dirbonds_set[1] if debug_flags.atom_debug: print "DIRBOND_ERROR noticed on", self return DIRBOND_ERROR, None, None pass def isThreePrimeEndAtom(self): """ Returns True if self is a three prime end atom of a DnaStrand. (This means the last non-bondpoint atom in the 5'->3' direction.) """ if self.is_singlet(): return False if not self.element.bonds_can_be_directional: return False # optimization nextatom = self.next_atom_in_bond_direction(1) if nextatom is None or nextatom.is_singlet(): return True # self is an end atom return False def isFivePrimeEndAtom(self): """ Returns True if self is a five prime end atom of a DnaStrand. (This means the last non-bondpoint atom in the 3'->5' direction.) """ if self.is_singlet(): return False if not self.element.bonds_can_be_directional: return False # optimization nextatom = self.next_atom_in_bond_direction(-1) if nextatom is None or nextatom.is_singlet(): return True # self is an end atom return False def strand_end_bond(self): #bruce 070415, revised 071016 ### REVIEW: rename? """ For purposes of possibly drawing self as an arrowhead, determine whether self is on the end of a chain of directional bonds (regardless of whether they have an assigned direction). But if self is a bondpoint attached to a chain of directional real bonds, treat it as not part of a bond chain, even if it's directional. [REVIEW: is that wise, if it has a direction set (which is probably an error)?] @return: None, or the sole directional bond on self (if it might be correct to use that for drawing self as an arrowhead). TODO: need a more principled separation of responsibilities between self and caller re "whether it might be correct to draw self as an arrowhead" -- what exactly is it our job to determine? REVIEW: also return an error code, for drawing red arrowheads in the case of certain errors? """ if not self.element.bonds_can_be_directional: return None # optimization statuscode, bond1, bond2 = self.directional_bond_chain_status() if statuscode == DIRBOND_CHAIN_END: assert bond1 assert bond2 is None return bond1 else: return None pass def directional_bonds(self): #bruce 070415 """ Return a list of our directional bonds. Its length might be 0, 1, or 2, or in the case of erroneous structures [or some legal ones as of mark 071014 changes], 3 or more. """ ### REVIEW: should this remain as a separate method, now that its result # can't be used naively? return filter(lambda bond: bond.is_directional(), self.bonds) def bond_directions_are_set_and_consistent(self): #bruce 071204 # REVIEW uses - replace with self._dna_updater__error?? """ Does self (a strand atom, base or base linker) have exactly two bond directions set, not inconsistently? @note: still used, but in some ways superceded by dna updater and the error code it can set. """ count_plus, count_minus = 0, 0 # for all bonds for bond in self.bonds: dir = bond.bond_direction_from(self) if dir == 1: count_plus += 1 elif dir == -1: count_minus += 1 return (count_plus, count_minus) == (1, 1) def desired_new_real_bond_direction(self): #bruce 080213 """ Something is planning to make a new real directional bond involving self, and will modify self's open bonds and/or their bond directions, but not self's existing real bonds or their bond directions. What bond direction (away from self) should it give the new real bond? """ count_plus, count_minus = 0, 0 # for real bonds only for bond in self.bonds: if not bond.is_open_bond(): # (could optim, doesn't matter) dir = bond.bond_direction_from(self) if dir == 1: count_plus += 1 elif dir == -1: count_minus += 1 counts = (count_plus, count_minus) if counts == (1, 0): return -1 elif counts == (0, 1): return 1 else: # Usually or always an error given how we are called, # but let the caller worry about this. # (Non-error cases are conceivable, e.g. within a dna generator, # or if isolated single bases are being bonded by the user.) return 0 pass def fix_open_bond_directions(self, bondpoint, want_dir): #bruce 080213 """ Something is planning to make a new real directional bond involving self, and will give it the bond direction want_dir (measured away from self). If bondpoint is not None, it plans to make this bond using that bondpoint (usually removing it). Otherwise... which bondpoint will it remove, if any?? ###REVIEW CALLERS If possible and advisable, fix all our open bond directions to best fit this situation -- i.e.: * make the direction from self to bondpoint equal want_dir (if bondpoint is not None); * if want_dir is not 0, remove equal directions from other open bonds of self to make the resulting situation consistent (do this even if you have to pick one arbitrarily? yes for now); * if want_dir is 0, do nothing (don't move any direction set on bondpoint to some other open bond of self, because this only happens on error (due to how we are called) and we should do as little as possible here, letting dna updater and/or user see and fix the error). @note: it's ok for this method to be slow. """ debug = debug_flags.atom_debug if debug: print "debug fyi: fix_open_bond_directions%r" % \ ((self, bondpoint, want_dir),) if bondpoint: bondpoint_bond = bondpoint.bonds[0] # redundant: assert bond.other(bondpoint) is self bondpoint_bond.set_bond_direction_from( self, want_dir) else: bondpoint_bond = None # print this until I see whether & how this happens: msg = "not sure what fix_open_bond_directions%r " \ "should do since bondpoint is None" % \ ((self, bondpoint, want_dir),) # not sure, because we might want to deduct want_dir from the # desired total direction we need to achieve on the other # bonds, depending on what caller will do. (If caller will # pick one arbitrarily, we need to know which one that is # now!) if debug_flags.atom_debug: print_compact_stack(msg + ": ") else: print msg + " (set debug flag to see stack)" if not self.bond_directions_are_set_and_consistent(): if debug: print "debug fyi: fix_open_bond_directions%r " \ "needs to fix other open bonds" % \ ((self, bondpoint, want_dir),) # Note: not doing the following would cause undesirable messages, # including history messages from the dna updater, but AFAIK would # cause no other harm when the dna updater is turned on (since the # updater would fix the errors itself if this could have fixed # them). if want_dir: num_fixed = [0] fixable_open_bonds = filter( lambda bond: bond.is_open_bond() and bond is not bondpoint_bond and bond.bond_direction_from(self) != 0 , self.bonds ) def number_with_direction( bonds, dir1 ): """ return the number of bonds in bonds which have the given direction from self """ return len( filter( lambda bond: bond.bond_direction_from(self) == dir1, bonds )) def fix_one( bonds, dir1): "fix one of those bonds (by clearing its direction)" bond_to_fix = filter( lambda bond: bond.bond_direction_from(self) == dir1, bonds )[0] if debug: print "debug fyi: fix_open_bond_directions(%r) " \ "clearing %r of direction %r" % \ (self, bond_to_fix, dir1) bond_to_fix.clear_bond_direction() num_fixed[0] += 1 assert num_fixed[0] <= len(self.bonds) # protect against infloop for dir_to_fix in (1, -1): # or, only fix bonds of direction want_dir? while ( number_with_direction( self.bonds, dir_to_fix ) > 1 and number_with_direction( fixable_open_bonds, dir_to_fix ) > 0 ): fix_one( fixable_open_bonds, dir_to_fix ) continue # if this didn't fix everything, let the dna updater complain # (i.e. we don't need to ourselves) pass return # from fix_open_bond_directions def next_atom_in_bond_direction(self, bond_direction): #bruce 071204 """ Assuming self is in a chain of directional bonds with consistently set directions, return the next atom (of any kind, including bondpoints) in that chain, in the given bond_direction. If the chain does not continue in the given direction, return None. If the assumptions are false, no error is detected, and no exception is raised, but either None or some neighbor atom might be returned. @note: does not verify that bond directions are consistent. Result is not deterministic if two bonds from self have same direction from self (depends on order of self.bonds). """ assert bond_direction in (-1, 1) for bond in self.bonds: dir = bond.bond_direction_from(self) if dir == bond_direction: return bond.other(self) # todo: could assert self is a termination atom or bondpoint, # or if not, that self.bond_directions_are_set_and_consistent() # (if we do, revise docstring) return None def bond_directions_to_neighbors(self): #bruce 080402 """ @return: a list of pairs (bond direction to neighbor, neighbor) for all neighbor atoms to which our bonds are directional (including bondpoints, strand sugar atoms, or Pl atoms). """ res = [] for bond in self.bonds: dir = bond.bond_direction_from(self) if dir: res.append( (dir, bond.other(self)) ) return res def axis_neighbor(self): #bruce 071203; bugfix 080117 for single strand """ Assume self is a PAM strand sugar atom; return the single neighbor of self which is a PAM axis atom, or None if there isn't one (indicating that self is part of a single stranded region). @note: before the dna updater is turned on by default, this may or may not return None for the single-stranded case, since there is no enforcement of one way of representing single strands. After it is turned on, it is likely that it will always return None for free- floating single strands, but this is not fully decided. For "sticky ends" it will return an axis atom, since they will be represented internally as double strands with one strand marked as unreal. """ axis_neighbors = filter( lambda atom: atom.element.role == 'axis', self.neighbors()) if axis_neighbors: assert len(axis_neighbors) == 1, \ "%r.axis_neighbor() finds more than one: %r" % \ (self, axis_neighbors) # stub, since the updater checks needed to ensure this # are NIM as of 071203 return axis_neighbors[0] return None def Pl_neighbors(self): #bruce 080122 """ Assume self is a PAM strand sugar atom; return the neighbors of self which are PAM Pl (pseudo-phosphate) atoms (or any variant thereof, which sometimes interposes between strand base sugar pseudoatoms). """ res = filter( lambda atom: atom.element is Pl5, self.neighbors()) return res def strand_base_neighbors(self): #bruce 071204 (nim, not yet needed; rename?) """ Assume self is a PAM strand sugar atom; return a list of the neighboring PAM strand sugar atoms (even if PAM5 linker atoms separate them from self). """ # review: should the return value also say in which direction each one lies, # whether in terms of bond_direction or base_index_direction? assert 0, "nim" def strand_next_baseatom(self, bond_direction = None): #bruce 071204 """ Assume self is a PAM strand sugar atom, and bond_direction is -1 or 1. Find the next PAM strand sugar atom (i.e. base atom) in the given bond direction, or None if it is missing (since the strand ends), or if any bond directions are unset or inconsistent, or if any other structural error causes difficulty, or if ._dna_updater__error is set in either self or in the atom we might otherwise return (even if that error was propogated from elsewhere in that atom's basepair, rather than being a problem with that atom itself). """ # note: API might be extended to permit passing a baseindex direction # instead, and working on either strand or axis baseatoms. assert bond_direction in (-1, 1) if self._dna_updater__error: #bruce 080131 new feature (part 1 of 3) return None atom1 = self.next_atom_in_bond_direction(bond_direction) # might be None or a bondpoint if atom1 is None: return None if atom1._dna_updater__error: #bruce 080131 new feature (part 2 of 3) return None # REVIEW: the following should probably test element.role == 'strand', # but that includes Se3 and Sh3 end-atoms, unlike this list. # Not an issue when dna updater is on and working, # but if it's disabled to work with an old file, that change # might cause bugs. But I don't feel fully comfortable with # making this depend at runtime on dna_updater_is_enabled() # (not sure why). So leave it alone for now. [bruce 080320] symbol = atom1.element.symbol # KLUGE -- should use another element attr, or maybe Atom subclass if symbol[0:2] not in ('Ss', 'Sj', 'Hp', 'Pl'): # base or base linker atoms (#todo: verify or de-kluge) return None if symbol.startswith('Pl'): # base linker atom # move one more atom to find the one to return atom1 = atom1.next_atom_in_bond_direction(bond_direction) # might be None or a bondpoint assert atom1 is not self # (false would imply one bond had two directions, # or two bonds between same two atoms) if atom1 is None: return None if atom1._dna_updater__error: #bruce 080131 new feature (part 3 of 3) return None if atom1.element.symbol[0:2] not in ('Ss', 'Sj', 'Hp'): return None pass return atom1 def _f_get_fake_Pl(self, direction): #bruce 080327 """ [friend method for PAM3+5 -> PAM5 conversion code] Assume self is a PAM3 or PAM5 strand sugar atom (either PAM model is possible, at least transiently, during PAM3+5 conversion). Find or make, and return, a cached fake Pl5 atom with a properly allocated and unchanging atom.key, to be used in the PAM5 form of self if self does not have a real Pl atom in the given bond_direction. The atom we return might be used only for mmp writing of a converted form, without making any changes in the model, or it might become a live atom and get used in the model. To make it a live atom, special methods must be called [nim] which remove it from being able to be returned by this method. If self does have a real such Pl atom, the atom we might otherwise return might still exist, but this method won't be called to ask for it. It may or may not detect that case. If it detects it, it will treat it as an error. However, it's not an error for the cached atom to survive during the time a live Pl atom takes it place, and to be reused if that live Pl atom ever dies. OTOH, the cached Pl atom might have formerly been a live Pl atom in the same place (i.e. bonded to self in the given bond_direction), killed when self was converted to PAM3. The 3d position (atom.posn()) of the returned atom is arbitrary, and can be changed by the caller for its own purposes. Absent structure changes to self, the identity, key, and 3d position of the returned atom won't be changed between calls of this method by the code that implements the service of which this method is part of the interface. """ fake_Pls = self._fake_Pls # None, or a 2-element list # Note: this list is NOT part of the undoable state, nor are the # fake atoms within it. # # REVIEW: should these have anything to do with storing Pl-posn # "+5" data? is that data in the undoable state? I think it is, # and these are not, so they probably don't help store it. if not fake_Pls: fake_Pls = self._fake_Pls = [None, None] assert direction in (1, -1) direction_index = (direction == 1) # arbitrary map from direction to [0,1] ### clean up dup code when we have some # review: change this to data_index, store a Fake_Gv in the same # array?? (don't know yet if writemmp pam5 conversion will need # one -- maybe not, or if it does it can probably be a flyweight, # but not yet known) Pl = fake_Pls[direction_index] if Pl is None: Pl = fake_Pls[direction_index] = Fake_Pl(self, direction) ## not: self.__class__(Pl5, V(0,0,0)) # obs cmt: maybe: let Pl be live, and if so, verify its bonding with self?? assert isinstance(Pl, Fake_Pl) # nonsense: ## assert Pl.killed() # for now return Pl # methods related to storing PAM3+5 Pl data on Ss # (Gv and Pl data share private helper methods) def _f_store_PAM3plus5_Pl_abs_position(self, direction, abspos, **opts): #bruce 080402, split 080409 """ [friend method for PAM3plus5 code, called from dna updater] If self is a PAM3 or PAM5 strand sugar atom, store "PAM3+5 Pl-position data" on self, for a hypothetical (or actual) Pl in the given bond direction from self, at absolute position abspos, converting this to a relative position using the baseframe info corresponding to self stored in self's DnaLadder, which must exist and be up to date (assumed, not checked). @warning: slow for non-end atoms of very long ladders, due to a linear search for self within the ladder. Could be optimized by passing an index hint if necessary. @note: this method might be inlined, when called by Pls between basepairs in the same DnaLadder, since much of it could be optimized in a loop over basepairs in order (in that case). """ assert direction in (1, -1) data_index = (direction == 1) self._store_PAM3plus5_abspos(data_index, abspos, **opts) # both Pl and Gv use this, with different data_index return def _store_PAM3plus5_abspos(self, data_index, abspos): """ [private helper, used for Pl and Gv methods] """ #bruce 080402, split 080409 if not self.element.symbol.startswith('Ss'): # kluge? needs to match Ss3 and Ss5, and not Pl. # not necessarily an error print "fyi: _store_PAM3plus5_abspos%r is a noop for this element" % \ (self, data_index, abspos) # remove when works if routine; leave in if never seen, to # notice bugs; current caller tries not to call in this case, # so this should not happen return if not self.can_make_up_Pl_abs_position(data_index and 1 or -1): # kluge: can store is the same as can make up, for now; # needed to fix bugs in pam conversion killing Pls next to # single strands [bruce 080413] return origin, rel_to_abs_quat, y_m_junk = _f_baseframe_data_at_baseatom(self) relpos = baseframe_abs_to_rel(origin, rel_to_abs_quat, abspos) if not self._PAM3plus5_Pl_Gv_data: self._PAM3plus5_Pl_Gv_data = [None, None, None] self._PAM3plus5_Pl_Gv_data[data_index] = relpos if debug_flags.DEBUG_DNA_UPDATER_VERBOSE: print "fyi, on %r for data_index %r stored relpos %r" % \ (self, data_index, relpos) # todo: use these prints to get constants for # default_Pl_relative_position (and Gv) [review: done now?] return def _f_recommend_PAM3plus5_Pl_abs_position(self, direction, make_up_position_if_necessary = True, **opts ): """ #doc; return None or a position @warning: this can return None even if make_up_position_if_necessary is passed, since some Ss atoms *can't* make up a position (e.g. those in DnaSingleStrandDomains, at least for now) """ #bruce 080402, split 080409 assert direction in (1, -1) data_index = (direction == 1) res = self._recommend_PAM3plus5_abspos(data_index, **opts) # both Pl and Gv use this, with different data_index if res is None and make_up_position_if_necessary: res = self._make_up_Pl_abs_position(direction) # might be None! # note: don't store this, even if not remove_data # (even though save in PAM5, then reload file, # *will* effectively store it, since in the file we don't # mark it as "made up, should be ignored"; if we did, # we'd want to erase that as soon as anything moved it, # and it'd be a kind of hidden info with mysterious effects, # so it's not a good idea to mark it that way) return res def _recommend_PAM3plus5_abspos(self, data_index, remove_data = False ): """ [private helper, used for Pl and Gv methods] #doc; return None or a position; never make up a position """ #bruce 080402, split 080409 data = None if self._PAM3plus5_Pl_Gv_data: data = self._PAM3plus5_Pl_Gv_data[data_index] # might be None if remove_data: self._PAM3plus5_Pl_Gv_data[data_index] = None # don't bother removing the [None, None, ...] list itself # (optim: should we remove it to conserve RAM??) # data is None or a relative Pl or Gv position if data is None: return None else: relpos = data origin, rel_to_abs_quat, y_m_junk = _f_baseframe_data_at_baseatom(self) return baseframe_rel_to_abs(origin, rel_to_abs_quat, relpos) pass # end of _recommend_PAM3plus5_abspos def _make_up_Pl_abs_position(self, direction): #bruce 080402 """ """ if not self.can_make_up_Pl_abs_position(direction): return None relpos = default_Pl_relative_position(direction) origin, rel_to_abs_quat, y_m_junk = _f_baseframe_data_at_baseatom(self) return baseframe_rel_to_abs(origin, rel_to_abs_quat, relpos) def can_make_up_Pl_abs_position(self, direction): #bruce 080413 try: ladder, whichrail, index = _f_find_new_ladder_location_of_baseatom(self) # not self.molecule.ladder, it finds the old invalid ladder instead except: print "bug? can_make_up_Pl_abs_position can't find new ladder of %r" % \ (self,) return True # should be False, but return True to mitigate new bugs caused # by this feature if not ladder: # probably can't happen print "bug? can_make_up_Pl_abs_position sees ladder of None for %r" % \ (self,) return True # should be False (see comment above) return ladder.can_make_up_Pl_abs_position_for(self, direction) # methods related to storing PAM3+5 Gv data on Ss # (Gv and Pl data share private helper methods) def _f_store_PAM3plus5_Gv_abs_position(self, abspos): #bruce 080409 """ [friend method for PAM3plus5 code, called from dna updater] If self is a PAM3 or PAM5 strand sugar atom, store "PAM3+5 Gv-position data" on self, for a hypothetical (or actual) Gv attached to self, at absolute position abspos, converting this to a relative position using the baseframe info corresponding to self stored in self's DnaLadder, which must exist and be up to date (assumed, not checked), BUT WHICH MAY NOT HAVE YET REMADE ITS CHUNKS (i.e. which may differ from atom. molecule.ladder for its atoms, if that even exists). @warning: slow for non-end atoms of very long ladders, due to a linear search for self within the ladder. Could be optimized by passing an index hint if necessary. @note: this method might be inlined, since much of it could be optimized in a loop over basepairs in order (in that case); or it might be passed a baseindex hint, or baseframe info. """ self._store_PAM3plus5_abspos( _GV5_DATA_INDEX, abspos) # both Pl and Gv use this, with different data_index return def _f_recommend_PAM3plus5_Gv_abs_position(self, make_up_position_if_necessary = True, **opts # e.g. remove_data = False ): """ #doc; return None or a position """ #bruce 080409 res = self._recommend_PAM3plus5_abspos( _GV5_DATA_INDEX, **opts) # both Pl and Gv use this, with different data_index if res is None and make_up_position_if_necessary: res = self._make_up_Gv_abs_position() return res def _make_up_Gv_abs_position(self): #bruce 080409 """ @note: self should be a strand sugar atom """ relpos = default_Gv_relative_position() # (a constant) # following code is the same as for the Pl version origin, rel_to_abs_quat, y_m_junk = _f_baseframe_data_at_baseatom(self) return baseframe_rel_to_abs(origin, rel_to_abs_quat, relpos) # == end of PAM strand atom methods # == PAM axis atom methods def strand_neighbors(self): #bruce 071203 """ Assume self is a PAM axis atom; return the neighbors of self which are PAM strand sugar atoms. There are always exactly one or two of these [NIM] after the dna updater has run. """ # [stub -- need more error checks in following (don't return Pl). # but this is correct if structures have no errors.] res = filter( lambda atom: atom.element.role == 'strand', self.neighbors()) ##assert len(res) in (1, 2), \ ## "error: axis atom %r has %d strand_neighbors (should be 1 or 2)"\ ## % (self, len(res)) # happens in mmkit - leave it as just a print at least until # we implem "delete bare atoms" - legal_nums = legal_numbers_of_strand_neighbors_on_axis() if not ( len(res) in legal_nums ): print "error: axis atom %r has %d strand_neighbors (should be %s)" % \ (self, len(res), " or ".join(map(str, legal_nums))) return res def axis_neighbors(self): #bruce 071204 # (used on axis atoms, not sure if used on strand atoms) return filter( lambda atom: atom.element.role == 'axis', self.neighbors()) # == end of PAM axis atom methods pass # end of class PAM_Atom_methods # end