RnaSequence.h

Go to the documentation of this file.

 
#ifndef INTARNA_RNASEQUENCE_H_
#define INTARNA_RNASEQUENCE_H_
 
#include <locale>
#include <string>
#include <vector>
 
#include "IntaRNA/general.h"
 
 
#ifndef VIENNA_RNA_PAIR_MAT_H
#define VIENNA_RNA_PAIR_MAT_H
extern "C" {
    #include <ViennaRNA/pair_mat.h>
}
#endif
 
namespace IntaRNA {
 
class RnaSequence {
 
public:
 
    typedef std::string String_type;
 
    typedef short Code_type;
 
    typedef std::vector<Code_type> CodeSeq_type;
 
 
    const static std::string SequenceAlphabetIUPAC;
 
    const static std::string SequenceAlphabet;
 
    static const size_t lastPos;
 
public:
 
    RnaSequence(const std::string& id
            , const std::string & seqString
            , const long idxPos0 = 1
            , const size_t seqNumber = 1 );
 
    virtual ~RnaSequence();
 
 
 
    const std::string&
    getId() const;
 
    const size_t
    getSeqNumber() const;
 
    size_t
    size() const;
 
    long
    getInOutIndex( const size_t i ) const;
 
    size_t
    getIndex( const long i ) const;
 
    size_t
    getReversedIndex( const size_t i ) const;
 
    const String_type&
    asString() const;
 
    const CodeSeq_type&
    asCodes() const;
 
    bool
    isAmbiguous() const;
 
    bool
    isAmbiguous( const size_t i ) const;
 
    const bool operator == ( const RnaSequence &rna2 ) const;
 
    friend std::ostream& operator<<(std::ostream& out, const RnaSequence& rna);
 
 
 
public:
 
 
    static
    String_type
    getUpperCase( const std::string & seqString );
 
 
    static
    CodeSeq_type
    getCodeForString( const String_type& seqString );
 
    static
    Code_type
    getCodeForChar( const char nucleotide );
 
 
    static
    bool
    isValidSequence( const std::string& sequence );
 
 
    static
    bool
    isValidSequenceIUPAC( const std::string& sequence );
 
 
    static
    bool
    areComplementary( const RnaSequence & s1, const RnaSequence & s2,
                    const size_t p1, const size_t p2 );
 
    static
    bool
    isGU( const RnaSequence & s1, const RnaSequence & s2,
                    const size_t p1, const size_t p2 );
 
protected:
 
 
 
    static std::locale codeLocale;
 
    static int bpGUcodes[];
 
    std::string id;
 
    String_type seqString;
 
    CodeSeq_type seqCode;
 
    bool ambiguous;
 
    long idxPos0;
 
