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

# Copyright 2004-2008 Nanorex, Inc.  See LICENSE file for details. 
"""
DnaDuplex.py -- DNA duplex generator helper classes, based on empirical data.

@author: Mark Sims
@version: $Id$
@copyright: 2004-2008 Nanorex, Inc.  See LICENSE file for details.

History:

Mark 2007-10-18:
- Created. Major rewrite of DnaGenHelper.py.
"""

import foundation.env as env
from utilities.debug import print_compact_stack
from platform_dependent.PlatformDependent import find_plugin_dir
from geometry.VQT import Q, cross, vlen
from utilities.Log      import orangemsg
from utilities.exception_classes import PluginBug
from simulation.sim_commandruns import adjustSinglet
from model.elements import PeriodicTable

Element_Ae3 = PeriodicTable.getElement('Ae3')
basepath_ok, basepath = find_plugin_dir("DNA")
if not basepath_ok:
    env.history.message(orangemsg("The cad/plugins/DNA directory is missing."))

RIGHT_HANDED = -1
LEFT_HANDED  =  1

from geometry.VQT import norm
from Numeric import dot
from dna.generators.B_Dna_Generator import B_Dna_Generator

class B_Dna_PAM3_Generator(B_Dna_Generator):
    """
    Provides a PAM3 reduced model of the B form of DNA.
    """
    model = "PAM3"

    strandA_atom_end1 = None

    def _strandAinfo(self, index):
        """
        Returns parameters needed to add the next base-pair to the duplex 
        build built.

        @param index: Base-pair index.
        @type  index: int
        """

        zoffset      =  0.0
        thetaOffset  =  0.0
        basename     =  "MiddleBasePair"
        basefile     =  self._baseFileName(basename)
        return (basefile, zoffset, thetaOffset)




    def _postProcess(self, baseList):
        """
        Final tweaks on the DNA chunk, including:

          - Transmute Ax3 atoms on each end into Ae3.
          - Adjust open bond singlets.

        @param baseList: List of basepair chunks that make up the duplex.
        @type  baseList: list

        @note: baseList must contain at least two base-pair chunks.
        """

        if len(baseList) < 1:
            print_compact_stack("bug? (ignoring) DnaDuplex._postProcess called but "\
                                "baseList is empty. Maybe dna_updater was "\
                                "run?: ")
            return

        start_basepair_atoms = baseList[0].atoms.values()
        end_basepair_atoms = baseList[-1].atoms.values()

        Ax_caps = filter( lambda atom: atom.element.symbol in ('Ax3'), 
                          start_basepair_atoms)
        Ax_caps += filter( lambda atom: atom.element.symbol in ('Ax3'), 
                           end_basepair_atoms)

        # Transmute Ax3 caps to Ae3 atoms. 
        # Note: this leaves two "killed singlets" hanging around,  
        #       one from each Ax3 cap.
        #
        # REVIEW: is it safe to simply not do this when dna_updater_is_enabled()?
        # [bruce 080320 question]
        for atom in Ax_caps:
            atom.Transmute(Element_Ae3)

        # X_List will contain 6 singlets, 2 of which are killed (non-bonded).
        # The other 4 are the 2 pair of strand open bond singlets.
        X_List = filter( lambda atom: atom.element.symbol in ('X'), 
                         start_basepair_atoms)
        X_List += filter( lambda atom: atom.element.symbol in ('X'), 
                          end_basepair_atoms)

        # Adjust the 4 open bond singlets.
        for singlet in X_List:
            if singlet.killed():
                # Skip the 2 killed singlets.
                continue
            adjustSinglet(singlet)


        return



    def _determine_axis_and_strandA_endAtoms_at_end_1(self, chunk):
        """
        Determine the axis end atom and the strand atom on strand 1 
        connected to it , at end1 . i.e. the first mouse click point from which 
        user started drawing the dna duplex (rubberband line)

        These are initially set to None. 

        The strand base atom of Strand-A, connected to the self.axis_atom_end1
        The vector between self.axis_atom_end1 and self.strandA_atom_end1 
        is used to determine the final orientation of the created duplex
        done in self._orient_to_position_first_strandA_base_in_axis_plane()

        This vector is aligned such that it is parallel to the screen. 

        i.e. StrandA end Atom and corresponding axis end atom are coplaner and 
        parallel to the screen.   

        @NOTE:The caller should make sure that the appropriate chunk is passed 
              as an argument to this method. This function itself only finds 
              and assigns the axis ('Ax3') and strand ('Ss3' atom and on StrandA)
              atoms it sees first to the  respective attributes. (and which pass
              other necessary tests)        

        @see: self.make() where this function is called
        @see: self._orient_to_position_first_strandA_base_in_axis_plane()
        """
        for atm in chunk.atoms.itervalues():
            if self.axis_atom_end1 is None:
                if atm.element.symbol == 'Ax3':
                    self.axis_atom_end1 = atm
            if self.strandA_atom_end1 is None:
                if atm.element.symbol == 'Ss3' and atm.getDnaBaseName() == 'a':
                    self.strandA_atom_end1 = atm

    def _orient_for_modify(self, end1, end2):
        """
        Orient the new dna to match up appropriately with the original dna 
        (being modified/resized)
        """

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

        chunkListForRotation = new_ladder.all_chunks()

        #This method fixes bug 2889. See that method for more comments.    
        self._redetermine_resizeEndStrand1Atom_and_strandA_atom_end1()

        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 _redetermine_resizeEndStrand1Atom_and_strandA_atom_end1(self):
        """
        @ATTENTION: The strandA endatom at end1 is modified in this method.
        It is originally computed in self._determine_axis_and_strandA_endAtoms()

        The recomputation is done to fix bug 2889 (for v1.1.0).  
        See B_Dna_PAM3_Generator.orient_for_modify() for details. 
        This NEEDS CLEANUP 
        @see: self._create_raw_duplex()
        @see: self.orient_for_modify()
        """
        #Method created just before codefreeze for v1.1.0 to fix newly discovered
        #bug 2889 -- Ninad 2008-06-02
        #Perhaps computing this strand atom  always be done at a later 
        #stage (like done in here). But not sure if this will cause any bugs. 
        #So not changing the original implementation .


        new_ladder =   self.axis_atom_end1.molecule.ladder       
        endBaseAtomList  = new_ladder.get_endBaseAtoms_containing_atom(self.axis_atom_end1)

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

        self.strandA_atom_end1 = None

        #Check if the resizeEndStrandAtom is a 5' or a 3' end. 
        #If it is a 5' or 3' end, then chose the strand end atom of the 
        #new duplex such that it is the opposite of it (i.e. if resizeEndStrandAtom 
        #is a 5' end, the strand end atom on new duplex should be chosen such that
        #its a 3' end. Using these references, we will orient the new duplex 
        #to match up correctly with the resizeEnd strand atom (and later
        #the old and new dnas will be bonded at these ends)

        #However, if the chosen resizeEndStrandAtom of original duplex 
        #doesn't have a 5' or 3' end, then we will choose the second 
        #resizeEndStrandEndAtom on the original duplex and do the same 
        #computations
        resizeEndStrandAtom_isFivePrimeEndAtom = False
        resizeEndStrandAtom_isThreePrimeEndAtom = False        
        if self._resizeEndStrand1Atom.isFivePrimeEndAtom():            
            resizeEndStrandAtom_isFivePrimeEndAtom = True
        elif self._resizeEndStrand1Atom.isThreePrimeEndAtom():
            resizeEndStrandAtom_isThreePrimeEndAtom = True

        if not (resizeEndStrandAtom_isFivePrimeEndAtom or \
                resizeEndStrandAtom_isThreePrimeEndAtom):
            if self._resizeEndStrand2Atom:
                if self._resizeEndStrand2Atom.isFivePrimeEndAtom():            
                    resizeEndStrandAtom_isFivePrimeEndAtom = True
                elif self._resizeEndStrand2Atom.isThreePrimeEndAtom():
                    resizeEndStrandAtom_isThreePrimeEndAtom = True

                if (resizeEndStrandAtom_isFivePrimeEndAtom or \
                    resizeEndStrandAtom_isThreePrimeEndAtom):
                    #Swap resizeEndStrand1Atom and resizeEndStrand2Atom
                    atm1, atm2 = self._resizeEndStrand1Atom, self._resizeEndStrand2Atom
                    self._resizeEndStrand1Atom, self._resizeEndStrand2Atom = atm2, atm1


        for atm in endStrandbaseAtoms:
            if atm is not None:
                if atm.isThreePrimeEndAtom() and \
                   resizeEndStrandAtom_isFivePrimeEndAtom:
                    self.strandA_atom_end1 = atm 
                    break 
                elif atm.isFivePrimeEndAtom() and \
                     resizeEndStrandAtom_isThreePrimeEndAtom:
                    self.strandA_atom_end1 = atm 
                    break 

        if self.strandA_atom_end1 is None:
            #As a fallback, set this atom to any atom in endStrandbaseAtoms
            #but, this may cause a bug in which bond directions are not 
            #properly set.
            for atm in endStrandbaseAtoms:
                if atm is not None:
                    self.strandA_atom_end1 = atm




    def _orient_to_position_first_strandA_base_in_axis_plane(self, baseList, end1, end2):
        """
        The self._orient method orients the DNA duplex parallel to the screen
        (lengthwise) but it doesn't ensure align the vector
        through the strand end atom on StrandA and the corresponding axis end 
        atom  (at end1) , parallel to the screen. 

        This function does that ( it has some rare bugs which trigger where it
        doesn't do its job but overall works okay )

        What it does: After self._orient() is done orienting, it finds a Quat 
        that rotates between the 'desired vector' between strand and axis ends at
        end1(aligned to the screen)  and the actual vector based on the current
        positions of these atoms.  Using this quat we rotate all the chunks 
        (as a unit) around a common center.

        @BUG: The last part 'rotating as a unit' uses a readymade method in 
        ops_motion.py -- 'rotateSpecifiedMovables' . This method itself may have
        some bugs because the axis of the dna duplex is slightly offset to the
        original axis. 

        @see: self._determine_axis_and_strandA_endAtoms_at_end_1()
        @see: self.make()
        """

        #the vector between the two end points. these are more likely
        #points clicked by the user while creating dna duplex using endpoints
        #of a line. In genral, end1 and end2 are obtained from self.make()
        b = norm(end2 - end1)        

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

        vectorAlongLadderStep =  cross(-self.assy.o.lineOfSight, b)
        unitVectorAlongLadderStep = norm(vectorAlongLadderStep)

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

        expected_vec = self.final_pos_strand_end_atom - self.axis_atom_end1.posn()

        q_new = Q(axis_strand_vector, expected_vec)

        if dot(axis_strand_vector, self.assy.o.lineOfSight) < 0:            
            q_new2 = Q(b, -q_new.angle)
        else:     
            q_new2 = Q(b, q_new.angle) 


        self.assy.rotateSpecifiedMovables(q_new2, baseList, 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 ladder 
        needs to be deleted. The latter is the case while resizing a 
        single strand of a Dna. 
        @see: self._remove_bases_from_dna() where this is called. 
        @see: Dna._strand_neighbors_to_delete() -- overridden here
        @see: B_Dna_PAM3_Generator_SingleStrand._strand_neighbors_to_delete() which 
              overrides this method
        """
        strand_neighbors_to_delete = axisAtom.strand_neighbors()
        return strand_neighbors_to_delete


    def _create_atomLists_for_regrouping(self, dnaGroup):
        """
        Creates and returns the atom lists that will be used to regroup the 
        chunks  within the DnaGroup. 

        @param dnaGroup: The DnaGroup whose atoms will be filtered and put into 
                         individual strand A or strandB or axis atom lists.
        @return: Returns a tuple containing three atom lists 
                 -- two atom lists for strand chunks and one for axis chunk. 
        @see: self._regroup()
        """
        _strandA_list  =  []
        _strandB_list  =  []
        _axis_list     =  []

        # Build strand and chunk atom lists.
        for m in dnaGroup.members:
            for atom in m.atoms.values():
                if atom.element.symbol in ('Ss3'):
                    if atom.getDnaBaseName() == 'a':
                        _strandA_list.append(atom)
                        #Now reset the DnaBaseName for the added atom 
                        # to 'unassigned' base i.e. 'X'
                        atom.setDnaBaseName('X')
                    elif atom.getDnaBaseName() == 'b':
                        _strandB_list.append(atom)
                        #Now reset the DnaBaseName for the added atom 
                        # to 'unassigned' base i.e. 'X'
                        atom.setDnaBaseName('X')
                    else:
                        msg = "Ss3 atom not assigned to strand 'a' or 'b'."
                        raise PluginBug(msg)
                elif atom.element.symbol in ('Ax3', 'Ae3'):
                    _axis_list.append(atom)

        return (_strandA_list, _strandB_list, _axis_list)