RefSeq.h

   1 #ifndef REFSEQ
   2 #define REFSEQ
   3
   4 #include<cassert>
   5 #include<fstream>
   6 #include<string>
   7 #include<vector>
   8
   9 #include "utils.h"
  10
  11 //Each Object can only be used once
  12 class RefSeq {
  13  public:
  14
  15   RefSeq() {
  16     fullLen = totLen = 0;
  17     name = ""; seq = "";
  18     fmasks.clear();
  19   }
  20
  21   //Constructor , seq : the forward strand of the reference
  22   //tag does not contain ">"
  23   //polyALen : length of polyA tail we add
  24   RefSeq(const std::string& name, const std::string& seq, int polyALen) {
  25     fullLen = seq.length();
  26     totLen = fullLen + polyALen;
  27     this->name = name;
  28     this->seq = seq;
  29     this->seq.append(polyALen, 'A');
  30
  31     assert(fullLen > 0 && totLen >= fullLen);
  32
  33     int len = (fullLen - 1) / NBITS + 1;
  34     fmasks.clear(); fmasks.resize(len, 0);
  35     // ask read to be at least OLEN long!
  36     for (int i = std::max(fullLen - OLEN + 1, 0); i < fullLen; i++) setMask(i);
  37   }
  38
  39   RefSeq(const RefSeq& o) {
  40     fullLen = o.fullLen;
  41     totLen = o.totLen;
  42     name = o.name;
  43     seq = o.seq;
  44     fmasks = o.fmasks;
  45   }
  46
  47   RefSeq& operator= (const RefSeq &rhs) {
  48     if (this != &rhs) {
  49       fullLen = rhs.fullLen;
  50       totLen = rhs.totLen;
  51       name = rhs.name;
  52       seq = rhs.seq;
  53       fmasks = rhs.fmasks;
  54     }
  55
  56     return *this;
  57   }
  58
  59   ~RefSeq() {
  60   }
  61
  62   bool read(std::ifstream&, int  = 0);
  63   void write(std::ofstream&);
  64
  65   int getFullLen() const { return fullLen; }
  66
  67   int getTotLen() const { return totLen; }
  68
  69   const std::string& getName() const { return name; }
  70
  71   std::string getSeq() const { return seq; }
  72
  73   std::string getRSeq() const {
  74     std::string rseq = "";
  75     for (int i = totLen - 1; i >= 0; i--) rseq.push_back(getCharacter(get_rbase_id(seq[i])));
  76     return rseq;
  77   }
  78
  79   //get the sequence  dir 0 : + 1 : -
  80   std::string getSeq(int dir) const {
  81     return (dir == 0 ? getSeq() : getRSeq());
  82   }
  83
  84
  85   int get_id(int pos, int dir) const {
  86     assert(pos >= 0 && pos < totLen);
  87     return (dir == 0 ? get_base_id(seq[pos]) : get_rbase_id(seq[totLen - pos - 1]));
  88   }
  89
  90   bool getMask(int seedPos) const {
  91     assert(seedPos >= 0 && seedPos < totLen);
  92     return fmasks[seedPos / NBITS] & mask_codes[seedPos % NBITS];
  93   }
  94
  95   void setMask(int seedPos) {
  96     assert(seedPos >= 0 && seedPos < totLen);
  97     fmasks[seedPos / NBITS] |= mask_codes[seedPos % NBITS];
  98   }
  99
 100  private:
 101   int fullLen; // fullLen : the original length of an isoform
 102   int totLen; // totLen : the total length, included polyA tails, if any
 103   std::string name; // the tag
 104   std::string seq; // the raw sequence, in forward strand
 105   std::vector<unsigned int> fmasks; // record masks for forward strand, each position occupies 1 bit
 106 };
 107
 108 //internal read; option 0 : read all 1 : do not read seqence and name
 109 bool RefSeq::read(std::ifstream& fin, int option) {
 110   std::string line;
 111
 112   if (!(fin>>fullLen>>totLen)) return false;
 113   assert(fullLen > 0 && totLen >= fullLen);
 114   getline(fin, line);
 115   if (!getline(fin, name)) return false;
 116   if (!getline(fin, seq)) return false;
 117
 118   int len = (fullLen - 1) / NBITS + 1; // assume each cell contains NBITS bits
 119   fmasks.resize(len, 0);
 120   for (int i = 0; i < len; i++)
 121     if (!(fin>>fmasks[i])) return false;
 122   getline(fin, line);
 123
 124   assert(option == 0 || option == 1);
 125   if (option == 1) { name = seq = ""; }
 126
 127   return true;
 128 }
 129
 130 //write to file in "internal" format
 131 void RefSeq::write(std::ofstream& fout) {
 132   fout<<fullLen<<" "<<totLen<<std::endl;
 133   fout<<name<<std::endl;
 134   fout<<seq<<std::endl;
 135
 136   int len = fmasks.size();
 137   for (int i = 0; i < len - 1; i++) fout<<fmasks[i]<<" ";
 138   fout<<fmasks[len - 1]<<std::endl;
 139 }
 140
 141 #endif