InteractionEnergy.h

Go to the documentation of this file.

 
#ifndef INTARNA_INTERACTIONENERGY_H_
#define INTARNA_INTERACTIONENERGY_H_
 
 
#include "IntaRNA/general.h"
#include "IntaRNA/Interaction.h"
#include "IntaRNA/Accessibility.h"
#include "IntaRNA/ReverseAccessibility.h"
 
namespace IntaRNA {
 
class InteractionEnergy {
 
public:
 
    enum ES_multi_mode {
        ES_multi_1only,
        ES_multi_2only,
        ES_multi_both,
    };
 
    struct EnergyContributions {
    public:
        E_type loops;
        E_type init;
        E_type ED1;
        E_type ED2;
        E_type dangleLeft;
        E_type dangleRight;
        E_type endLeft;
        E_type endRight;
        E_type energyAdd;
    };
 
 
    InteractionEnergy( const Accessibility & accS1
            , const ReverseAccessibility & accS2
            , const size_t maxInternalLoopSize1
            , const size_t maxInternalLoopSize2
            , const E_type energyAdd
            , const bool energyWithDangle
            , const bool internalLoopGU
            );
 
    virtual ~InteractionEnergy();
 
 
    virtual
    E_type
    getE( const size_t i1, const size_t j1
            , const size_t i2, const size_t j2
            , const E_type hybridE ) const;
 
    virtual
    E_type
    getE( const Z_type Z ) const;
 
    virtual
    EnergyContributions
    getE_contributions( const Interaction & interaction ) const;
 
    virtual
    bool
    areComplementary( const size_t i1, const size_t i2 ) const;
 
    virtual
    bool
    isGU( const size_t i1, const size_t i2 ) const;
 
    virtual
    size_t
    size1() const;
 
    virtual
    size_t
    size2() const;
 
    virtual
    E_type
    getED1( const size_t i1, const size_t j1 ) const;
 
    virtual
    E_type
    getED2( const size_t i2, const size_t j2 ) const;
 
    virtual
    bool
    isAccessible1( const size_t i ) const;
 
    virtual
    bool
    isAccessible2( const size_t i ) const;
 
 
    virtual
    E_type
    getE_multi(  const size_t i1, const size_t j1
                , const size_t i2, const size_t j2
                , const ES_multi_mode ES_mode ) const;
 
    virtual
    E_type
    getES1( const size_t i1, const size_t j1 ) const = 0;
 
    virtual
    E_type
    getES2( const size_t i2, const size_t j2 ) const = 0;
 
    virtual
    E_type
    getE_multiUnpaired( const size_t numUnpaired ) const = 0;
 
    virtual
    E_type
    getE_multiHelix( const size_t j1, const size_t j2 ) const = 0;
 
    virtual
    E_type
    getE_multiClosing() const = 0;
 
    virtual
    E_type
    getE_init( ) const = 0;
 
