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         }
  22
  23         bool read(int argc, std::istream* argv[], int flags = 7);
  24         void write(int argc, std::ostream* argv[]);
  25
  26         const int getReadLength() const { return len; /*readseq.length();*/ } // If need memory and .length() are guaranteed O(1), use statement in /* */
  27         const std::string& getReadSeq() const { return readseq; }
  28
  29         void calc_lq(bool, int); // calculate if this read is low quality. Without calling this function, isLowQuality() will always be false
  30
  31 private:
  32         int len; // read length
  33         std::string readseq; // read sequence
  34 };
  35
  36 //If return false, you should not trust the value of any member
  37 bool SingleRead::read(int argc, std::istream* argv[], int flags) {
  38         std::string line;
  39
  40         assert(argc == 1);
  41         if (!getline((*argv[0]), line)) return false;
  42         if (line[0] != '>') { fprintf(stderr, "Read file does not look like a FASTA file!"); exit(-1); }
  43         name = "";
  44         if (flags & 4) { name = line.substr(1); }
  45         if (!getline((*argv[0]), readseq)) return false;
  46         len = readseq.length(); // set read length
  47         if (!(flags & 1)) { readseq = ""; }
  48
  49         return true;
  50 }
  51
  52 void SingleRead::write(int argc, std::ostream* argv[]) {
  53         assert(argc == 1);
  54         (*argv[0])<<">"<<name<<std::endl<<readseq<<std::endl;
  55 }
  56
  57 //calculate if this read is low quality
  58 void SingleRead::calc_lq(bool hasPolyA, int seedLen) {
  59         low_quality = false;
  60         if (len < seedLen) { low_quality = true; return; }
  61
  62         // if no polyA, no need to do the following calculation
  63         if (!hasPolyA) return;
  64
  65         assert(readseq != "");
  66
  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         threshold_1 = int(0.9 * len - 1.5 * sqrt(len * 1.0) + 0.5);
  71         threshold_2 = (OLEN - 1) / 2 + 1;
  72         for (int i = 0; i < len; i++) {
  73                 if (readseq[i] == 'A') {
  74                         ++numA;
  75                         if (i < OLEN) ++numAO;
  76                 }
  77                 if (readseq[i] == 'T') {
  78                         ++numT;
  79                         if (i >= len - OLEN) ++numTO;
  80                 }
  81         }
  82
  83         if (numA >= threshold_1) {
  84                 low_quality = (numAO >= threshold_2);
  85         }
  86         else if (numT >= threshold_1) {
  87                 low_quality = (numTO >= threshold_2);
  88         }
  89         else low_quality = false;
  90 }
  91
  92 #endif