14 class SingleReadQ : public Read {
16 SingleReadQ() { readseq = qscore = ""; len = 0; }
17 SingleReadQ(const std::string& name, const std::string& readseq, const std::string& qscore) {
19 this->readseq = readseq;
20 this->qscore = qscore;
21 this->len = readseq.length();
24 bool read(int argc, std::istream* argv[], int flags = 7);
25 void write(int argc, std::ostream* argv[]);
27 int getReadLength() const { return len; }
28 const std::string& getReadSeq() const { return readseq; }
29 const std::string& getQScore() const { return qscore; }
31 void calc_lq(bool, int); // calculate if this read is low quality. Without calling this function, isLowQuality() will always be false
34 int len; // read length
35 std::string readseq, qscore; // qscore : quality scores
38 bool SingleReadQ::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 FASTQ file!\n"); exit(-1); }
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 = ""; }
57 void SingleReadQ::write(int argc, std::ostream* argv[]) {
59 (*argv[0])<<"@"<<name<<std::endl<<readseq<<std::endl<<"+\n"<<qscore<<std::endl;
62 //calculate if this read is low quality
63 void SingleReadQ::calc_lq(bool hasPolyA, int seedLen) {
65 if (len < seedLen) { low_quality = true; return; }
67 // if no polyA, no need to do the following calculation
68 if (!hasPolyA) return;
70 assert(readseq != "");
72 int numA = 0, numT = 0, numAO = 0, numTO = 0; // numAO : number of A in overlap seed region
73 int threshold_1, threshold_2;
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') {
80 if (i < OLEN) ++numAO;
82 if (readseq[i] == 'T') {
84 if (i >= len - OLEN) ++numTO;
88 if (numA >= threshold_1) {
89 low_quality = (numAO >= threshold_2);
91 else if (numT >= threshold_1) {
92 low_quality = (numTO >= threshold_2);
94 else low_quality = false;