    virtual
    E_type
    getE_interLeft( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const = 0;
 
 
    virtual
    E_type
    getE_danglingLeft( const size_t i1, const size_t i2 ) const = 0;
 
 
    virtual
    E_type
    getE_danglingRight( const size_t j1, const size_t j2 ) const = 0;
 
    virtual
    E_type
    getE_endLeft( const size_t i1, const size_t i2 ) const = 0;
 
    virtual
    E_type
    getE_endRight( const size_t j1, const size_t j2 ) const = 0;
 
    virtual
    Z_type
    getPr_danglingLeft( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const;
 
    virtual
    Z_type
    getPr_danglingRight( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const;
 
    virtual
    const Accessibility &
    getAccessibility1() const;
 
    virtual
    const ReverseAccessibility &
    getAccessibility2() const;
 
    const size_t getMaxInternalLoopSize1() const {
        return maxInternalLoopSize1;
    }
 
    const size_t getMaxInternalLoopSize2() const {
        return maxInternalLoopSize2;
    }
 
    virtual
    Z_type
    getRT() const = 0;
 
    virtual
    Z_type
    getBoltzmannWeight( const E_type energy ) const ;
 
    virtual
    Z_type
    getBoltzmannWeight( const Z_type energy ) const ;
 
 
    virtual
    Interaction::BasePair
    getBasePair( const size_t i1, const size_t i2 ) const;
 
 
    virtual
    size_t
    getIndex1( const Interaction::BasePair & bp ) const;
 
 
    virtual
    size_t
    getIndex2( const Interaction::BasePair & bp ) const;
 
    virtual
    E_type
    getEnergyAdd() const;
 
    virtual
    bool
    isValidInternalLoop( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const;
 
    bool
    isInternalLoopGUallowed() const;
 
    virtual
    E_type
    getEall1() const = 0;
 
    virtual
    E_type
    getEall2() const = 0;
 
protected:
 
    const Accessibility & accS1;
 
    const ReverseAccessibility & accS2;
 
    const size_t maxInternalLoopSize1;
 
    const size_t maxInternalLoopSize2;
 
    const E_type energyAdd;
 
    const bool energyWithDangles;
 
    const bool internalLoopGU;
 
    static
    bool
    isAllowedLoopRegion( const RnaSequence& seq, const size_t i, const size_t j, const size_t maxInternalLoopSize );
 
};
 
 
 
 
 
 
inline
InteractionEnergy::InteractionEnergy( const Accessibility & accS1
                , const ReverseAccessibility & accS2
                , const size_t maxInternalLoopSize1
                , const size_t maxInternalLoopSize2
                , const E_type energyAdd
                , const bool energyWithDangles
                , const bool internalLoopGU
        )
  :
    accS1(accS1)
    , accS2(accS2)
    , maxInternalLoopSize1(maxInternalLoopSize1)
    , maxInternalLoopSize2(maxInternalLoopSize2)
    , energyAdd(energyAdd)
    , energyWithDangles(energyWithDangles)
    , internalLoopGU(internalLoopGU)
{
}
 
 
inline
InteractionEnergy::~InteractionEnergy()
{
}
 
 
inline
bool
InteractionEnergy::
isAllowedLoopRegion( const RnaSequence& seq, const size_t i, const size_t j, const size_t maxInternalLoopSize )
{
    // ensure index and loop size validity
    return     i < seq.size()
            && j < seq.size()
            && seq.asString().at(i) != 'N'
            && seq.asString().at(j) != 'N'
            && i <= j
            && (j-i) <= (1+maxInternalLoopSize);
 
}
 
 
inline
bool
InteractionEnergy::
areComplementary( const size_t i1, const size_t i2 ) const
{
    return RnaSequence::areComplementary( accS1.getSequence(), accS2.getSequence(), i1, i2)
        && isAccessible1(i1)
        && isAccessible2(i2);
}
 
 
inline
bool
InteractionEnergy::
isGU( const size_t i1, const size_t i2 ) const
{
    return RnaSequence::isGU( accS1.getSequence(), accS2.getSequence(), i1, i2);
}
 
 
inline
size_t
InteractionEnergy::
size1() const
{
    return getAccessibility1().getSequence().size();
}
 
 
inline
size_t
InteractionEnergy::
size2() const
{
    return getAccessibility2().getSequence().size();
}
 
 
inline
bool
InteractionEnergy::
isValidInternalLoop( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const
{
    return
           (j1-i1>0 && j2-i2>0)
        && areComplementary( i1, i2)
        && areComplementary( j1, j2)
        && InteractionEnergy::isAllowedLoopRegion(accS1.getSequence(), i1, j1, maxInternalLoopSize1)
        && InteractionEnergy::isAllowedLoopRegion(accS2.getSequence(), i2, j2, maxInternalLoopSize2)
        && ( internalLoopGU || (i1+1==j1 && i2+1==j2) || (!isGU(i1,i2) && !isGU(j1,j2)) ) // GU-allowed or stacking or no GU
        ;
}
 
 
inline
bool
InteractionEnergy::
isInternalLoopGUallowed() const
{
    return internalLoopGU;
}
 
 
inline
const Accessibility &
InteractionEnergy::
getAccessibility1() const
{
    return accS1;
}
 
 
inline
const ReverseAccessibility &
InteractionEnergy::
getAccessibility2() const
{
    return accS2;
}
 
 
inline
E_type
InteractionEnergy::
getEnergyAdd() const
{
    return energyAdd;
}
 
 
inline
Z_type
InteractionEnergy::
getBoltzmannWeight( const E_type e ) const
{
    // TODO can be optimized when using exp-energies from VRNA
    return Z_exp( - E_2_Z(e) / getRT() );
}
 
 
inline
Z_type
InteractionEnergy::
getBoltzmannWeight( const Z_type e ) const
{
    // TODO can be optimized when using exp-energies from VRNA
    return Z_exp( - e / getRT() );
}
 
 
inline
Interaction::BasePair
InteractionEnergy::
getBasePair( const size_t i1, const size_t i2 ) const
{
    return Interaction::BasePair( i1, getAccessibility2().getReversedIndex(i2) );
}
 
 
inline
size_t
InteractionEnergy::
getIndex1( const Interaction::BasePair & bp ) const
{
    return bp.first;
}
 
 
inline
size_t
InteractionEnergy::
getIndex2( const Interaction::BasePair & bp ) const
{
    return getAccessibility2().getReversedIndex( bp.second );
}
 
 
inline
E_type
InteractionEnergy::
getED1( const size_t i1, const size_t j1 ) const
{
    return getAccessibility1().getED( i1, j1 );
}
 
 
inline
E_type
InteractionEnergy::
getED2( const size_t i2, const size_t j2 ) const
{
    return getAccessibility2().getED( i2, j2 );
}
 
 
inline
bool
InteractionEnergy::
isAccessible1( const size_t i ) const
{
    return
            (!getAccessibility1().getSequence().isAmbiguous(i))
            && getAccessibility1().getAccConstraint().isAccessible(i)
            ;
}
 
 
inline
bool
InteractionEnergy::
isAccessible2( const size_t i ) const
{
    return
            (!getAccessibility2().getSequence().isAmbiguous(i))
            && getAccessibility2().getAccConstraint().isAccessible(i);
}
 
 
inline
Z_type
InteractionEnergy::
getPr_danglingLeft( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const
{
    // initial probabilities
    Z_type probDangle1 = 1.0, probDangle2 = 1.0;
 
    // if dangle1 possible
    if (i1>0)  {
        // Pr( i1-1 is unpaired | i1..j1 unpaired )
        probDangle1 =
            std::max( (Z_type)0.0
                    , std::min( (Z_type)1.0
                            , getBoltzmannWeight( getED1(i1-1,j1)-getED1(i1,j1) )
                            )
                    )
            ;
    }
    // if dangle2 possible
    if (i2>0)  {
        // Pr( i2-1 is unpaired | i2..j2 unpaired )
        probDangle2 =
            std::max( (Z_type)0.0
                    , std::min( (Z_type)1.0
                            , getBoltzmannWeight( getED2(i2-1,j2)-getED2(i2,j2) )
                            )
                    )
            ;
    }
 
    // get overall probability
    return probDangle1 * probDangle2;
}
 
 
inline
Z_type
InteractionEnergy::
getPr_danglingRight( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const
{
    // initial probabilities
    Z_type probDangle1 = 1.0, probDangle2 = 1.0;
 
    // if dangle1 possible
    if (j1+1<size1())  {
        // Pr( j1+1 is unpaired | i1..j1 unpaired )
        probDangle1 =
            std::max( (Z_type)0.0
                    , std::min( (Z_type)1.0
                            , getBoltzmannWeight( getED1(i1,j1+1)-getED1(i1,j1) )
                            )
                    )
            ;
    }
    // if dangle2 possible
    if (j2+1<size2())  {
        // Pr( j2+1 is unpaired | i2..j2 unpaired )
        probDangle2 =
            std::max( (Z_type)0.0
                    , std::min( (Z_type)1.0
                            , getBoltzmannWeight( getED2(i2,j2+1)-getED2(i2,j2) )
                            )
                    )
            ;
    }
 
    // get overall probability
    return probDangle1 * probDangle2;
}
 
 
inline
E_type
InteractionEnergy::
getE( const size_t i1, const size_t j1
        , const size_t i2, const size_t j2
        , const E_type hybridE ) const
{
    // check if hybridization energy and EDs are not infinite
    if ( E_isNotINF(hybridE)
            && (getED1( i1, j1 ) < Accessibility::ED_UPPER_BOUND)
            && (getED2( i2, j2 ) < Accessibility::ED_UPPER_BOUND))
    {
        // compute overall interaction energy
        return hybridE
                // accessibility penalty
                + getED1( i1, j1 )
                + getED2( i2, j2 )
                // dangling end penalty
                // weighted by the probability that ends are unpaired
                + (energyWithDangles ? Z_2_E(E_2_Z(getE_danglingLeft( i1, i2 ))*getPr_danglingLeft(i1,j1,i2,j2)) : E_type(0))
                + (energyWithDangles ? Z_2_E(E_2_Z(getE_danglingRight( j1, j2 ))*getPr_danglingRight(i1,j1,i2,j2)) : E_type(0))
                // helix closure penalty
                + getE_endLeft( i1, i2 )
                + getE_endRight( j1, j2 )
                + getEnergyAdd()
                ;
    } else {
        // hybridE is infinite, thus overall energy is infinity as well
        return E_INF;
    }
}
 
 
inline
E_type
InteractionEnergy::
getE( const Z_type Z ) const
{
    // convert partition function to ensemble energy
    // convert to E_type
    return Z_2_E( - getRT() * Z_log( Z ) );
}
 
 
inline
E_type
InteractionEnergy::
getE_multi(  const size_t i1, const size_t j1
            , const size_t i2, const size_t j2
            , const ES_multi_mode ES_mode ) const
{
#if INTARNA_IN_DEBUG_MODE
    if (i1 >= j1 ) throw std::runtime_error("InteractionEnergy::getE_multi() : i1>=j1 : "+toString(i1)+" "+toString(j1));
    if (i2 >= j2 ) throw std::runtime_error("InteractionEnergy::getE_multi() : i2>=j2 : "+toString(i2)+" "+toString(j2));
    if (j1 >= size1()) throw std::runtime_error("InteractionEnergy::getE_multi() : j1>=size1() : "+toString(j1)+" "+toString(size1()));
    if (j2 >= size2()) throw std::runtime_error("InteractionEnergy::getE_multi() : j2>=size2() : "+toString(j2)+" "+toString(size2()));
    if (! areComplementary(i1,i2)) throw std::runtime_error("InteractionEnergy::getE_multi() : not complementary : "+toString(i1)+" "+toString(i2));
    if (! areComplementary(j1,j2)) throw std::runtime_error("InteractionEnergy::getE_multi() : not complementary : "+toString(j1)+" "+toString(j2));
#endif
 
    return
            // intramolecular structure contributions
              (ES_mode != ES_multi_2only ? getES1(i1,j1) : 0)
            + (ES_mode != ES_multi_1only ? getES2(i2,j2) : 0)
            // dangling end treatments of helix ends
            + (energyWithDangles ? getE_danglingRight(i1,i2) : E_type(0))
            + (energyWithDangles ? getE_danglingLeft(j1,j2) : E_type(0))
            // helix closure penalties
            + getE_endRight(i1,i2)
            + getE_endLeft(j1,j2)
            // multiloop unpaired contributions
            + getE_multiUnpaired(
                    (ES_mode == ES_multi_2only ? j1-i1-1 : 0 )
                    + (ES_mode == ES_multi_1only ? j2-i2-1 : 0 )
                    )
            // multiloop helix contribution (right side interaction site)
            + getE_multiHelix( j1, j2 )
            // multiloop closure
            + getE_multiClosing()
            ;
}
 
 
 
 
} // namespace
 
 
#endif /* INTERACTIONENERGY_H_ */