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