Refs.h

   1 #ifndef REFS
   2 #define REFS
   3
   4 #include<cstdio>
   5 #include<cstring>
   6 #include<cstdlib>
   7 #include<cassert>
   8 #include<string>
   9 #include<fstream>
  10 #include<vector>
  11
  12 #include "utils.h"
  13 #include "RefSeq.h"
  14 #include "RefSeqPolicy.h"
  15 #include "PolyARules.h"
  16
  17
  18 class Refs {
  19  public:
  20   Refs() {
  21     M = 0;
  22     seqs.clear();
  23     has_polyA = false;
  24   }
  25
  26   ~Refs() {}
  27
  28   void makeRefs(char*, RefSeqPolicy&, PolyARules&);
  29   void loadRefs(char*, int = 0);
  30   void saveRefs(char*);
  31
  32   int getM() { return M; } // get number of isoforms
  33
  34   //int getNS() { return M + 1; } // get number of parameters, I do not think we need it here.
  35
  36   RefSeq& getRef(int sid) { return seqs[sid]; } // get a particular reference
  37
  38   std::vector<RefSeq>& getRefs() { return seqs; } // may be slow, for copying the whole thing
  39
  40   bool hasPolyA() { return has_polyA; } // if any of sequence has poly(A) tail added
  41
  42   //lim : >=0 If mismatch > lim , return; -1 find all mismatches
  43   int countMismatch(const std::string& seq, int pos, const std::string& readseq, int LEN, int lim = -1) {
  44     int nMis = 0; // number of mismatches
  45
  46     for (int i = 0; i < LEN; i++) {
  47       char rc = toupper(readseq[i]);
  48       if (seq[i + pos] == 'N' || rc == 'N' || seq[i + pos] != rc) nMis++;
  49
  50       // a speed up tech
  51       if (lim >= 0 && nMis > lim) return nMis;
  52     }
  53
  54     return nMis;
  55   }
  56
  57   bool isValid(int sid, int dir, int pos, const std::string& readseq, int LEN, int C) {
  58     if (sid <= 0 || sid > M || (dir != 0 && dir != 1) ||  pos < 0 || pos + LEN > seqs[sid].getTotLen() || LEN > (int)readseq.length()) return false;
  59     const std::string& seq = seqs[sid].getSeq(dir);
  60     return countMismatch(seq, pos, readseq, LEN, C) <= C;
  61   }
  62
  63   // get segment from refs
  64   std::string getSegment(int sid, int dir, int pos, int LEN) {
  65     if (pos < 0 || pos + LEN > seqs[sid].getTotLen()) return "fail";
  66
  67     const std::string& seq = seqs[sid].getSeq(dir);
  68     std::string seg = "";
  69
  70     for (int i = 0; i < LEN; i++)
  71       seg.append(1,  seq[pos + i]);
  72
  73     return seg;
  74   }
  75
  76  private:
  77   int M; // # of isoforms, id starts from 1
  78   std::vector<RefSeq> seqs;  // reference sequences, starts from 1; 0 is for noise gene
  79   bool has_polyA; // if at least one sequence has polyA added, the value is true; otherwise, the value is false
  80 };
  81
  82 //inpF in fasta format
  83 void Refs::makeRefs(char *inpF, RefSeqPolicy& policy, PolyARules& rules) {
  84   //read standard fasta format here
  85   std::ifstream fin;
  86   std::string tag, line, rawseq;
  87   void* pt; // istream& is indeed a pointer, that's why I can use void* here
  88
  89   seqs.clear();
  90   seqs.push_back(RefSeq()); // noise isoform
  91
  92   M = 0;
  93   has_polyA = false;
  94
  95   fin.open(inpF);
  96   if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
  97   pt = getline(fin, line);
  98   while (pt != 0 && line[0] == '>') {
  99     tag = line.substr(1);
 100     rawseq = "";
 101     while((pt = getline(fin, line)) && line[0] != '>') {
 102       rawseq += line;
 103     }
 104     if (rawseq.size() <= 0) {
 105       printf("Warning: Fasta entry %s has an empty sequence! It is omitted!\n", tag.c_str());
 106       continue;
 107     }
 108     ++M;
 109     seqs.push_back(RefSeq(tag, policy.convert(rawseq), rules.getLenAt(tag)));
 110     has_polyA = has_polyA || seqs[M].getFullLen() < seqs[M].getTotLen();
 111   }
 112   fin.close();
 113
 114   if (verbose) { printf("Refs.makeRefs finished!\n"); }
 115 }
 116
 117 //inpF in fasta format, with sequence all in one line together
 118 //option 0 read all, 1 do not read sequences
 119 void Refs::loadRefs(char *inpF, int option) {
 120   std::ifstream fin;
 121   RefSeq seq;
 122
 123   fin.open(inpF);
 124   if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
 125   seqs.clear();
 126   seqs.push_back(RefSeq());
 127
 128   M = 0;
 129   has_polyA = false;
 130
 131   bool success;
 132   do {
 133     success = seq.read(fin, option);
 134     if (success) {
 135         seqs.push_back(seq);
 136         ++M;
 137         has_polyA = has_polyA || seq.getFullLen() < seq.getTotLen();
 138     }
 139   } while (success);
 140
 141   fin.close();
 142
 143   assert(M + 1 == (int)seqs.size());
 144
 145   if (verbose) { printf("Refs.loadRefs finished!\n"); }
 146 }
 147
 148 void Refs::saveRefs(char* outF) {
 149   std::ofstream fout;
 150
 151   fout.open(outF);
 152   for (int i = 1; i <= M; i++) {
 153     seqs[i].write(fout);
 154   }
 155   fout.close();
 156
 157   if (verbose) { printf("Refs.saveRefs finished!\n"); }
 158 }
 159
 160 #endif