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
  88   seqs.clear();
  89   seqs.push_back(RefSeq()); // noise isoform
  90
  91   M = 0;
  92   has_polyA = false;
  93
  94   fin.open(inpF);
  95   if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
  96   getline(fin, line);
  97   while ((fin) && (line[0] == '>')) {
  98     tag = line.substr(1);
  99     rawseq = "";
 100     while((getline(fin, line)) && (line[0] != '>')) {
 101       rawseq += line;
 102     }
 103     if (rawseq.size() <= 0) {
 104       printf("Warning: Fasta entry %s has an empty sequence! It is omitted!\n", tag.c_str());
 105       continue;
 106     }
 107     ++M;
 108     seqs.push_back(RefSeq(tag, policy.convert(rawseq), rules.getLenAt(tag)));
 109     has_polyA = has_polyA || seqs[M].getFullLen() < seqs[M].getTotLen();
 110   }
 111   fin.close();
 112
 113   if (verbose) { printf("Refs.makeRefs finished!\n"); }
 114 }
 115
 116 //inpF in fasta format, with sequence all in one line together
 117 //option 0 read all, 1 do not read sequences
 118 void Refs::loadRefs(char *inpF, int option) {
 119   std::ifstream fin;
 120   RefSeq seq;
 121
 122   fin.open(inpF);
 123   if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", inpF); exit(-1); }
 124   seqs.clear();
 125   seqs.push_back(RefSeq());
 126
 127   M = 0;
 128   has_polyA = false;
 129
 130   bool success;
 131   do {
 132     success = seq.read(fin, option);
 133     if (success) {
 134         seqs.push_back(seq);
 135         ++M;
 136         has_polyA = has_polyA || seq.getFullLen() < seq.getTotLen();
 137     }
 138   } while (success);
 139
 140   fin.close();
 141
 142   assert(M + 1 == (int)seqs.size());
 143
 144   if (verbose) { printf("Refs.loadRefs finished!\n"); }
 145 }
 146
 147 void Refs::saveRefs(char* outF) {
 148   std::ofstream fout;
 149
 150   fout.open(outF);
 151   for (int i = 1; i <= M; i++) {
 152     seqs[i].write(fout);
 153   }
 154   fout.close();
 155
 156   if (verbose) { printf("Refs.saveRefs finished!\n"); }
 157 }
 158
 159 #endif