    size_t seqNumber;
 
};
 
 
 
 
inline
RnaSequence::RnaSequence(
        const std::string & id
        , const std::string & seqString
        , const long idxPos0
        , const size_t seqNumber )
 :
    id(id)
    , seqString(getUpperCase(seqString))
    , seqCode(getCodeForString(this->seqString))
    , ambiguous(this->seqString.contains('N'))
    , idxPos0(idxPos0)
    , seqNumber(seqNumber)
{
#if INTARNA_IN_DEBUG_MODE
    if (id.size() == 0) {
        throw std::runtime_error("RnaSequence::RnaSequence : id empty");
    }
    if (seqString.size() == 0) {
        throw std::runtime_error("RnaSequence::RnaSequence : seqString empty");
    }
#endif
 
}
 
 
inline
RnaSequence::~RnaSequence()
{
}
 
 
inline
const std::string&
RnaSequence::
getId() const
{
    return id;
}
 
 
inline
const size_t
RnaSequence::
getSeqNumber() const
{
    return seqNumber;
}
 
 
inline
size_t
RnaSequence::
size() const
{
    return seqString.size();
}
 
 
inline
long
RnaSequence::
getInOutIndex( const size_t i ) const
{
#if INTARNA_IN_DEBUG_MODE
    if (i >= size()) {
        throw std::runtime_error("RnaSequence::getInOutIndex : index "+toString(i)+" >= length "+toString(size()));
    }
#endif
    // get in/out index
    long p = idxPos0 + (long)i;
    // check for -+1 index transition
    if (idxPos0<0 && p>=0) {
        p++;
    }
    // final in/out index
    return p;
}
 
 
inline
size_t
RnaSequence::
getIndex( const long i ) const
{
    // check lower bounds
    if (i < idxPos0) {
        throw std::runtime_error("RnaSequence::getIndex : index "+toString(i)+" < idxPos0 "+toString(idxPos0));
    }
    // shift to internal indexing
    size_t p = (size_t)(i - idxPos0);
    // check for -+1 index transition
    if (idxPos0<0 && i>=0) {
        assert(i>0);
        p--;
    }
    // check upper bound
    if (p >= size()) {
        throw std::runtime_error("RnaSequence::getIndex : index "+toString(i)+" relates to "+toString(p)+" >= length "+toString(size()));
    }
    // return internal index
    return p;
}
 
 
inline
size_t
RnaSequence::
getReversedIndex( const size_t i ) const
{
#if INTARNA_IN_DEBUG_MODE
    if (i >= size()) {
        throw std::runtime_error("RnaSequence::getReversedIndex : index "+toString(i)+" >= seq.length "+toString(size()));
    }
#endif
 
    return this->size() -i -1;
}
 
 
inline
const
RnaSequence::
String_type&
RnaSequence::
asString() const
{
    return seqString;
}
 
 
inline
const
RnaSequence::
CodeSeq_type&
RnaSequence::
asCodes() const
{
    return seqCode;
}
 
 
inline
bool
RnaSequence::
isAmbiguous() const
{
    return ambiguous;
}
 
 
inline
bool
RnaSequence::
isAmbiguous( const size_t i ) const
{
    // check for ambiguous nucleotide encoding
    return this->seqString.at(i) == 'N';
}
 
 
inline
RnaSequence::
String_type
RnaSequence::
getUpperCase( const std::string & seqString )
{
#if INTARNA_IN_DEBUG_MODE
    if (!isValidSequenceIUPAC(seqString)) throw std::runtime_error("RnaSequence::getUpperCase() : the given sequence contains non-IUPAC codes : '"+seqString+"'");
#endif
 
    // create container to fill
    String_type seqRet(seqString.size(),'_');
 
    for (size_t i=0; i<seqString.size(); ++i)
    {
        // get upper case characters
        seqRet[i] = std::toupper(seqString.at(i),codeLocale);
        if (seqRet[i]=='T') {
            seqRet[i] = 'U';
        }
        // overwrite non-ACGU characters with N = ambiguous (no distinction needed)
        if (!SequenceAlphabet.contains(seqRet[i])) {
            seqRet[i] = 'N';
        }
    }
 
    return seqRet;
}
 
 
inline
RnaSequence::
CodeSeq_type
RnaSequence::
getCodeForString( const String_type& seqString )
{
    // create container to fill and init with 'X'
    CodeSeq_type seqCode(seqString.size());
 
    for (size_t i=0; i<seqString.size(); ++i)
    {
        seqCode[i] = getCodeForChar( seqString.at(i) );
    }
 
    return seqCode;
}
 
 
inline
bool
RnaSequence::
isValidSequence( const std::string & sequence )
{
    // check whether or not the string contains unsupported characters
    return sequence.find_first_not_of(SequenceAlphabet) == std::string::npos;
}
 
 
inline
bool
RnaSequence::
isValidSequenceIUPAC( const std::string & sequence )
{
    // check whether or not the string contains unsupported characters
    return sequence.find_first_not_of(SequenceAlphabetIUPAC) == std::string::npos;
}
 
 
inline
std::ostream&
operator<<(std::ostream& out, const RnaSequence& rna)
{
    // add ID(SEQUENCE) to stream
    out <<rna.id <<'(' <<rna.asString() <<')';
    return out;
}
 
 
inline
RnaSequence::
Code_type
RnaSequence::
getCodeForChar( const char nucleotide )
{
#if INTARNA_IN_DEBUG_MODE
    // check if nucleotide character is NOT supported
    if (!SequenceAlphabet.contains(nucleotide))
        throw std::runtime_error("RnaSequence::getCodeForChar() : unsupported nucleotide character '"+toString(nucleotide)+"' in sequence");
#endif
 
    // otherwise get encoding:
    // use nucleotide character encoding from Vienna RNA package
    return (Code_type)encode_char(nucleotide);
}
 
 
inline
bool
RnaSequence::
areComplementary( const RnaSequence & s1, const RnaSequence & s2,
                    const size_t p1, const size_t p2 )
{
#if INTARNA_IN_DEBUG_MODE
    // check if valid positions
    if (p1>=s1.size() || p2>=s2.size())
        throw std::runtime_error("RnaSequence::areComplementary : index positions p1/p2 ("
                + toString(p1)+"/"+toString(p2)
                + ") are out of bounds s1/s2 ("
                + toString(s1.size())+"/"+toString(s2.size())
                +")"
                );
#endif
 
    // check via VRNA util
    return BP_pair[s1.seqCode.at(p1)][s2.seqCode.at(p2)] > 0;
}
 
 
inline
bool
RnaSequence::
isGU( const RnaSequence & s1, const RnaSequence & s2,
                    const size_t p1, const size_t p2 )
{
#if INTARNA_IN_DEBUG_MODE
    // check if valid positions
    if (p1>=s1.size() || p2>=s2.size())
        throw std::runtime_error("RnaSequence::areComplementary : index positions p1/p2 ("
                + toString(p1)+"/"+toString(p2)
                + ") are out of bounds s1/s2 ("
                + toString(s1.size())+"/"+toString(s2.size())
                +")"
                );
#endif
 
    // get bp code
    int bpCode = BP_pair[s1.seqCode.at(p1)][s2.seqCode.at(p2)];
 
    // check if GU pair
    return bpCode == bpGUcodes[0] || bpCode == bpGUcodes[1] ;
}
 
 
inline
const bool
RnaSequence::
operator == ( const RnaSequence &rna2 ) const
{
    return  // check if same object (pointer)
            this == &rna2
            || (
                // ensure same lengths
                this->size() == rna2.size()
                // ids are identical (most likely shorter, i.e. faster check)
                && this->id == rna2.id
                // sequences identical
                && this->seqString == rna2.seqString
            );
}
 
 
} // namespace
 
 
#endif /* RNASEQUENCE_H_ */