InteractionEnergyVrna.h

Go to the documentation of this file.

 
#ifndef INTARNA_INTERACTIONENERGYVIENNA_H_
#define INTARNA_INTERACTIONENERGYVIENNA_H_
 
#include "IntaRNA/InteractionEnergy.h"
#include "IntaRNA/VrnaHandler.h"
 
extern "C" {
    #include <ViennaRNA/utils.h>
    #include <ViennaRNA/fold_vars.h>
    #include <ViennaRNA/model.h>
    #include <ViennaRNA/params.h>
    #include <ViennaRNA/loop_energies.h>
}
#ifndef VIENNA_RNA_PAIR_MAT_H
#define VIENNA_RNA_PAIR_MAT_H
extern "C" {
    #include <ViennaRNA/pair_mat.h>
}
#endif
 
#include "IntaRNA/Matrix.h"
 
#define Evrna_2_E( e ) ( static_cast<E_type>(e) )
 
namespace IntaRNA {
 
// http://www.tbi.univie.ac.at/RNA/ViennaRNA/doc/RNAlib-2.3.0.pdf
 
class InteractionEnergyVrna: public InteractionEnergy {
 
public:
 
 
 
    InteractionEnergyVrna( const Accessibility & accS1
                    , const ReverseAccessibility & accS2
                    , VrnaHandler &vrnaHandler
                    , const size_t maxInternalLoopSize1 = 16
                    , const size_t maxInternalLoopSize2 = 16
                    , const bool initES = false
                    , const E_type energyAdd = Ekcal_2_E(0.0)
                    , const bool energyWithDangles = true
                    , const bool internalLoopGU = true
                );
 
    virtual ~InteractionEnergyVrna();
 
 
    virtual
    E_type
    getES1( const size_t i1, const size_t j1 ) const;
 
    virtual
    E_type
    getES2( const size_t i2, const size_t j2 ) const;
 
    virtual
    E_type
    getE_multiUnpaired( const size_t numUnpaired ) const;
 
    virtual
    E_type
    getE_multiHelix( const size_t j1, const size_t j2 ) const;
 
    virtual
    E_type
    getE_multiClosing() const;
 
    virtual
    E_type
    getE_init() const;
 
