fixed bug where excessive 5-prime low qualities caused incorrect trimming

Author: najoshi

Makefile

@@ -1,5 +1,5 @@
 PROGRAM_NAME = sickle
-VERSION = 1.21
+VERSION = 1.29
 CC = gcc
 CFLAGS = -Wall -pedantic -DVERSION=$(VERSION)
 DEBUG = -g

src/sliding.c

@@ -95,13 +95,6 @@ cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int
 			for (j=window_start; j<window_start+window_size; j++) {
 				if (get_quality_num (fqrec->qual.s[j], qualtype, fqrec, j) < qual_threshold) {
 					three_prime_cut = j;
-
-					/* if cutting length is less than threshold then return -1 for both */
-					/* to indicate that the read should be discarded */
-					if (three_prime_cut - five_prime_cut < length_threshold) {
-						three_prime_cut = -1;
-						five_prime_cut = -1;
-					}
 					break;
 				}
 			}
@@ -117,12 +110,16 @@ cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int
 		window_start++;
 	}
 
-    /* If you never find a five prime cut site, then discard whole read */
-    if (found_five_prime == 0) {
+
+    /* if cutting length is less than threshold then return -1 for both */
+    /* to indicate that the read should be discarded */
+    /* Also, if you never find a five prime cut site, then discard whole read */
+    if ((found_five_prime == 0) || (three_prime_cut - five_prime_cut < length_threshold)) {
         three_prime_cut = -1;
         five_prime_cut = -1;
     }
 
+
 	retvals = (cutsites*) malloc (sizeof(cutsites));
 	retvals->three_prime_cut = three_prime_cut;
 	retvals->five_prime_cut = five_prime_cut;

src/trim_paired.c

@@ -324,6 +324,9 @@ int paired_main (int argc, char *argv[]) {
     p1cut = sliding_window (fqrec1, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, discard_n);
     p2cut = sliding_window (fqrec2, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, discard_n);
 
+    if (debug) {printf ("p1cut: %d,%d\n", p1cut->five_prime_cut, p1cut->three_prime_cut);}
+    if (debug) {printf ("p2cut: %d,%d\n", p2cut->five_prime_cut, p2cut->three_prime_cut);}
+
     /* The sequence and quality print statements below print out the sequence string starting from the 5' cut */
     /* and then only print out to the 3' cut, however, we need to adjust the 3' cut */
     /* by subtracting the 5' cut because the 3' cut was calculated on the original sequence */

src/trim_single.c

@@ -170,6 +170,8 @@ int single_main (int argc, char *argv[]) {
 
 		p1cut = sliding_window (fqrec, qualtype, single_length_threshold, single_qual_threshold, no_fiveprime, discard_n);
 
+        if (debug) {printf ("P1cut: %d,%d\n", p1cut->five_prime_cut, p1cut->three_prime_cut);}
+
 		/* if sequence quality and length pass filter then output record, else discard */
 		if (p1cut->three_prime_cut >= 0) {
 			fprintf (outfile, "@%s", fqrec->name.s);