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         RefSeq() {
  15                 fullLen = totLen = 0;
  16                 name = ""; seq = "";
  17                 fmasks.clear();
  18         }
  19
  20         //Constructor , seq : the forward strand of the reference
  21         //tag does not contain ">"
  22         //polyALen : length of polyA tail we add
  23         RefSeq(const std::string& name, const std::string& seq, int polyALen) {
  24                 fullLen = seq.length();
  25                 totLen = fullLen + polyALen;
  26                 this->name = name;
  27                 this->seq = seq;
  28                 this->seq.append(polyALen, 'A');
  29
  30                 assert(fullLen > 0 && totLen >= fullLen);
  31
  32                 int len = (fullLen - 1) / NBITS + 1;
  33                 fmasks.assign(len, 0);
  34                 // set mask if poly(A) tail is added
  35                 if (polyALen > 0) {
  36                         for (int i = std::max(fullLen - OLEN + 1, 0); i < fullLen; i++) setMask(i);
  37                 }
  38   }
  39
  40         RefSeq(const RefSeq& o) {
  41                 fullLen = o.fullLen;
  42                 totLen = o.totLen;
  43                 name = o.name;
  44                 seq = o.seq;
  45                 fmasks = o.fmasks;
  46         }
  47
  48         RefSeq& operator= (const RefSeq &rhs) {
  49                 if (this != &rhs) {
  50                         fullLen = rhs.fullLen;
  51                         totLen = rhs.totLen;
  52                         name = rhs.name;
  53                         seq = rhs.seq;
  54                         fmasks = rhs.fmasks;
  55                 }
  56
  57                 return *this;
  58         }
  59
  60         ~RefSeq() {}
  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         int get_id(int pos, int dir) const {
  85                 assert(pos >= 0 && pos < totLen);
  86                 return (dir == 0 ? get_base_id(seq[pos]) : get_rbase_id(seq[totLen - pos - 1]));
  87         }
  88
  89         bool getMask(int seedPos) const {
  90                 assert(seedPos >= 0 && seedPos < totLen);
  91                 return fmasks[seedPos / NBITS] & mask_codes[seedPos % NBITS];
  92         }
  93
  94         void setMask(int seedPos) {
  95                 assert(seedPos >= 0 && seedPos < totLen);
  96                 fmasks[seedPos / NBITS] |= mask_codes[seedPos % NBITS];
  97         }
  98
  99 private:
 100         int fullLen; // fullLen : the original length of an isoform
 101         int totLen; // totLen : the total length, included polyA tails, if any
 102         std::string name; // the tag
 103         std::string seq; // the raw sequence, in forward strand
 104         std::vector<unsigned int> fmasks; // record masks for forward strand, each position occupies 1 bit
 105 };
 106
 107 //internal read; option 0 : read all 1 : do not read seqences
 108 bool RefSeq::read(std::ifstream& fin, int option) {
 109         std::string line;
 110
 111         if (!(fin>>fullLen>>totLen)) return false;
 112         assert(fullLen > 0 && totLen >= fullLen);
 113         getline(fin, line);
 114         if (!getline(fin, name)) return false;
 115         if (!getline(fin, seq)) return false;
 116
 117         int len = (fullLen - 1) / NBITS + 1; // assume each cell contains NBITS bits
 118         fmasks.assign(len, 0);
 119         for (int i = 0; i < len; i++)
 120             if (!(fin>>fmasks[i])) return false;
 121         getline(fin, line);
 122
 123         assert(option == 0 || option == 1);
 124         if (option == 1) { seq = ""; }
 125
 126         return true;
 127 }
 128
 129 //write to file in "internal" format
 130 void RefSeq::write(std::ofstream& fout) {
 131         fout<<fullLen<<" "<<totLen<<std::endl;
 132         fout<<name<<std::endl;
 133         fout<<seq<<std::endl;
 134
 135         int len = fmasks.size();
 136         for (int i = 0; i < len - 1; i++) fout<<fmasks[i]<<" ";
 137         fout<<fmasks[len - 1]<<std::endl;
 138 }
 139
 140 #endif