    virtual
    E_type
    getE_interLeft( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const;
 
 
    virtual
    E_type
    getE_danglingLeft( const size_t i1, const size_t i2 ) const;
 
 
    virtual
    E_type
    getE_danglingRight( const size_t j1, const size_t j2 ) const;
 
 
    virtual
    E_type
    getE_endLeft( const size_t i1, const size_t i2 ) const;
 
    virtual
    E_type
    getE_endRight( const size_t j1, const size_t j2 ) const;
 
    virtual
    E_type
    getEall1() const;
 
    virtual
    E_type
    getEall2() const;
 
    virtual
    Z_type
    getRT() const;
 
protected:
 
 
    vrna_md_t foldModel;
 
    vrna_param_t * foldParams;
 
    Z_type RT;
 
    const int bpCG;
 
    const int bpGC;
 
    typedef UpperTriangularMatrix<E_type> EsMatrix;
 
    EsMatrix * esValues1;
 
    EsMatrix * esValues2;
 
    mutable E_type Eall1;
 
    mutable E_type Eall2;
 
    bool
    isGC( const size_t i1, const size_t i2 ) const;
 
    void
    computeES( const Accessibility & acc, EsMatrix & esToFill );
 
    E_type
    computeIntraEall( const Accessibility & acc ) const;
};
 
 
 
inline
E_type
InteractionEnergyVrna::
getE_init() const
{
    // init term is sequence independent
    return Evrna_2_E(foldParams->DuplexInit);
}
 
 
inline
E_type
InteractionEnergyVrna::
getE_endLeft( const size_t i1, const size_t i2 ) const
{
    // VRNA non-GC penalty
    return Evrna_2_E( isGC(i1,i2) ? 0.0 : foldParams->TerminalAU );
}
 
 
inline
E_type
InteractionEnergyVrna::
getE_endRight( const size_t j1, const size_t j2 ) const
{
    // VRNA non-GC penalty
    return Evrna_2_E( isGC(j1,j2) ? 0.0 : foldParams->TerminalAU );
}
 
 
inline
Z_type
InteractionEnergyVrna::
getRT() const
{
    return RT;
}
 
 
inline
E_type
InteractionEnergyVrna::
getE_interLeft( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const
{
    // if valid internal loop
    if ( isValidInternalLoop(i1,j1,i2,j2) ) {
        assert( i1!=j1 && i2!=j2 );
        // Vienna RNA : compute internal loop / stacking energy for base pair [i1,i2]
        return Evrna_2_E(E_IntLoop( (int)j1-i1-1    // unpaired region 1
                            , (int)j2-i2-1  // unpaired region 2
                            , BP_pair[accS1.getSequence().asCodes().at(i1)][accS2.getSequence().asCodes().at(i2)]   // type BP (i1,i2)
                            , BP_pair[accS2.getSequence().asCodes().at(j2)][accS1.getSequence().asCodes().at(j1)]   // type BP (j2,j1)
                            , accS1.getSequence().asCodes().at(i1+1)
                            , accS2.getSequence().asCodes().at(i2+1)
                            , accS1.getSequence().asCodes().at(j1-1)
                            , accS2.getSequence().asCodes().at(j2-1)
                            , foldParams))
                ;
    } else {
        return E_INF;
    }
}
 
 
inline
E_type
InteractionEnergyVrna::
getE_danglingLeft( const size_t i1, const size_t i2 ) const
{
    // Vienna RNA : dangling end contribution
    return Evrna_2_E(vrna_E_ext_stem( BP_pair[accS1.getSequence().asCodes().at(i1)][accS2.getSequence().asCodes().at(i2)]
                              , ( i1==0 ? -1 : accS1.getSequence().asCodes().at(i1-1) )
                              , ( i2==0 ? -1 : accS2.getSequence().asCodes().at(i2-1) )
                              , foldParams
                              ))
            // substract closing penalty
            - getE_endLeft(i1,i2);
}
 
 
inline
E_type
InteractionEnergyVrna::
getE_danglingRight( const size_t j1, const size_t j2 ) const
{
    // Vienna RNA : dangling end contribution (reverse base pair to be sequence end conform)
    return Evrna_2_E(vrna_E_ext_stem( BP_pair[accS2.getSequence().asCodes().at(j2)][accS1.getSequence().asCodes().at(j1)]
                              , ( j2+1>=accS2.getSequence().size() ? -1 : accS2.getSequence().asCodes().at(j2+1) )
                              , ( j1+1>=accS1.getSequence().size() ? -1 : accS1.getSequence().asCodes().at(j1+1) )
                              , foldParams
                              ))
            // substract closing penalty
            - getE_endRight(j1,j2);
}
 
 
inline
E_type
InteractionEnergyVrna::
getES1( const size_t i1, const size_t j1 ) const
{
#if INTARNA_IN_DEBUG_MODE
    // sanity check
    if (i1>j1) throw std::runtime_error("InteractionEnergy::getES1(i1="+toString(i1)+" > j1="+toString(j1));
    if (j1>=size1()) throw std::runtime_error("InteractionEnergy::getES1() : j1="+toString(j1)+" >= size1()="+toString(size1()));
    if (esValues1 == NULL) throw std::runtime_error("InteractionEnergy::getES1() : ES values not initialized");
#endif
 
    // return computed value
    return (*esValues1)(i1,j1);
}
 
 
inline
E_type
InteractionEnergyVrna::
getES2( const size_t i2, const size_t j2 ) const
{
#if INTARNA_IN_DEBUG_MODE
    // sanity check
    if (i2>j2) throw std::runtime_error("InteractionEnergy::getES2(i2="+toString(i2)+" > j2="+toString(j2));
    if (j2>=size2()) throw std::runtime_error("InteractionEnergy::getES2() : j2="+toString(j2)+" >= size2()="+toString(size2()));
    if (esValues2 == NULL) throw std::runtime_error("InteractionEnergy::getES2() : ES values not initialized");
#endif
 
    // return computed value
    return (*esValues2)(i2,j2);
}
 
 
inline
bool
InteractionEnergyVrna::
isGC( const size_t i1, const size_t i2 ) const
{
    const int bpType = BP_pair[accS1.getSequence().asCodes().at(i1)][accS2.getSequence().asCodes().at(i2)];
    return (bpType==bpCG || bpType==bpGC);
}
 
 
inline
E_type
InteractionEnergyVrna::
getE_multiUnpaired( const size_t numUnpaired ) const
{
    return E_type(numUnpaired) * (Evrna_2_E(foldParams->MLbase));
}
 
 
inline
E_type
InteractionEnergyVrna::
getE_multiHelix( const size_t j1, const size_t j2 ) const
{
    return (Evrna_2_E(foldParams->MLintern[
                                         BP_pair[accS2.getSequence().asCodes().at(j2)]
                                                 [accS1.getSequence().asCodes().at(j1)]
                                        ]));
}
 
 
inline
E_type
InteractionEnergyVrna::
getE_multiClosing() const
{
    return (Evrna_2_E(foldParams->MLclosing));
}
 
 
inline
E_type
InteractionEnergyVrna::
getEall1() const
{
    // compute Z if needed
    if (E_isINF(Eall1)) {
        Eall1 = computeIntraEall( accS1 );
    }
    return Eall1;
}
 
 
inline
E_type
InteractionEnergyVrna::
getEall2() const
{
    // compute Z if needed
    if (E_isINF(Eall2)) {
        Eall2 = computeIntraEall( accS2.getAccessibilityOrigin() );
    }
    return Eall2;
}
 
 
 
} // namespace
 
#endif /* INTERACTIONENERGYVIENNA_H_ */