-void SingleReadQ::calc_lq() {
- int numA = 0, numT = 0, numAO = 0, numTO = 0; // numAO : number of A in overlap seed region
- int threshold_1, threshold_2;
-
- if (len < OLEN) { low_quality = true; return; }
-
- threshold_1 = int(0.9 * len - 1.5 * sqrt(len * 1.0) + 0.5);
- threshold_2 = (OLEN - 1) / 2 + 1;
- for (int i = 0; i < len; i++) {
- if (readseq[i] == 'A') {
- ++numA;
- if (i < OLEN) ++numAO;
- }
- if (readseq[i] == 'T') {
- ++numT;
- if (i >= len - OLEN) ++numTO;
- }
- }
-
- if (numA >= threshold_1) {
- low_quality = (numAO >= threshold_2);
- }
- else if (numT >= threshold_1) {
- low_quality = (numTO >= threshold_2);
- }
- else low_quality = false;
+//calculate if this read is low quality
+void SingleReadQ::calc_lq(bool hasPolyA, int seedLen) {
+ low_quality = false;
+ if (len < seedLen) { low_quality = true; return; }
+
+ // if no polyA, no need to do the following calculation
+ if (!hasPolyA) return;
+
+ assert(readseq != "");
+
+ int numA = 0, numT = 0, numAO = 0, numTO = 0; // numAO : number of A in overlap seed region
+ int threshold_1, threshold_2;
+
+ threshold_1 = int(0.9 * len - 1.5 * sqrt(len * 1.0) + 0.5);
+ threshold_2 = (OLEN - 1) / 2 + 1;
+ for (int i = 0; i < len; i++) {
+ if (readseq[i] == 'A') {
+ ++numA;
+ if (i < OLEN) ++numAO;
+ }
+ if (readseq[i] == 'T') {
+ ++numT;
+ if (i >= len - OLEN) ++numTO;
+ }
+ }
+
+ if (numA >= threshold_1) {
+ low_quality = (numAO >= threshold_2);
+ }
+ else if (numT >= threshold_1) {
+ low_quality = (numTO >= threshold_2);
+ }
+ else low_quality = false;