path: root/cad/src/outtakes/BuildDuplex_old/DnaGenHelper.py

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

WARNING: this file has been mostly superseded by DnaDuplex.py.
It is used only in DnaGenerator.py, also superseded.

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

History:

Jeff 2007-06-13 (created)
  - Transfered class Dna, A_Dna, A_Dna_PAM5, B_Dna, B_Dna_PAM5,
    Z_Dna abd Z_Dna_PAM5 from Will Ware's DNAGenerator.py.
  - Added accessor methods getBaseRise/setBaseRise

Mark 2007-08-17:
  - make() creates a duplex structure only (no single strands).
  - Created standard directory structure for DNA MMP files.
  - PAM5 MMP files reworked and simplfied.
  - A single PAM5 base-pair duplex is made correctly.
  - Moved Dna constants to Dna_Constants.py.
  - Changed default base to "X" (unassigned) from "N" (aNy base).
  - Implemented Standard IUB Codes (see Dna_Constants.py)

Mark 2007-10-18:
- Did a major rewrite of this module, superseding it --
  DnaDuplex.py.

Bruce 2008-1-1:
- renamed this from Dna.py to DnaGenHelper.py so it's
  not in a way of a proposed "dna" package.
"""

# To do:
#
# 1) Atomistic and PAM5 generated models should have exact orientation
# (i.e. rotational origin).

import foundation.env as env
import os
import re

from math    import sin, cos, pi

from utilities.debug import print_compact_traceback

from platform_dependent.PlatformDependent import find_plugin_dir
from files.mmp.files_mmp import readmmp
from geometry.VQT import V
from commands.Fuse.fusechunksMode import fusechunksBase
from utilities.Log import orangemsg
from utilities.exception_classes import PluginBug

from dna.model.Dna_Constants import basesDict, dnaDict

atompat = re.compile("atom (\d+) \((\d+)\) \((-?\d+), (-?\d+), (-?\d+)\)")
numberPattern = re.compile(r"^\s*(\d+)\s*$")

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

class Dna:
    """
    Dna base class. It is inherited by B_Dna and Z_Dna subclasses.

    @ivar baseRise: The rise (spacing) between base-pairs along the helical
                    (Z) axis.
    @type baseRise: float

    @ivar handedness: Right-handed (B and A forms) or left-handed (Z form).
    @type handedness: int

    @ivar model: The model representation, where:
                          - "Atomistic" = atomistic model
                          - "PAM5" = PAM5 reduced model.
    @type model: str

    @ivar sequence: The sequence (of strand A of the duplex).
    @type sequence: str
    """

    def make(self, group, inSequence, basesPerTurn, position = V(0, 0, 0)):
        """
        Makes a DNA duplex with the sequence I{inSequence}.
        The duplex is oriented with its central axis coincident
        to the Z axis.

        @param assy: The assembly (part).
        @type  assy: L{assembly}

        @param group: The group node object containing the DNA. The caller
                      is responsible for creating an empty group and passing
                      it here. When finished, this group will contain the DNA
                      model.
        @type  group: L{Group}

        @param basesPerTurn: The number of bases per helical turn.
        @type  basesPerTurn: float

        @param position: The position in 3d model space at which to create
                         the DNA strand. This should always be 0, 0, 0.
        @type position:  position
        """

        self.sequence  =  inSequence
        assy           =  group.assy
        baseList       =  []

        def insertBaseFromMmp(filename, subgroup, tfm, position = position):
            """
            Insert the atoms for a nucleic acid base from an MMP file into
            a single chunk.
             - If atomistic, the atoms for each base are in a separate chunk.
             - If PAM5, the pseudo atoms for each base-pair are together in a
               chunk.

            @param filename: The mmp filename containing the base
                             (or base-pair).
            @type  filename: str

            @param subgroup: The part group to add the atoms to.
            @type  subgroup: L{Group}

            @param tfm: Transform applied to all new base atoms.
            @type  tfm: V

            @param position: The origin in space of the DNA duplex, where the
                             3' end of strand A is 0, 0, 0.
            @type  position: L{V}
            """
            try:
                ok, grouplist = readmmp(assy, filename, isInsert = True)
            except IOError:
                raise PluginBug("Cannot read file: " + filename)
            if not grouplist:
                raise PluginBug("No atoms in DNA base? " + filename)

            viewdata, mainpart, shelf = grouplist

            for member in mainpart.members:
                # 'member' is a chunk containing a set of base-pair atoms.
                for atm in member.atoms.values():
                    atm._posn = tfm(atm._posn) + position
                    if atm.element.symbol in ('Se3', 'Ss3', 'Ss5'):
                        if atm.getDnaBaseName() == "a":
                            baseLetter = currentBaseLetter
                        else:
                            try:
                                baseLetter = \
                                    basesDict[currentBaseLetter]['Complement']
                            except:
                                # If complement not found, just assign same
                                # letter.
                                baseLetter = currentBaseLetter
                        if 0:
                            print "Ss(%r) being set to %r." \
                                  % (atm.getDnaBaseName(), baseLetter)
                        atm.setDnaBaseName(baseLetter)

                member.name = currentBaseLetter
                subgroup.addchild(member)
                baseList.append(member)

            # Clean up.
            del viewdata
            shelf.kill()

        def rotateTranslateXYZ(inXYZ, theta, z):
            """
            Returns the new XYZ coordinate rotated by I{theta} and
            translated by I{z}.

            @param inXYZ: The original XYZ coordinate.
            @type  inXYZ: V

            @param theta: The base twist angle.
            @type  theta: float

            @param z: The base rise.
            @type  z: float

            @return: The new XYZ coordinate.
            @rtype:  V
            """
            c, s = cos(theta), sin(theta)
            x = c * inXYZ[0] + s * inXYZ[1]
            y = -s * inXYZ[0] + c * inXYZ[1]
            return V(x, y, inXYZ[2] + z)

        # Begin making DNA.
        subgroup = group
        subgroup.open = False

        # Calculate the twist per base in radians.
        twistPerBase = (self.handedness * 2 * pi) / basesPerTurn
        theta        = 0.0
        z            = 0.5 * self.baseRise * (len(self.sequence) - 1)

        # Create PAM5 or PAM3 duplex.
        for i in range(len(self.sequence)):
            currentBaseLetter = self.sequence[i]
            basefile, zoffset, thetaOffset = \
                self._strandAinfo(currentBaseLetter, i)
            def tfm(v, theta = theta + thetaOffset, z1 = z + zoffset):
                return rotateTranslateXYZ(v, theta, z1)

            insertBaseFromMmp(basefile, subgroup, tfm)

            theta -= twistPerBase
            z     -= self.baseRise

        # Fuse the base chunks together into continuous strands.
        fcb = fusechunksBase()
        fcb.tol = 1.5
        for i in range(len(baseList) - 1):
            fcb.find_bondable_pairs([baseList[i]], [baseList[i + 1]])
            fcb.make_bonds(assy)

        try:
            self._postProcess(baseList)
        except:
            if env.debug():
                print_compact_traceback(
                    "debug: exception in %r._postProcess(baseList = %r) \
                    (reraising): " % (self, baseList,))
            raise

        return

    def _postProcess(self, baseList):
        return

    def _baseFileName(self, basename):
        """
        Returns the full pathname to the mmp file containing the atoms
        of a nucleic acid base (or base-pair).

        Example: If I{basename} is "MidBasePair" and this is a PAM5 model of
        B-DNA, this returns:

          - "C:$HOME\cad\plugins\DNA\B-DNA\PAM5-bases\MidBasePair.mmp"

        @param basename: The basename of the mmp file without the extention
                         (i.e. "adenine", "MidBasePair", etc.).
        @type  basename: str

        @return: The full pathname to the mmp file.
        @rtype:  str
        """
        form    = self.form             # A-DNA, B-DNA or Z-DNA
        model   = self.model + '-bases' # Atomistic or PAM5
        return os.path.join(basepath, form, model, '%s.mmp' % basename)

    def getBaseRise( self ):
        """
        Get the base rise (spacing) between base-pairs.
        """
        return float( self.baseRise )

    def setBaseRise( self, inBaseRise ):
        """
        Set the base rise (spacing) between base-pairs.

        @param inBaseRise: The base rise in Angstroms.
        @type  inBaseRise: float
        """
        self.baseRise  =  inBaseRise

    pass

class A_Dna(Dna):
    """
    Provides an atomistic model of the A form of DNA.

    The geometry for A-DNA is very twisty and funky. We need to to research
    the A form since it's not a simple helix (like B) or an alternating helix
    (like Z).

    @attention: This class is not implemented yet.
    """
    form       =  "A-DNA"
    model      =  "Atomistic"
    baseRise   =  dnaDict['A-DNA']['DuplexRise']
    handedness =  RIGHT_HANDED

    def _strandAinfo(self, baseLetter, index):
        """
        Raise exception since A-DNA is not support.
        """
        raise PluginBug("A-DNA is not yet implemented.")

    def _strandBinfo(self, sequence, index):
        """
        Raise exception since A-DNA is not support.
        """
        raise PluginBug("A-DNA is not yet implemented.")
    pass

class A_Dna_Atomistic(A_Dna):
    """
    Provides an atomistic model of the A form of DNA.

    @attention: This class is not implemented yet.
    """
    pass

class A_Dna_PAM5(A_Dna):
    """
    Provides a PAM5 reduced model of the B form of DNA.

    @attention: This class is not implemented yet.
    """
    pass

class B_Dna(Dna):
    """
    Provides an atomistic model of the B form of DNA.
    """
    form       =  "B-DNA"
    baseRise   =  dnaDict['B-DNA']['DuplexRise']
    handedness =  RIGHT_HANDED
    pass

class B_Dna_Atomistic(B_Dna):
    """
    Provides an atomistic model of the B form of DNA.
    """
    model      =  "Atomistic"

    def _strandAinfo(self, baseLetter, index):
        """
        Returns parameters needed to add a base to strand A.

        @param baseLetter: The base letter.
        @type  baseLetter: str

        @param index: Index in base sequence. This is unused.
        @type  index: int
        """
        zoffset      =  0.0
        thetaOffset  =  0.0
        basename     =  basesDict[baseLetter]['Name']
        basefile     =  self._baseFileName(basename)
        return (basefile, zoffset, thetaOffset)

    def _strandBinfo(self, baseLetter, index):
        """
        Returns parameters needed to add a base to strand B.

        @param baseLetter: The base letter.
        @type  baseLetter: str

        @param index: Index in base sequence. This is unused.
        @type  index: int
        """
        zoffset      =  0.0
        thetaOffset  =  210 * (pi / 180)
        basename     =  basesDict[baseLetter]['Name']
        basefile     =  self._baseFileName(basename)
        return (basefile, zoffset, thetaOffset)

    pass

class B_Dna_PAM5(B_Dna):
    """
    Provides a PAM5 reduced model of the B form of DNA.
    """
    model      =  "PAM5"

    def _isStartPosition(self, index):
        """
        Returns True if I{index} points the first base position (5') of
        self's sequence.

        @param index: Index in base sequence.
        @type  index: int

        @return: True if index is zero.
        @rtype : bool
        """
        if index == 0:
            return True
        else:
            return False

    def _isEndPosition(self, index):
        """
        Returns True if I{index} points the last base position (3') of self's
        sequence.

        @param index: Index in base sequence.
        @type  index: int

        @return: True if index is zero.
        @rtype : bool
        """
        if index ==  len(self.sequence) - 1:
            return True
        else:
            return False

    def _strandAinfo(self, baseLetter, index):
        """
        Returns parameters needed to add a base to strand A.

        @param baseLetter: The base letter.
        @type  baseLetter: str

        @param index: Index in base sequence.
        @type  index: int
        """
        zoffset      =  0.0
        thetaOffset  =  0.0
        basename     =  "MiddleBasePair"

        if self._isStartPosition(index):
            basename = "StartBasePair"

        if self._isEndPosition(index):
            basename = "EndBasePair"

        if len(self.sequence) == 1:
            basename = "SingleBasePair"

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

    def _postProcess(self, baseList): # bruce 070414
        """
        Set bond direction on the backbone bonds.

        @param baseList: List of basepair chunks that make up the duplex.
        @type  baseList: list
        """
        # This implem depends on the specifics of how the end-representations
        # are terminated. If that's changed, it might stop working or it might
        # start giving wrong results. In the current representation,
        # baseList[0] (a chunk) has two bonds whose directions we must set,
        # which will determine the directions of their strands:
        #   Ss5 -> Sh5, and Ss5 <- Pe5.
        # Just find those bonds and set the strand directions (until such time
        # as they can be present to start with in the end1 mmp file).
        # (If we were instead passed all the atoms, we could be correct if we
        # just did this to the first Pe5 and Sh5 we saw, or to both of each if
        # setting the same direction twice is allowed.)
        atoms = baseList[0].atoms.values()
        Pe_list = filter( lambda atom: atom.element.symbol in ('Pe5'), atoms)
        Sh_list = filter( lambda atom: atom.element.symbol in ('Sh5'), atoms)

        if len(Pe_list) == len(Sh_list) == 1:
            for atom in Pe_list:
                assert len(atom.bonds) == 1
                atom.bonds[0].set_bond_direction_from(atom, 1, propogate = True)
            for atom in Sh_list:
                assert len(atom.bonds) == 1
                atom.bonds[0].set_bond_direction_from(atom, -1, propogate = True)
        else:
            #bruce 070604 mitigate bug in above code when number of bases == 1
            # by not raising an exception when it fails.
            msg = "Warning: strand not terminated, bond direction not set \
            (too short)"
            env.history.message( orangemsg( msg))

            # Note: It turns out this bug is caused by a bug in the rest of the
            # generator (which I didn't try to diagnose) -- for number of
            # bases == 1 it doesn't terminate the strands, so the above code
            # can't find the termination atoms (which is how it figures out
            # what to do without depending on intimate knowledge of the base
            # mmp file contents).

            # print "baseList = %r, its len = %r, atoms in [0] = %r" % \
            #       (baseList, len(baseList), atoms)
            ## baseList = [<molecule 'unknown' (11 atoms) at 0xb3d6f58>],
            ## its len = 1, atoms in [0] = [Ax1, X2, X3, Ss4, Pl5, X6, X7, Ss8, Pl9, X10, X11]

            # It would be a mistake to fix this here (by giving it that
            # intimate knowledge) -- instead we need to find and fix the bug
            # in the rest of generator when number of bases == 1.
        return
    pass

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

    def _postProcess(self, baseList): # bruce 070414
        """
        Set bond direction on the backbone bonds.

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

        @note: baseList must contain at least two base-pair chunks.
        """
        # This implem depends on the specifics of how the end-representations
        # are terminated. If that's changed, it might stop working or it might
        # start giving wrong results. In the current representation,
        # baseList[0] (a chunk) is the starting end of the duplex. It has two
        # bonds whose directions we should set, which will determine the
        # directions of their strands: Se3 -> Ss3, and Sh3 <- Ss3.
        # Just find those bonds and set the strand directions.

        assert len(baseList) >= 2
        basepair1_atoms = baseList[0].atoms.values() # StartBasePair.MMP
        basepair2_atoms = baseList[1].atoms.values() # MiddleBasePair or EndBasePair.MMP
        Se3_list = filter( lambda atom: atom.element.symbol in ('Se3'), basepair1_atoms)
        Sh3_list = filter( lambda atom: atom.element.symbol in ('Sh3'), basepair1_atoms)
        # To set the direction of the Se3 -> Ss3 bond, we need the Ss3 atom
        # from the second base-pair (i.e. baseList[1]) that is connected to
        # the Se3 atom in the first base-pair.
        # Ss3_list will have only one Ss3 atom if the duplex is only two
        # base-pairs long. Otherwise, Ss3_list will have two Ss3 atoms.
        Ss3_list = filter( lambda atom: atom.element.symbol in ('Ss3'), basepair2_atoms)
        if len(Se3_list) == len(Sh3_list) == 1:
            for atom in Se3_list:
                assert len(atom.bonds) == 2
                # This is fragile since it is dependent on the bond order in the
                # PAM3 StartBasePair.MMP file. If the order gets changed in that
                # file, this will fail.
                # A more robust approach would be to check for the Se3 -> Ss3
                # bond, which is the one we want here. Mark 2007-09-27.
                #atom.bonds[1].set_bond_direction_from(atom, 1, propogate = True)
                # This implements the more robust stategy mentioned
                # above. I'd like Bruce to review it and confirm it is good.
                # Mark 2007-09-27.
                #atom.bonds[1].set_bond_direction_from(Ss3_list[0], -1, propogate = True)
                try: # Try the first Ss3 atom in Ss3_list.
                    atom.bonds[1].set_bond_direction_from(Ss3_list[0], -1, propogate = True)
                except: # That wasn't it, so it must be the second Ss3 atom.
                    atom.bonds[1].set_bond_direction_from(Ss3_list[1], -1, propogate = True)
            for atom in Sh3_list:
                assert len(atom.bonds) == 1
                atom.bonds[0].set_bond_direction_from(atom, -1, propogate = True)
        else:
            #bruce 070604 mitigate bug in above code when number of bases == 1
            # by not raising an exception when it fails.
            msg = "Warning: strand not terminated, bond direction not set \
            (too short)"
            env.history.message( orangemsg( msg))
        return

class Z_Dna(Dna):
    """
    Provides an atomistic model of the Z form of DNA.
    """

    form       =  "Z-DNA"
    baseRise   =  dnaDict['Z-DNA']['DuplexRise']
    handedness =  LEFT_HANDED

class Z_Dna_PAM5(Z_Dna):
    """
    Provides an atomistic model of the Z form of DNA.

    @attention: This class is not implemented yet.
    """
    model      =  "Atomistic"

    def _strandAinfo(self, baseLetter, index):
        """
        Returns parameters needed to add a base to strand A.

        @param baseLetter: The base letter.
        @type  baseLetter: str

        @param index: Index in base sequence.
        @type  index: int
        """

        thetaOffset  =  0.0
        basename  =  basesDict[baseLetter]['Name']

        if (index & 1) != 0:
            # Index is odd.
            basename  +=  "-outer"
            zoffset    =  2.045
        else:
            # Index is even.
            basename  +=  '-inner'
            zoffset    =  0.0

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

    def _strandBinfo(self, baseLetter, index):
        """
        Returns parameters needed to add a base to strand B.

        @param baseLetter: The base letter.
        @type  baseLetter: str

        @param index: Index in base sequence. This is unused.
        @type  index: int
        """

        thetaOffset  =  0.5 * pi
        basename     =  basesDict[baseLetter]['Name']

        if (index & 1) != 0:
            basename  +=  '-inner'
            zoffset    =  -0.055
        else:
            basename  +=  "-outer"
            zoffset    =  -2.1

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

class Z_Dna_PAM5(Z_Dna):
    """
    Provides a PAM5 reduced model of the Z form of DNA.

    @attention: This class is not implemented yet.
    """
    pass