main page namespaces classes files
ELL home

src/ell/rna/RNAFreeEnergy_TB.cc

Go to the documentation of this file.

#include <biu/assertbiu.hh>
#include "ell/rna/RNAFreeEnergy_TB.hh"

extern "C" {
    #include <ViennaRNA/fold_vars.h> // for dangles parameter
    #include <ViennaRNA/pair_mat.h>  // for encode_char function
//  #include <ViennaRNA/fold.h>
//  #include <ViennaRNA/utils.h>
}
extern "C"  float energy_of_struct(const char *string, const char *structure);
extern "C"  void update_fold_params ();
extern "C"  int encode_char(char c);
extern "C"  short* make_pair_table(const char *structure);
extern "C"  int loop_energy(short * ptable, short *s, short *s1, int i);

void dummyCall() {
    make_pair_matrix(); // needed to avoid warning due to <ViennaRNA/pair_mat.h>
}

#include <cstring>  //For the string functions
#include <limits.h>

namespace ell
{

    const biu::Alphabet RNAFreeEnergy_TB::COMMON_ALPHABET("ACGU", 1);
    const biu::AllowedBasePairs RNAFreeEnergy_TB::COMMON_BPAIRS(
                                &COMMON_ALPHABET, "AU,CG,GU");
    const double RNAFreeEnergy_TB::ENERGY_INF((double)UINT_MAX);
    
    bool RNAFreeEnergy_TB::calculateEnergyChanges = true;
    
      // whether or not the fold parameters of the Vienna package are 
      // already set or not
    bool ViennaFoldParamsSet = false;

    
// RNAFreeEnergy_TB

    RNAFreeEnergy_TB::RNAFreeEnergy_TB(const std::string& rnaSeqStr
                            , const std::string& rnaStructBracketDotStr
                            , const double energy_ )
     :  biu::RNAStructure_TB(rnaSeqStr, rnaStructBracketDotStr, &COMMON_BPAIRS)
        , energy(ENERGY_INF)
        , seqVienna(NULL)
        , pairTable(NULL)
    {
          // ViennaRNA package parameter for dangling ends
        setViennaFoldParams();
          // inform that molecule was initialised
        moleculeChanged();
          // set energy if given
        if (energy_ != ENERGY_INF) {
            setEnergy(energy_);
        }
    }
    
    RNAFreeEnergy_TB::RNAFreeEnergy_TB(biu::Sequence* rnaSeq
                                , const biu::Structure* const rnaStructBracketDot
                                , const bool seqIsShared
                                , const double energy_ )
     :  biu::RNAStructure_TB(rnaSeq, rnaStructBracketDot, &COMMON_BPAIRS, seqIsShared)
        , energy(ENERGY_INF)
        , seqVienna(NULL)
        , pairTable(NULL)
    {
          // ViennaRNA package parameter for dangling ends
        setViennaFoldParams();
          // inform that molecule was initialised
        moleculeChanged();
          // set energy if given
        if (energy_ != ENERGY_INF) {
            setEnergy(energy_);
        }
    }
    
    RNAFreeEnergy_TB::RNAFreeEnergy_TB(const RNAFreeEnergy_TB& rnaFreeEnerg)
     :  biu::RNAStructure_TB(rnaFreeEnerg)
        , energy(rnaFreeEnerg.energy)
        , seqVienna(NULL)
        , pairTable(NULL)
    {
    }

    RNAFreeEnergy_TB::~RNAFreeEnergy_TB()   
    {
        if (seqVienna != NULL) delete seqVienna;
        if (pairTable != NULL) delete pairTable;
    }

