SingleReadQ.h

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