path: root/cad/src/dna/generators/B_Dna_PAM3_SingleStrand_Generator.py

# Copyright 2008 Nanorex, Inc.  See LICENSE file for details.
"""
B_Dna_PAM3_Generator_SingleStrand allows resizing a single strand of a DNA.
See self.modify(). Note that self.make() to create a DnaStrnad from scratch
is not implemented as of 2008-04-02.

@author: Ninad
@copyright: 2008 Nanorex, Inc.  See LICENSE file for details.
@version:$Id$

History:
Created on 2008-04-02

TODO as of 2008-04-07
- add more documentation. may undergo major refactoring/renaming changes.
- post FNANO allow resizing of a single strand such that: if axis atoms
on original dna match with the ones on the newly created single stranded dna,
delete those atoms and fuse the bare strand with existing atoms
- May be revise the dna creation process thoughout. Example, see how to generate
 positions of strand-axis atoms without using the old way of inserting bases from
 mmp and performing transforms on them.
- move this file to dna/model
- Methods such as the following need REVIEW post FNANO08 (not urgent). May be
 revise algorithm for optimization or supporting some features such as
 'extend and bond
 self._fuse_new_dna_with_original_duplex
 self._replace_overlapping_axisAtoms_of_new_dna

"""

from geometry.VQT import Q, norm, vlen, cross
from Numeric import dot
from utilities.debug import print_compact_stack
from model.bonds import bond_at_singlets

from dna.generators.B_Dna_PAM3_Generator import B_Dna_PAM3_Generator

