SingleRead.h

   1 #ifndef SINGLEREAD
   2 #define SINGLEREAD
   3
   4 #include<cmath>
   5 #include<cstdio>
   6 #include<cstdlib>
   7 #include<cassert>
   8 #include<iostream>
   9 #include<string>
  10
  11 #include "utils.h"
  12 #include "Read.h"
  13
  14 class SingleRead : public Read {
  15  public:
  16   SingleRead() { readseq = ""; len = 0; }
  17   SingleRead(const std::string& name, const std::string& readseq) {
  18           this->name = name;
  19           this->readseq = readseq;
  20           this->len = readseq.length();
  21           calc_lq();
  22   }
  23
  24   bool read(int argc, std::istream* argv[], int flags = 7);
  25   void write(int argc, std::ostream* argv[]);
  26
  27   const int getReadLength() const { return len; /*readseq.length();*/ } // If need memory and .length() are guaranteed O(1), use statement in /* */
  28   const std::string& getReadSeq() const { return readseq; }
  29
  30  private:
  31   int len; // read length
  32   std::string readseq; // read sequence
  33
  34   void calc_lq();
  35 };
  36
  37 //If return false, you should not trust the value of any member
  38 bool SingleRead::read(int argc, std::istream* argv[], int flags) {
  39   std::string line;
  40
  41   assert(argc == 1);
  42   if (!getline((*argv[0]), line)) return false;
  43   if (line[0] != '>') { fprintf(stderr, "Read file does not look like a FASTA file!"); exit(-1); }
  44   name = "";
  45   if (flags & 4) { name = line.substr(1); }
  46   if (!getline((*argv[0]), readseq)) return false;
  47   len = readseq.length(); // set read length
  48   if (!(flags & 1)) { readseq = ""; }
  49
  50   if (flags & 1) calc_lq();
  51
  52   return true;
  53 }
  54
  55 void SingleRead::write(int argc, std::ostream* argv[]) {
  56   assert(argc == 1);
  57   (*argv[0])<<">"<<name<<std::endl<<readseq<<std::endl;
  58 }
  59
  60 void SingleRead::calc_lq() {
  61   int numA = 0, numT = 0, numAO = 0, numTO = 0; // numAO : number of A in overlap seed region
  62   int threshold_1, threshold_2;
  63
  64   if (len < OLEN) { low_quality = true; return; }
  65
  66   threshold_1 = int(0.9 * len - 1.5 * sqrt(len * 1.0) + 0.5);
  67   threshold_2 = (OLEN - 1) / 2 + 1;
  68   for (int i = 0; i < len; i++) {
  69     if (readseq[i] == 'A') {
  70       ++numA;
  71       if (i < OLEN) ++numAO;
  72     }
  73     if (readseq[i] == 'T') {
  74       ++numT;
  75       if (i >= len - OLEN) ++numTO;
  76     }
  77   }
  78
  79   if (numA >= threshold_1) {
  80     low_quality = (numAO >= threshold_2);
  81   }
  82   else if (numT >= threshold_1) {
  83     low_quality = (numTO >= threshold_2);
  84   }
  85   else low_quality = false;
  86 }
  87
  88 #endif