14 class SingleRead : public Read {
16 SingleRead() { readseq = ""; len = 0; }
17 SingleRead(const std::string& name, const std::string& readseq) {
19 this->readseq = readseq;
20 this->len = readseq.length();
24 bool read(int argc, std::istream* argv[], int flags = 7);
25 void write(int argc, std::ostream* argv[]);
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; }
31 int len; // read length
32 std::string readseq; // read sequence
37 //If return false, you should not trust the value of any member
38 bool SingleRead::read(int argc, std::istream* argv[], int flags) {
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); }
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 = ""; }
50 if (flags & 1) calc_lq();
55 void SingleRead::write(int argc, std::ostream* argv[]) {
57 (*argv[0])<<">"<<name<<std::endl<<readseq<<std::endl;
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;
64 threshold_1 = int(0.9 * len - 1.5 * sqrt(len * 1.0) + 0.5);
65 threshold_2 = (OLEN - 1) / 2 + 1;
66 for (int i = 0; i < len; i++) {
67 if (readseq[i] == 'A') {
69 if (i < OLEN) ++numAO;
71 if (readseq[i] == 'T') {
73 if (i >= len - OLEN) ++numTO;
77 if (numA >= threshold_1) {
78 low_quality = (numAO >= threshold_2);
80 else if (numT >= threshold_1) {
81 low_quality = (numTO >= threshold_2);
83 else low_quality = false;