fixed bug where disabling 5-prime trimming discarded all reads

najoshi · najoshi · commit 8bec85654e36 · 2014-07-22T19:03:28.000-07:00
diff --git a/src/sickle.h b/src/sickle.h
@@ -96,6 +96,6 @@ typedef struct __cutsites_ {
 /* Function Prototypes */
 int single_main (int argc, char *argv[]);
 int paired_main (int argc, char *argv[]);
-cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int qual_threshold, int no_fiveprime, int discard_n);
+cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int qual_threshold, int no_fiveprime, int trunc_n, int debug);
 
 #endif /*SICKLE_H*/
diff --git a/src/sliding.c b/src/sliding.c
@@ -32,7 +32,7 @@ int get_quality_num (char qualchar, int qualtype, kseq_t *fqrec, int pos) {
 }
 
 
-cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int qual_threshold, int no_fiveprime, int trunc_n) {
+cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int qual_threshold, int no_fiveprime, int trunc_n, int debug) {
 
 	int window_size = (int) (0.1 * fqrec->seq.l);
 	int i,j;
@@ -65,10 +65,14 @@ cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int
 
 		window_avg = (double)window_total / (double)window_size;
 
+        if (debug) printf ("no_fiveprime: %d, found 5prime: %d, window_avg: %f\n", no_fiveprime, found_five_prime, window_avg);
+
 		/* Finding the 5' cutoff */
 		/* Find when the average quality in the window goes above the threshold starting from the 5' end */
 		if (no_fiveprime == 0 && found_five_prime == 0 && window_avg >= qual_threshold) {
 
+        if (debug) printf ("inside 5-prime cut\n");
+
 			/* at what point in the window does the quality go above the threshold? */
 			for (j=window_start; j<window_start+window_size; j++) {
 				if (get_quality_num (fqrec->qual.s[j], qualtype, fqrec, j) >= qual_threshold) {
@@ -77,14 +81,16 @@ cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int
 				}
 			}
 
+            if (debug) printf ("five_prime_cut: %d\n", five_prime_cut);
+
 			found_five_prime = 1;
 		}
 
 		/* Finding the 3' cutoff */
 		/* if the average quality in the window is less than the threshold */
 		/* or if the window is the last window in the read */
 		if ((window_avg < qual_threshold || 
-			window_start+window_size > fqrec->qual.l) && found_five_prime == 1) {
+			window_start+window_size > fqrec->qual.l) && (found_five_prime == 1 || no_fiveprime)) {
 
 			/* at what point in the window does the quality dip below the threshold? */
 			for (j=window_start; j<window_start+window_size; j++) {
@@ -115,11 +121,14 @@ cutsites* sliding_window (kseq_t *fqrec, int qualtype, int length_threshold, int
     /* 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)) {
+    if ((found_five_prime == 0 && !no_fiveprime) || (three_prime_cut - five_prime_cut < length_threshold)) {
         three_prime_cut = -1;
         five_prime_cut = -1;
+
+        if (debug) printf("%s\n", fqrec->name.s);
     }
 
+    if (debug) printf ("\n\n");
 
 	retvals = (cutsites*) malloc (sizeof(cutsites));
 	retvals->three_prime_cut = three_prime_cut;
diff --git a/src/trim_paired.c b/src/trim_paired.c
@@ -373,8 +373,8 @@ int paired_main(int argc, char *argv[]) {
             break;
         }
 
-        p1cut = sliding_window(fqrec1, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, trunc_n);
-        p2cut = sliding_window(fqrec2, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, trunc_n);
+        p1cut = sliding_window(fqrec1, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, trunc_n, debug);
+        p2cut = sliding_window(fqrec2, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, trunc_n, debug);
 
         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);
diff --git a/src/trim_single.c b/src/trim_single.c
@@ -191,7 +191,7 @@ int single_main(int argc, char *argv[]) {
 
     while ((l = kseq_read(fqrec)) >= 0) {
 
-        p1cut = sliding_window(fqrec, qualtype, single_length_threshold, single_qual_threshold, no_fiveprime, trunc_n);
+        p1cut = sliding_window(fqrec, qualtype, single_length_threshold, single_qual_threshold, no_fiveprime, trunc_n, debug);
 
         if (debug) printf("P1cut: %d,%d\n", p1cut->five_prime_cut, p1cut->three_prime_cut);
 

Original file line number	Diff line number	Diff line change
`@@ -373,8 +373,8 @@ int paired_main(int argc, char *argv[]) {`
`373`	`373`	`break;`
`374`	`374`	`}`
`375`	`375`
`376`		`- p1cut = sliding_window(fqrec1, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, trunc_n);`
`377`		`- p2cut = sliding_window(fqrec2, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, trunc_n);`
	`376`	`+ p1cut = sliding_window(fqrec1, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, trunc_n, debug);`
	`377`	`+ p2cut = sliding_window(fqrec2, qualtype, paired_length_threshold, paired_qual_threshold, no_fiveprime, trunc_n, debug);`
`378`	`378`
`379`	`379`	`if (debug) printf("p1cut: %d,%d\n", p1cut->five_prime_cut, p1cut->three_prime_cut);`
`380`	`380`	`if (debug) printf("p2cut: %d,%d\n", p2cut->five_prime_cut, p2cut->three_prime_cut);`