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
  24         bool read(int argc, std::istream* argv[], int flags = 7);
  25         void write(int argc, std::ostream* argv[]);
  26
  27         int getReadLength() const { return len; }
  28         const std::string& getReadSeq() const { return readseq; }
  29         const std::string& getQScore() const { return qscore; }
  30
  31         void calc_lq(bool, int); // calculate if this read is low quality. Without calling this function, isLowQuality() will always be false
  32
  33 private:
  34         int len; // read length
  35         std::string readseq, qscore; // qscore : quality scores
  36 };
  37
  38 bool SingleReadQ::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 FASTQ file!\n"); exit(-1); }
  44         name = "";
  45         if (flags & 4) { name = line.substr(1); }
  46         if (!getline((*argv[0]), readseq)) return false;
  47         len = readseq.length();
  48         if (!(flags & 1)) { readseq = ""; }
  49         if (!getline((*argv[0]), line)) return false;
  50         if (line[0] != '+') { fprintf(stderr, "Read file does not look like a FASTQ file!\n"); exit(-1); }
  51         if (!getline((*argv[0]), qscore)) return false;
  52         if (!(flags & 2)) { qscore = ""; }
  53
  54         return true;
  55 }
  56
  57 void SingleReadQ::write(int argc, std::ostream* argv[]) {
  58         assert(argc == 1);
  59         (*argv[0])<<"@"<<name<<std::endl<<readseq<<std::endl<<"+\n"<<qscore<<std::endl;
  60 }
  61
  62 //calculate if this read is low quality
  63 void SingleReadQ::calc_lq(bool hasPolyA, int seedLen) {
  64         low_quality = false;
  65         if (len < seedLen) { low_quality = true; return; }
  66
  67         // if no polyA, no need to do the following calculation
  68         if (!hasPolyA) return;
  69
  70         assert(readseq != "");
  71
  72         int numA = 0, numT = 0, numAO = 0, numTO = 0; // numAO : number of A in overlap seed region
  73         int threshold_1, threshold_2;
  74
  75         threshold_1 = int(0.9 * len - 1.5 * sqrt(len * 1.0) + 0.5);
  76         threshold_2 = (OLEN - 1) / 2 + 1;
  77         for (int i = 0; i < len; i++) {
  78                 if (readseq[i] == 'A') {
  79                         ++numA;
  80                         if (i < OLEN) ++numAO;
  81                 }
  82                 if (readseq[i] == 'T') {
  83                         ++numT;
  84                         if (i >= len - OLEN) ++numTO;
  85                 }
  86         }
  87
  88         if (numA >= threshold_1) {
  89                 low_quality = (numAO >= threshold_2);
  90         }
  91         else if (numT >= threshold_1) {
  92                 low_quality = (numTO >= threshold_2);
  93         }
  94         else low_quality = false;
  95 }
  96
  97 #endif