    void        
    RNAFreeEnergy_TB::applySingleMoveInPlace(const SingleMove& move) {
        
          // check if energy should be updated
        if (calculateEnergyChanges && energy != ENERGY_INF) {
            
            size_t around = move.first;
            while (around > 0) {
                around--;
                if ( validTreeStruc->at(around).pair != INVALID_INDEX ) {
                    if ( validTreeStruc->at(around).pair > move.first ) {
                        break;
                    } else {
                          // jump over base pair
                        around = validTreeStruc->at(around).pair;
                    }
                }
            }
            const bool nothingAround = move.first == 0 || (around == 0 && validTreeStruc->at(around).pair == INVALID_INDEX);
              // create temporary data structure
            if (seqVienna == NULL || pairTable == NULL) {
                updateSeqVienna();
            }
            
              // temp vars for energy contributions
            double aroundE = 0.0;
            double moveE = 0.0;
            double stackE = 0.0;
            
//std::cerr <<" E before = " <<energy 
//          <<" eosPT = "
//          <<energy_of_struct_pt(COMMON_ALPHABET.getString(*rnaSeq).c_str(),pairTable, seqVienna, seqVienna)
//          <<"  eos = "
//          <<energy_of_struct(COMMON_ALPHABET.getString(*rnaSeq).c_str(),bracketStructStr.c_str());
//std::cerr <<"\n";
            // remove bond between move.first and move.second
            if (validTreeStruc->at(move.first).pair == move.second) {
if (around == 0)
//std::cerr <<"\n###  DELETION\n";
                stackE = loop_energy(pairTable, seqVienna, seqVienna, (nothingAround?0:around+1)); // around
                

                

                moveE = loop_energy(pairTable, seqVienna, seqVienna, (move.first+1));  // move
                deleteBond( move.first, move.second );
                  // update pair table if necessary
                pairTable[move.first+1] = 0;
                pairTable[move.second+1] = 0;
                aroundE = loop_energy(pairTable, seqVienna, seqVienna, (nothingAround?0:around+1)); // around
                  // calculate new energy
                energy += double(aroundE - stackE - moveE) / 100.0;
            } 
            // add bond
            else {
//std::cerr <<"\n###  INSERTION\n";
                assertbiu(isValidInsertMove( move.first, move.second ),
                    "can't close bond with single move, because it is not valid");
                aroundE = loop_energy(pairTable, seqVienna, seqVienna, (nothingAround?0:around+1)); // around
                insertBond( move.first, move.second );
                  // update pair table if necessary
                pairTable[move.first+1] = (short)move.second+1;
                pairTable[move.second+1] = (short)move.first+1;
                moveE = loop_energy(pairTable, seqVienna, seqVienna, (move.first+1)); // move
                stackE = loop_energy(pairTable, seqVienna, seqVienna, (nothingAround?0:around+1)); // around
                  // calculate new energy
                energy += double(moveE + stackE - aroundE) / 100.0;
            }
//std::cerr <<" E after = " <<energy 
//          <<" eosPT = "
//          <<energy_of_struct_pt(COMMON_ALPHABET.getString(*rnaSeq).c_str(),pairTable, seqVienna, seqVienna)
//          <<"  eos = "
//          <<energy_of_struct(COMMON_ALPHABET.getString(*rnaSeq).c_str(),bracketStructStr.c_str());

            assertbiu( energy_of_struct( 
                            COMMON_ALPHABET.getString(*rnaSeq).c_str()
                            , bracketStructStr.c_str() ) - energy < 0.0001
                , "updated energy differs from direct calculation with Vienna RNA package");

        }
          // energy should be recalculated
        else {
            // remove bond between move.first and move.second
            if (validTreeStruc->at(move.first).pair == move.second) {
                deleteBond( move.first, move.second );
                  // update pair table if necessary
                if (pairTable != NULL) {
                    pairTable[move.first+1] = 0;
                    pairTable[move.second+1] = 0;
                }
            } 
            // add bond
            else {
                assertbiu(isValidInsertMove( move.first, move.second ),
                    "can't close bond with single move, because it is not valid");
                insertBond( move.first, move.second );
                  // update pair table if necessary
                if (pairTable != NULL) {
                    pairTable[move.first+1] = (short)move.second+1;
                    pairTable[move.second+1] = (short)move.first+1;
                }
            }
            
              // inform that structure has changed
            moleculeChanged();
        }
    }

    RNAFreeEnergy_TB&   
    RNAFreeEnergy_TB::operator= (const RNAFreeEnergy_TB& rnaStruct2) {
        
        if (this != &rnaStruct2) {
              // call assignment operator of super class
            biu::RNAStructure_TB::operator=(rnaStruct2);
            
              // apply assignments for this class
            energy = rnaStruct2.energy;
            if (seqVienna != NULL) delete seqVienna;
            seqVienna = NULL;
            if (pairTable != NULL) delete pairTable;
            pairTable = NULL;
        }
        return *this;
    }
    
    void 
    RNAFreeEnergy_TB
	::updateSeqVienna() 
    {
        const size_t length = bracketStructStr.size();
        if (seqVienna == NULL) {
            seqVienna = new short[length+2];
//          S1= (short *) space(sizeof(short)*(l+2));
            /* S1 exists only for the special X K and I bases and energy_set!=0 */
            seqVienna[0] = (short) length;

            const std::string seq(COMMON_ALPHABET.getString(*rnaSeq));
            for (size_t i=1; i<=length; i++) { /* make numerical encoding of sequence */
                seqVienna[i] = (short) encode_char(seq.at(i-1));
//              S1[i] = alias[seqVienna[i]];   /* for mismatches of nostandard bases */
            }
            /* for circular folding add first base at position n+1 and last base at
            position 0 in S1    */
            seqVienna[length+1] = seqVienna[1]; 
//          S1[l+1]=S1[1]; S1[0] = S1[l];
            

        }
        if (pairTable == NULL) {
              // create pair_table
//          pairTable = new short[length+2];
//          pairTable[0] = (short)length;
//          for (size_t i=1; i<=length; i++) {
//              if (validTreeStruc->at(i-1).pair != INVALID_INDEX) {
//                  pairTable[i] = (short)validTreeStruc->at(i-1).pair +1;
//              } else {
//                  pairTable[i] = 0;
//              }
//          }
            pairTable = make_pair_table(bracketStructStr.c_str());
        }

//std::cerr <<"\n#### cur seqVienna: l = " <<seqVienna[0] <<"\n";
//for (short i=0; i < seqVienna[0]; i++) {
//  std::cerr <<seqVienna[i] <<" ";
//}
//std::cerr <<"\n#### cur pairTable: l = " <<pairTable[0] <<"\n";
//for (short i=0; i < pairTable[0]; i++) {
//  std::cerr <<pairTable[i] <<" ";
//}
//std::cerr <<"\n";

    }

    
    void 
    RNAFreeEnergy_TB
	::setViennaFoldParams(void) {
          // set dangling end handling
        dangles = 2;
          // force the update in the vienna lib
        update_fold_params();
    }

    
} // namespace ell