geneticcode.hh

/*
 * geneticcode.hh
 *
 * License: Artistic License, see file LICENSE.TXT or 
 *          https://opensource.org/licenses/artistic-license-1.0
 */
 
#ifndef _GENETIC_CODE_HH
#define _GENETIC_CODE_HH
 
// project includes
#include "types.hh"
#include "projectio.hh"
 
// standard C/C++ includes
#include <string>
#include <cstring>
#include <iostream>
#include <fstream>
 
using namespace std;
 
 
#define NUM_AA 20
#define NUM_TRANSTABS 24
 
/*
 * === am I at a possible splice site or translation start site ===
 */
#define DECLARE_ON(NAME, PATTERN, COUNT)        \
    inline bool NAME(const char* dna) {         \
    return strncmp(dna, PATTERN, COUNT) == 0;   \
    }
// 'gt'
#define DSS_SEQUENCE "gt"
#define RDSS_SEQUENCE "ac"
DECLARE_ON(onDSS,     DSS_SEQUENCE, 2)
DECLARE_ON(onRDSS,    RDSS_SEQUENCE, 2)
 
// 'gc'
#define ALT_DSS_SEQUENCE "gc"
#define ALT_RDSS_SEQUENCE "gc"
DECLARE_ON(onAltDSS,  ALT_DSS_SEQUENCE, 2)
DECLARE_ON(onAltRDSS, ALT_RDSS_SEQUENCE, 2)
 
// 'gt' or 'gc'
inline bool onGenDSS(const char* dna) {
    return 
    onDSS(dna) || (Constant::dss_gc_allowed && onAltDSS(dna));
}
inline bool onGenRDSS(const char* dna) {
    return 
    onRDSS(dna) || (Constant::dss_gc_allowed && onAltRDSS(dna));
}
 
// 'ag'
#define ASS_SEQUENCE "ag"
#define RASS_SEQUENCE "ct"
DECLARE_ON(onASS,     ASS_SEQUENCE, 2)
DECLARE_ON(onRASS,    RASS_SEQUENCE, 2)
 
// 'atg'
#define STARTCODON "atg"
#define RCSTARTCODON "cat"
DECLARE_ON(onStart,   STARTCODON, 3)
DECLARE_ON(onRStart,  RCSTARTCODON, 3)
 
// stop codons
DECLARE_ON(onOchre,   OCHRECODON, 3)
DECLARE_ON(onAmber,   AMBERCODON, 3)
DECLARE_ON(onOpal,    OPALCODON, 3)
DECLARE_ON(onROchre,  RCOCHRECODON, 3)
DECLARE_ON(onRAmber,  RCAMBERCODON, 3)
DECLARE_ON(onROpal,   RCOPALCODON, 3)
 
/*
 * wcComplement
 *
 * Watson/Crick complement of character. Preserves case.
 * Doesn't do anything if character is non-acgt
 *
 */
inline char wcComplement(char c) {
    switch (c) {
        case 'a': return 't';
        case 'c': return 'g';
        case 'g': return 'c';
        case 't': return 'a';
        case 'A': return 'T';
        case 'C': return 'G';
        case 'G': return 'C';
        case 'T': return 'A';
        default:
        return c;
    }
}
 
 
/*
 * putReverseComplement
 *
 * copies the reverse complement of [dna, dna+len) to result
 * result may be a char* or a string iterator
 *
 * no validation is performed and no null character is appended
 */
template <class Iterator>
inline void putReverseComplement(Iterator result, const char* dna, int len) {
    const char* s = dna + len;
    Iterator t = result;
    while (--s >= dna) 
    *t++ = wcComplement(*s);
}
 
/*
 * reverseComplement
 * 
 * allocates a new character array and fills it with the reverse complement of 'dna'
 * eg agcngt -> acngct
 * doesn't change uppercase or lowercase property
 */
inline char* reverseComplement(const char* dna) {
    if (dna == NULL)
    return NULL;
    int len = strlen(dna);
    char* result = new char[len+1];
    putReverseComplement(result, dna, len);
    result[len] = '\0';
    return result;
}
 
inline void reverseComplementString(string &text) {
    int n = text.length();
    for (int i=0; i < n/2; i++) {
        char c;
        c = text[i];
        text[i] = wcComplement(text[n-i-1]);
        text[n-i-1] = wcComplement(c);
    }
    if (n%2) // n odd, must complement the middle character
    text[n/2] = wcComplement(text[n/2]);
}
 
 
inline void reverseString(string &text) {
    int i = 0;
    int n = text.length();
    while (i < (n/2)) {
        char c;
        c = (text[i]);
        text[i] = text[n-i-1];
        text[n-i-1] = c;
        i++;
    }
}
 
 
class Seq2Int {
public:
    Seq2Int(int s) : size(s) {}
    int operator() (const char* s) const {
    int erg=0;
    for( int i = 0; i < size; i++ ){
        erg <<= 2;
        erg |= base2int(s[i]);
    }
    return erg;
    }
    int rc(const char* s) const {
    int erg=0;
    for ( int i=0; i<size; i++) 
        erg |= base2int(wcComplement(s[i])) << (2*i);
    return erg;
    }
    int rev(const char* s) const {
    int erg=0;
    for ( int i=0; i<size; i++) 
        erg |= base2int(s[i]) << (2*i);
    return erg;
    }
    string inv(int pn) const {
    string result(size, 'n');
    for(int i = size-1; i>=0; i--) {
        result[i] = int2base(pn%4);
        pn >>= 2; // == k /= 4;
    }
    return result;
    }
    string INV(int pn) const {
    string result(size, 'N');
    for(int i = size-1; i>=0; i--) {
        result[i] = int2BASE(pn%4);
        pn >>= 2; // == k /= 4;
    }
    return result;
    }
    