class B_Dna_PAM3_SingleStrand_Generator(B_Dna_PAM3_Generator):
    """
    B_Dna_PAM3_SingleStrand_Generator allows resizing a single strand of a DNA.
    See self.modify(). Note that self.make() to create a DnaStrnad from scratch
    is not implemented as of 2008-04-02.
    @see: superclass  B_Dna_PAM3_Generator for more documentation,
    """
    def _strandAinfo(self, index):
        """
        Returns parameters needed to add the next base of the
        @param index: Base-pair index.
        @type  index: int
        """
        #unused param index. Just inheriting from the superclass..
        del index

        zoffset      =  0.0
        thetaOffset  =  0.0
        basename     =  "PAM3-SingleStrand/MiddleBasePair"
        basefile     =  self._baseFileName(basename)
        return (basefile, zoffset, thetaOffset)

    def _orient_for_modify(self, end1, end2):
        """
        Do the final orientation of the newly generated dna, around its own
        axis so that the axis and strand atoms on this new dna are
        fusable with the resize end on the original dna.
        @param end1: end1 of new dna (segment) generated
        @param end2: end2 of new dna (segment) generated
        """

        b = norm(end2 - end1)
        new_ladder =   self.axis_atom_end1.molecule.ladder

        chunkListForRotation = new_ladder.all_chunks()

        endBaseAtomList  = new_ladder.get_endBaseAtoms_containing_atom(self.axis_atom_end1)

        endStrandbaseAtoms = (endBaseAtomList[0], endBaseAtomList[2])

        self.strandA_atom_end1 = None
        #TODO: REVIEW / REFACTOR THIS and DOC
        for atm in endStrandbaseAtoms:
            if atm is not None:
                self.strandA_atom_end1 = atm
                self._set_bond_direction_on_new_strand(atm)


        axis_strand_vector = (self.strandA_atom_end1.posn() - \
                              self.axis_atom_end1.posn())


        vectorAlongLadderStep =  self._resizeEndStrand1Atom.posn() - \
                              self._resizeEndAxisAtom.posn()

        unitVectorAlongLadderStep = norm(vectorAlongLadderStep)

        self.final_pos_strand_end_atom = \
            self.axis_atom_end1.posn() + \
            vlen(axis_strand_vector)*unitVectorAlongLadderStep

        q_new = Q(axis_strand_vector, vectorAlongLadderStep)

        if dot(axis_strand_vector, cross(vectorAlongLadderStep, b)) < 0:
            q_new2 = Q(b, -q_new.angle)
        else:
            q_new2 = Q(b, q_new.angle)

        self.assy.rotateSpecifiedMovables(q_new2, chunkListForRotation, end1)
        
        
    def _strand_neighbors_to_delete(self, axisAtom):
        """
        Returns a list of strand neighbors of the given axis atom to delete
        from the original dna being resized (and resizing will result in
        removing bases/ basepairs from the dna). This method determines
        whether both the strand neigbors of this axisAtom need to be deleted
        or is it just a single strand neighbor on a specific Dna Strand Group
        needs to be deleted. The latter is the case while resizing a
        single strand of a Dna. i.e. in this class, it will find a strand neighbor
        to the given axis atom on the strand being resized.
        @see: self._remove_bases_from_dna() where this is called.
        @see: B_Dna_PAM3_Generator._strand_neighbors_to_delete() -- overridden here

        """
        #Note that sicne we are resizing a single strand, the list this
        #method returns will have only one strand atom or will be empty if
        #strand atom on the strand being resized is not found.
        #

        strand_neighbors_to_delete = []

        strand_neighbors  = axisAtom.strand_neighbors()

        if self._resizeEndStrand1Atom:
            mol = self._resizeEndStrand1Atom.molecule
            resize_strand = mol.parent_node_of_class(self.assy.DnaStrand)
            for atm in strand_neighbors:
                strand = atm.molecule.parent_node_of_class(self.assy.DnaStrand)
                if strand is resize_strand:
                    strand_neighbors_to_delete = [atm]
                    break
        return strand_neighbors_to_delete


    def _set_bond_direction_on_new_strand(self,
                                          strandEndAtom_of_new_strand
                                          ):
        """
        Set the bond direction on the new strand. The bond direction will be set
        such that the new strand can be fused with the strand being resized
        Example: if the last strand atom of strand being resized is a 5' end atom,
        this method will make sure that the first strand end atom of the
        new strand we created is a 3' end.
        @param strandEndAtom_of_new_strand: The strand end base atom on the new
               strand. The bond direction will be set from this atom and will
               be propogated further along the strand atom chain this atom
               belongs to (of the new strand)
        @type strandEndAtom_of_new_strand: B{Atom}
        """
        atm = strandEndAtom_of_new_strand

        strand_bond = None
        for bnd in atm.bonds:
            if bnd.is_directional():
                strand_bond = bnd
                break
        if strand_bond:
            mol = self._resizeEndStrand1Atom.molecule

            orig_dnaStrand = mol.parent_node_of_class(self.assy.DnaStrand)

            five_prime_orig, three_prime_orig = \
                           orig_dnaStrand.get_strand_end_base_atoms()

            if self._resizeEndStrand1Atom is five_prime_orig:
                #+1 indicates the absolute bond direction i.e. from 5' to 3' end
                strand_bond.set_bond_direction_from(atm, +1, propogate = True)
            elif self._resizeEndStrand1Atom is three_prime_orig:
                strand_bond.set_bond_direction_from(atm, - 1, propogate = True)

    def _replace_overlapping_axisAtoms_of_new_dna(self, new_endBaseAtomList):
        """
        Checks if the new dna (which is a single strand in this class)
        has any axis atoms that overlap the axis atoms of the original dna.
        If it finds such atoms, the such overlapping atoms of the new dna are
        replaced with that on the original dna. Because of this operation, the
        strand atoms of the new dna are left without axis atoms. So, this method
        then calls appropriate method to create bonds between new strand atoms
        and corresponding axis atoms of the original dna.

        Also, the replacement operation could leave out some neighboring axis
        atoms of the *new dna* without bonds. So those need to be bonded with
        the axis atom of the original dna which replaced their neighbor.
        This is done by calling self._bond_axisNeighbors_with_orig_axisAtoms()

        @see self._fuse_new_dna_with_original_duplex() for  a detail example

        @see: self._bond_bare_strandAtoms_with_orig_axisAtoms()
        @see self._bond_axisNeighbors_with_orig_axisAtoms()

        @BUG: Bug or unsupported feature:
              If the axis atoms of the original dna have a broken bond in
              between, then the wholechain will stop at the broken bond
              and thus, overlapping axis atoms of new dna after that bond
              won't be removed --- in short, the extended strand won't be
              properly fused. We need to come up with well defined rules
              ..example -- when the strand should decide it needs to remove
              overlapping atoms? ...even when it is overlapping with a
              different dna? (not the one we are resizing ) etc.
        """
        #new dna generated by self.modify
        endAxisAtom_new_dna = new_endBaseAtomList[1]
        axis_wholechain_new_dna = endAxisAtom_new_dna.molecule.wholechain
        atomlist_with_overlapping_atoms = axis_wholechain_new_dna.get_all_baseatoms()

        #dna being resized (the original structure which is being resized)
        axis_wholechain_orig_dna = self._resizeEndAxisAtom.molecule.wholechain
        atomlist_to_keep = axis_wholechain_orig_dna.get_all_baseatoms()
        axisEndBaseAtoms_orig_dna = axis_wholechain_orig_dna.end_baseatoms()

        overlapping_atoms = \
                          self._find_overlapping_axisAtomPairs(atomlist_to_keep,
                                                               atomlist_with_overlapping_atoms)


        axis_and_strand_atomPairs_to_bond = []
        axis_and_axis_atomPairs_to_bond = []
        fusableAxisAtomPairsDict = {}

        for atm_to_keep, atm_to_delete in overlapping_atoms:
            #Make sure that atm_to_keep (the axis atom on original dna)
            #has only one strand neighbor, OTHERWISE , if we delete
            #the overlapping axis atoms on the new dna, the bare strand of
            #the new dna can not be bonded with the old dna axis! and we
            #will endup having a bare strand with no axis atom!
            #-Ninad 2008-04-04

            strand_neighbors_of_atm_to_keep = atm_to_keep.strand_neighbors()

            #Before deleting the overlapping axis atom of the new dna,
            #make sure that the corresponding old axis atom has only
            #one strand neighbor. Otherwise, there will be no axis atom
            #available for the bare strand atom that will result because
            #of this delete operation!
            if len(strand_neighbors_of_atm_to_keep) == 1 and atm_to_delete:
                #We know that the new dna is a single strand. So the
                #axis atom will ONLY have a single strand atom atatched.
                #If not, it will be a bug!
                strand_atom_new_dna = atm_to_delete.strand_neighbors()[0]

                #We will fuse this strand atom to the old axis atom
                axis_and_strand_atomPairs_to_bond.append((atm_to_keep,
                                                     strand_atom_new_dna))

                fusableAxisAtomPairsDict[atm_to_keep] = atm_to_delete
                ##fusableAxisAtomPairs.append((atm_to_keep, atm_to_delete))


        for atm_to_keep in axisEndBaseAtoms_orig_dna:
            #This means that we are at the end of the chain of the
            #original dna. There could be some more axis atoms on the
            #new dna that go beyond the original dna axis chain -- [A]
            #Thus, after we replace the overlapping axis atoms of the
            #new dna with the original axis atoms, we should make sure
            #to bond them with the atoms [A] mentioned above
            if fusableAxisAtomPairsDict.has_key(atm_to_keep):
                atm_to_delete = fusableAxisAtomPairsDict[atm_to_keep]
                for neighbor in atm_to_delete.axis_neighbors():
                    if neighbor is not None and \
                       neighbor not in fusableAxisAtomPairsDict.values():
                        axis_and_axis_atomPairs_to_bond.append((atm_to_keep,
                                                           neighbor))

        for atm_to_delete in fusableAxisAtomPairsDict.values():
            try:
                #Now delete the overlapping axis atom on new dna
                atm_to_delete.kill()
            except:
                print_compact_stack("Strand resizing: Error killing axis atom")

        if axis_and_strand_atomPairs_to_bond:
            self._bond_bare_strandAtoms_with_orig_axisAtoms(axis_and_strand_atomPairs_to_bond)

        if axis_and_axis_atomPairs_to_bond:
            self._bond_axisNeighbors_with_orig_axisAtoms(axis_and_axis_atomPairs_to_bond)

    def _bond_bare_strandAtoms_with_orig_axisAtoms(self,
                                                   axis_and_strand_atomPairs_to_bond):
        """
        Create bonds between the bare strand atoms of the new dna with the
        corresponding axis atoms of the original dna. This method should be
        called ONLY from self._replace_overlapping_axisAtoms_of_new_dna() where
        we remove the axis atoms of new dna that overlap the ones in old dna.
        We need to re-bond the bare strand atoms as a result of this replacment,
        with the axis atoms of the old dna.
        @param axis_and_strand_atomPairs_to_bond: A list containing pairs of the
               axis and strand atoms that will be bonded together. It is of the
               format [(atm1, atm2) , ....]
               Where,
               atm1 is always axisAtom of original dna
               atm2 is always bare strandAtom of new dna

        @type  axis_and_strand_atomPairs_to_bond: list
        @see: self._replace_overlapping_axisAtoms_of_new_dna() which calls this.
        """

        #@REVIEW/ TODO: Make sure that we are not bondingan axis atom with a
        #5' or 3' bondpoint of the strand atom.Skip the pair if its one and
        #the same atom (?) This may not be needed because of the other
        #code but is not obvious , so better to make sure in this method
        #-- Ninad 2008-04-11


        for atm1, atm2 in axis_and_strand_atomPairs_to_bond:
            #Skip the pair if its one and the same atom.
            if atm1 is not atm2:
                for s1 in atm1.singNeighbors():
                    if atm2.singNeighbors():
                        s2 = atm2.singNeighbors()[0]
                        bond_at_singlets(s1, s2, move = False)
                        #REVIEW 2008-04-10 if 'break' is needed -
                        break


        #Alternative-- Use find bondable pairs using class fuseChunksBase
        # Fuse the base-pair chunks together into continuous strands.
        #But this DOESN't WORK if bondpoints are disoriented or beyond the
        #tolerance limit. So its always better to bond atoms explicitely as done
        #above
        ##from commands.Fuse.fusechunksMode import fusechunksBase
        ##fcb = fusechunksBase()
        ##fcb.tol = 2.0
        ##fcb.find_bondable_pairs_in_given_atompairs(axis_and_strand_atomPairs_to_bond)
        ##fcb.make_bonds(self.assy)

    def _bond_axisNeighbors_with_orig_axisAtoms(self,
                                                axis_and_axis_atomPairs_to_bond
                                                ):
        """
        The operation that replaces the overlapping axis atoms of the new dnas
        with the corresponding axis atoms of the original dna could leave out
        some neighboring axis atoms of the *new dna* without bonds. So those
        need to be bonded with the axis atom of the original dna which replaced
        their neighbor. This method does that job.
        @param axis_and_axis_atomPairs_to_bond: A list containing pairs of the
               axis atom of original dna and axis atom on the new dna (which was
               a neighbor of the atom deleted in the replacement operation and
               was bonded to it) . Thies eatoms will be bonded with each other.
               It is of the format [(atm1, atm2) , ....]
               Where,
               atm1 is always axisAtom of original dna which replaced the
                     overlapping axis atom of new dna (and thereby broke
                     bond between that atom's axis neighbors )

               atm2 is always axis Atom of new dna , that was a neighbor to
                    an axis atom, say 'A', replaced by original dna axis atom
                    (because 'A' was overlapping) and was previously
                    bonded to 'A'.
        """
        for atm1, atm2 in axis_and_axis_atomPairs_to_bond:
            #Skip the pair if its one and the same atom.
            if atm1 is atm2:
                continue
            #the following check (atm1.singNeighbors()) < 1....)doesn't work in
            #some cases! So explicitly checking the length of the bondpoint
            #neighbors of the atom in the for loop afterwords.
            ##if len(atm1.singNeighbors()) < 1 or len(atm2.singNeighbors()) < 1:
                ##continue

            #Loop through bondpoints of atm1
            for s1 in atm1.singNeighbors():
                #bond point neigbors of atm2
                if atm2.singNeighbors():
                    s2 = atm2.singNeighbors()[0]
                    bond_at_singlets(s1, s2, move = False)
                    break

    def _fuse_new_dna_with_original_duplex(self,
                                           new_endBaseAtomList,
                                           endBaseAtomList):
        """
        First it fuses the end axis and strand atoms of the new dna with the
        corresponding end base atoms of the original dna.

        Then, it checks if the new dna (which is a single strand in this class)
        has any axis atoms that overlap the axis atoms of the original dna.
        If it finds such atoms, the such overlapping atoms of the new dna are
        replaced with that on the original dna. (This is done by calling
        self._replace_overlapping_axisAtoms_of_new_dna())


        EXAMPLE:

        Figure A:
        Figure A shows the original dna -- * denote the axis atoms
        and we will be resizing strand (2)

        (1) 3' <--<--<--<--<--<--<--<--<--<--<--<--< 5'
               |  |  |  |  |  |  |  |  |  |  |  |  |
        (A)    *--*--*--*--*--*--*--*--*--*--*--*--*
               |  |  |  |  |  |  |  |  |  |  |  |  |
        (2) 5' >-->-->-->-->-->3'

        Figure B:
        Lets assume that we extend strand (2) by 4 bases. So,
        the new dna that we will create in self.modify() is indicated by the
        strand (2) with symbols 'X' for strand base atoms and symbol 'o' for
        the axis base atoms. Clearly, the axis atoms 'o' overlap the axis atoms
        of the original dna. So we will remove those. As a result of this
        replacement, strand atoms 'X' will be without any axis atoms , so we
        will bond them with the original axis atoms '*'. The overall result is
        the original duplex's strand 2 gets extended by 4 bases and it still
        remains connected with the corresponding bases in strand (1) (the
        complementary strand )
        (1) 3' <--<--<--<--<--<--<--<--<--<--<--<--< 5'
               |  |  |  |  |  |  |  |  |  |  |  |  |
        (A)    *--*--*--*--*--*--o--o--o--o--*--*--*
               |  |  |  |  |  |  |  |  |  |  |  |  |
        (2) 5' >-->-->-->-->-->==X==X==X==X>

        @see self._replace_overlapping_axisAtoms_of_new_dna() for more details

        TODO:
         - method needs to be renamed The original dna may be a single stranded
        dna or a duplex. Until 2008-04-02 ,it was possible to create or modify
        only a duplex and thats why the name 'duplex'
        -  Refactor further. May be use superclass method with some mods.

        @see: self.modify()
        @see: B_Dna_PAM3_Generator._fuse_new_dna_with_original_duplex()


        """

        #FUSE new dna strand with the original duplex
        chunkList1 = \
                   [ new_endBaseAtomList[0].molecule,
                     self._resizeEndStrand1Atom.molecule]

        chunkList2 = \
                   [ new_endBaseAtomList[1].molecule,
                      self._resizeEndAxisAtom.molecule]


        #Set the chunk color and chunk display of the new duplex such that
        #it matches with the original duplex chunk color and display
        #Actually, fusing the chunks should have taken care of this, but
        #for some unknown reasons, its not happening. May be because
        #chunks are not 'merged'?  ... Setting display and color for new
        #duplex chunk is explicitely done below. Fixes bug 2711
        for chunkPair in (chunkList1, chunkList2):
            display = chunkPair[1].display
            color   = chunkPair[1].color
            chunkPair[0].setDisplayStyle(display)
            if color:
                chunkPair[0].setcolor(color)

            #Original implementation which relied on  on fuse chunks for finding
            #bondable atom pairs within a tolerance limit. This is no longer
            #used and can be removed after more testing of explicit bonding
            #done in self._bond_atoms_in_atomPairs() (called below)
            #-- Ninad 2008-04-11
            ##self.fuseBasePairChunks(chunkPair)

        endAtomPairsToBond = [ (new_endBaseAtomList[0],
                                self._resizeEndStrand1Atom),

                                (new_endBaseAtomList[1],
                                 self._resizeEndAxisAtom)]

        #Create explicit bonds between the end base atoms
        #(like done in self._bond_bare_strandAtoms_with_orig_axisAtoms())
        #instead of relying on fuse chunks (which relies on finding
        #bondable atom pairs within a tolerance limit. This fixes bug 2798
        self._bond_atoms_in_atomPairs(endAtomPairsToBond)
        
        

        #Make sure to call this AFTER the endatoms of new and
        #original DNAs are joined. Otherwise (if you replace the overlapping
        #atoms first), if the endAxisAtom of the new dna is an overlapping
        #atom, it will get deleted and in turn it will invalidate the
        #new_endBaseAtomList! (and thus, even the end strand atom won't
        #be fused with original dna end strand atom)
        self._replace_overlapping_axisAtoms_of_new_dna(new_endBaseAtomList)
        


    def _find_overlapping_axisAtomPairs(self,
                                         atomlist_to_keep,
                                         atomlist_with_overlapping_atoms,
                                         tolerance = 1.1
                                         ):
        """
        @param atomlist_to_keep: Atomlist which will be used as a reference atom
               list. The atoms in this list will be used to find the atoms
               in the *other list* which overlap atom positions in *this list*.
               Thus, the atoms in 'atomlist_to_keep' will be preserved (and thus
               won't be appended to self.overlapping_atoms)
        @type atomlist_to_keep: list
        @param atomlist_with_overlapping_atoms: This list will be checked with
               the first list (atom_list_to_keep) to find overlapping atoms.
               The atoms in this list that overlap with the atoms from the
               original list will be appended to self.overlapping_atoms
               (and will be eventually deleted)

        """

        overlapping_atoms_to_delete = []

        # Remember: chunk = a selected chunk  = atomlist to keep
        # mol = a non-selected -- to find overlapping atoms from


        # Loop through all the atoms in the selected chunk.
        # Use values() if the loop ever modifies chunk or mol--
        for a1 in atomlist_to_keep:
            # Singlets can't be overlapping atoms. SKIP those
            if a1.is_singlet():
                continue

            # Loop through all the atoms in atomlist_with_overlapping_atoms.
            for a2 in atomlist_with_overlapping_atoms:
                # Only atoms of the same type can be overlapping.
                # This also screens singlets, since a1 can't be a
                # singlet.
                if a1.element is not a2.element:
                    continue

                # Compares the distance between a1 and a2.
                # If the distance
                # is <= tol, then we have an overlapping atom.
                # I know this isn't a proper use of tol,
                # but it works for now.   Mark 050901
                if vlen (a1.posn() - a2.posn()) <= tolerance:
                    # Add this pair to the list--
                    overlapping_atoms_to_delete.append( (a1,a2) )
                    # No need to check other atoms in this chunk--
                    break

        return overlapping_atoms_to_delete