Sequence_AMOS.hh

Go to the documentation of this file.

#ifndef __Sequence_AMOS_HH
#define __Sequence_AMOS_HH 1

#include "Universal_AMOS.hh"
#include <string>




namespace AMOS {

//================================================ Sequence_t ==================
//==============================================================================
class Sequence_t : public Universal_t
{

protected:

  uint8_t * seq_m;      
  uint8_t * qual_m;     
  Size_t length_m;      


  static const uint8_t COMPRESS_BIT  = 0x1;   
  static const uint8_t ADENINE_BITS  = 0x0;   
  static const uint8_t CYTOSINE_BITS = 0x40;  
  static const uint8_t GUANINE_BITS  = 0x80;  
  static const uint8_t THYMINE_BITS  = 0xC0;  
  static const uint8_t SEQ_BITS      = 0xC0;  
  static const uint8_t QUAL_BITS     = 0x3F;  


  //--------------------------------------------------- compress ---------------
  static inline uint8_t compress (char seqchar, char qualchar)
  {
    //-- Force quality score into its bits
    qualchar = Char2Qual (qualchar);

    if ( qualchar & SEQ_BITS )
      return 0;

    switch ( toupper (seqchar) )
      {
      case 'A': return (uint8_t)qualchar | ADENINE_BITS;
      case 'C': return (uint8_t)qualchar | CYTOSINE_BITS;
      case 'G': return (uint8_t)qualchar | GUANINE_BITS;
      case 'T': return (uint8_t)qualchar | THYMINE_BITS;
      case 'N': return 0;
      default:
        return 0;
      }
  }


  //--------------------------------------------------- uncompress -------------
  static inline std::pair<char, char> uncompress (uint8_t byte)
  {
    std::pair<char, char> retval;

    switch ( byte & SEQ_BITS )
      {
      case ADENINE_BITS:  retval . first = 'A'; break;
      case CYTOSINE_BITS: retval . first = 'C'; break;
      case GUANINE_BITS:  retval . first = 'G'; break;
      case THYMINE_BITS:  retval . first = 'T'; break;
      }

    byte &= QUAL_BITS;
    if ( byte == 0 )
      retval . first = 'N';

    retval . second = Qual2Char (byte);

    return retval;
  }


  //--------------------------------------------------- readRecord -------------
  virtual void readRecord (std::istream & fix, std::istream & var);

  //--------------------------------------------------- readRecordFix ----------
  virtual void readRecordFix (std::istream & fix);

  //--------------------------------------------------- writeRecord ------------
  virtual void writeRecord (std::ostream & fix, std::ostream & var) const;


public:

  static const NCode_t NCODE;


  //--------------------------------------------------- Sequence_t -------------
  Sequence_t ( )
  {
    seq_m = qual_m = NULL;
    length_m = 0;
  }


  //--------------------------------------------------- Sequence_t -------------
  Sequence_t (const Sequence_t & source)
  {
    seq_m = qual_m = NULL;
    *this = source;
  }


  //--------------------------------------------------- ~Sequence_t ------------
  virtual ~Sequence_t ( )
  {
    free (seq_m);
    free (qual_m);
  }


  //--------------------------------------------------- clear ------------------
  virtual void clear ( );


  //--------------------------------------------------- compress ---------------
  void compress ( );


  //--------------------------------------------------- getBase ----------------
  std::pair<char, char> getBase (Pos_t index) const
  {
    if (seq_m == NULL)
      AMOS_THROW_ARGUMENT("No sequence data");
    
    if ( index < 0 || index >= length_m )
      AMOS_THROW_ARGUMENT ("Requested sequence index is out of range");

    if ( isCompressed( ) )
      return uncompress (seq_m [index]);
    else
      return std::make_pair ((char)(seq_m [index]), (char)(qual_m [index]));
  }



  //--------------------------------------------------- getGCContent --------------
  virtual double getGCContent (const Range_t & rng) const
  {
    int gc = 0;
    int all = 0;

    int start = rng.getLo();
    int stop = rng.getHi();

    // skip ambiguities, gaps

    for (int i = start; i < stop; i++)
    {
      switch(getBase(i).first)
      {
        case 'A':
        case 'a': 
        case 'T':
        case 't': all++; break;

        case 'C':
        case 'c':
        case 'G':
        case 'g': all++; gc++; break;
      };
    }

    return (all) ? ((double)gc)/all : 0.0;
  }

  //--------------------------------------------------- getGCContent --------------
  virtual double getGCContent ( ) const
  {
    return getGCContent(Range_t(0, length_m));
  }



  //--------------------------------------------------- getLength --------------
  Size_t getLength ( ) const
  {
    return length_m;
  }


  //--------------------------------------------------- getNCode ---------------
  virtual NCode_t getNCode ( ) const
  {
    return Sequence_t::NCODE;
  }


  //--------------------------------------------------- getQualString ----------
  std::string getQualString ( ) const
  {
    return getQualString (Range_t (0, length_m));
  }


  //--------------------------------------------------- getQualString ----------
  std::string getQualString (Range_t range) const;


  //--------------------------------------------------- getSeqString -----------
  std::string getSeqString ( ) const
  {
    return getSeqString (Range_t (0, length_m));
  }


  //--------------------------------------------------- getSeqString -----------
  std::string getSeqString (Range_t range) const;


  //--------------------------------------------------- isCompressed -----------
  bool isCompressed ( ) const
  {
    return flags_m . nibble & COMPRESS_BIT;
  }


  //--------------------------------------------------- readMessage ------------
  virtual void readMessage (const Message_t & msg);


  //--------------------------------------------------- setBase ----------------
  void setBase (char seqchar, char qualchar, Pos_t index)
  {  
    if (seq_m == NULL)
      AMOS_THROW_ARGUMENT("No sequence data");

    if ( index < 0 || index >= length_m )
      AMOS_THROW_ARGUMENT ("Requested sequence index is out of range");

    if ( isCompressed( ) )
      seq_m [index] = compress (seqchar, qualchar);
    else
      {
        seq_m  [index] = seqchar;
        qual_m [index] = qualchar;
      }
  }


  //--------------------------------------------------- setSequence ------------
  void setSequence (const char * seq, const char * qual);


  //--------------------------------------------------- setSequence ------------
  void setSequence (const std::string & seq, const std::string & qual);


  //--------------------------------------------------- uncompress -------------
  void uncompress ( );


  //--------------------------------------------------- operator= --------------
  Sequence_t & operator= (const Sequence_t & source);


  //--------------------------------------------------- writeMessage -----------
  virtual void writeMessage (Message_t & msg) const;

};

} // namespace AMOS

#endif // #ifndef __Sequence_AMOS_HH

---