    int read(istream& strm) const {
    char* buf = new char[size];
    // Read the character from the stream until the internal
    // representation size is reached, or a character
    // cannot be recognised.
    //------------------------------------------------------------
    for( int i = 0; i < size; i++ ){
        char c;
        strm.get(c);
        c = toupper(c);
        if( c == 'A' || c == 'C' || c == 'G' || c == 'T' )
        buf[i] = c;
    }
    int result = (*this)(buf);
    delete[] buf;
    return result;
    }
 
private:
    static int base2int(char c);
    static char int2base(int i);
    static char int2BASE(int i);
    int size;
};
 
inline int Seq2Int::base2int(char c) {
    switch (c) {
        case 'a': case 'A':
            return 0;
        case 'c': case 'C':
            return 1;
        case 'g': case 'G':
            return 2;
        case 't': case 'T':
            return 3;
        default:
        throw InvalidNucleotideError(c);
    }
}
 
inline char Seq2Int::int2base(int i) {
    switch (i) {
        case 0:
            return 'a';
        case 1:
            return 'c';
        case 2:
            return 'g';
        case 3:
            return 't';
        default:
            throw ProjectError("Seq2Int::int2base: internal error: i=" + itoa(i));
    }
}
 
inline char Seq2Int::int2BASE(int i) {
    switch (i) {
        case 0:
            return 'A';
        case 1:
            return 'C';
        case 2:
            return 'G';
        case 3:
            return 'T';
        default:
            throw ProjectError("Seq2Int::int2base: internal error: i=" + itoa(i));
    }
}
 
class ORF {
public:
    ORF() { start = end = -1; complete5prime = complete3prime = 0; strand = plusstrand; }
    int len() {return abs(end - start) + 1;}
    int start;
    int end;
    Strand strand;
    bool complete5prime;
    bool complete3prime;
};
 
class GeneticCode {
public:
     static void init() {
    chooseTranslationTable(1);
    }
    static void printReverseGeneticMap();
    
private:
    GeneticCode() {}
    static void reverseMap();
    static const char aa_symbols_with_stop[];
    static bool start_codons[];
    static Double start_codon_probs[];
    static const char * const TranslationTables[];
    static const char * const StartCodons[];
    static Seq2Int codon;
    static int translationtable;
    static int numStartCodons;
public:
    static void chooseTranslationTable( int );
    static const char* const aa_symbols;
    static const char* const aa_names[];
    static int map[];
    static int **syncodons;
    static int *codonsOfAA;
    static int get_aa_from_symbol(char c) {
    // this one returns -1 on '*' and string::npos-1 on invalid symbols
    // if (int)string::npos is -1, we get -2 here for invalid
    return string(GeneticCode::aa_symbols_with_stop).find(c)-1;
    }
    static char translate(int n) {
    return aa_symbols_with_stop[map[n]+1];
    }
    static char translate(const char* t);
    static char revtranslate(const char* t);
    static bool isStopcodon(const char* t) {
    return translate(t)=='*';
    }
    static bool isStartcodon(const char* t, bool rc=false){
    try {
        if (rc)
        return start_codons[codon.rc(t)];   
        else
        return start_codons[codon(t)];
    } catch (InvalidNucleotideError &e){}
    return false;
    }
    static bool isStartcodon(int pn){
    return start_codons[pn];
    }
    static Double startCodonProb(const char* t, bool rc=false){
    try {
        int pn = rc? codon.rc(t) : codon(t);
        if (start_codons[pn])
        return start_codon_probs[pn];
    } catch (InvalidNucleotideError &e){}
    return 0;
    }
    static Double startCodonProb(int pn){
    return start_codon_probs[pn];
    }
    static bool isRCStopcodon(const char* t) {
    return revtranslate(t)=='*';
    }
    static bool containsInFrameStopcodon(const char*, int, int, bool, int);
    static void printStartCodons();
    static void trainStartCodonProbs(int startcounts[]);
    static void writeStart(ofstream &out); // write start codon probs to file
    static void readStart(ifstream &in); // read start codon probs from file
    static bool is_purine(int b){
    return (b==0 || b==2); // 0=a, 1=c, 2=g, 3=t, a and g are purines
    }
    static ORF longestORF(const char* dna);
};
 
// getSampledCDS: sample a coding region from emission probabilities
char *getSampledCDS(vector<Double> *emiprobs, int k, int numCodons);
 
#endif    //  _GENETIC_CODE_HH