Fix thread count crash (#12)

* hange reading procedure with respect to multithreading * Refactor code to use new reading procedure * Fix array out of bound bug * Get max number of threads from input * Set max threads count 1024 * Fix non-overlaping kmer case * Fix thread count crash. * Delete max thread count * Change best distance to INF Co-authored-by: Reza Soltani <reza@host21-107.vpn.ubc.ca> Co-authored-by: Reza Soltani <reza@host17-107.vpn.ubc.ca> Co-authored-by: Reza Soltani <reza@host18-107.vpn.ubc.ca> Co-authored-by: Reza Soltani <reza@host7-107.vpn.ubc.ca> Co-authored-by: Reza Soltani <reza@host24-107.vpn.ubc.ca>
vpc-ccg · Oct 16, 2020 · 1de79ac · 1de79ac
1 parent 5e6086c
commit 1de79ac
Show file tree

Hide file tree

Showing 6 changed files with 40 additions and 12 deletions.
diff --git a/src/circminer.cpp b/src/circminer.cpp
@@ -24,7 +24,7 @@ using namespace std;
 
 char versionNumberMajor[10] = "0";
 char versionNumberMinor[10] = "4";
-char versionNumberPatch[10] = "4";
+char versionNumberPatch[10] = "5";
 
 pthread_mutex_t write_lock;
 pthread_mutex_t pmap_lock;
@@ -56,6 +56,10 @@ int main(int argc, char **argv) {
 
     genome_packer.init(referenceFilename);
 
+    fq_parser1.init();
+    if (pairedEnd)
+        fq_parser2.init();
+
     /****************************************************
      * INDEXING
      ***************************************************/

diff --git a/src/fastq_parser.cpp b/src/fastq_parser.cpp
@@ -5,11 +5,9 @@
 #include "fastq_parser.h"
 
 FASTQParser::FASTQParser(void) {
-    init();
 }
 
 FASTQParser::FASTQParser(char *filename) {
-    init();
     reset(filename);
 }
 

diff --git a/src/filter.cpp b/src/filter.cpp
@@ -276,7 +276,7 @@ int FilterRead::process_mates(int thid, const chain_list &forward_chain, const R
 // 		uint32_t reverse_start = mate_pairs[i].reverse.frags[0].rpos;
 // 		uint32_t reverse_end   = mate_pairs[i].reverse.frags[mate_pairs[i].reverse.chain_len-1].rpos + mate_pairs[i].reverse.frags[mate_pairs[i].reverse.chain_len-1].len - 1;
 
-        bool is_forward_left = is_left_chain(mate_pairs[i].forward, mate_pairs[i].reverse);
+        bool is_forward_left = is_left_chain(mate_pairs[i].forward, mate_pairs[i].reverse, forward_rec->seq_len);
 
 // 		if (forward_start <= reverse_end) {
         if (is_forward_left) {

diff --git a/src/process_circ.cpp b/src/process_circ.cpp
@@ -134,7 +134,9 @@ void ProcessCirc::sort_fq_internal(char *fqname) {
     sprintf(sorted_file, "%s.srt", fqname);
     sprintf(wmode, "%c", 'w');
 
-    FASTQParser fq_parser(fqname);
+    FASTQParser fq_parser;
+    fq_parser.init();
+    fq_parser.reset(fqname);
     vector <RecordStr> all_records;
     Record *current_record;
 
@@ -259,10 +261,13 @@ void ProcessCirc::do_process(void) {
 
     bool is_pe = pairedEnd;
 
-    FASTQParser fq_parser1(fq_file1);
+    FASTQParser fq_parser1;
+    fq_parser1.init();
+    fq_parser1.reset(fq_file1);
     Record *current_record1 = NULL;
 
     FASTQParser fq_parser2;
+    fq_parser2.init();
     Record *current_record2 = NULL;
     if (is_pe) {
         fq_parser2.reset(fq_file2);

diff --git a/src/utils.cpp b/src/utils.cpp
@@ -824,7 +824,7 @@ void reverse_str(char *s, int n, char *revs) {
 }
 
 // is a on the left side?
-bool is_left_chain(chain_t a, chain_t b) {
+bool is_left_chain(chain_t a, chain_t b, int read_length) {
     uint32_t a_beg = a.frags[0].rpos;
     uint32_t b_beg = b.frags[0].rpos;
     uint32_t a_end = a.frags[a.chain_len - 1].rpos + a.frags[a.chain_len - 1].len - 1;
@@ -838,27 +838,48 @@ bool is_left_chain(chain_t a, chain_t b) {
         return a_beg < b_beg;
     } else {
         uint32_t i = 0, j = 0;
+        int best_distance = INF;
+        int best_i = -1;
+        int best_j = -1;
         while (i < a.chain_len and j < b.chain_len) {
             uint32_t bj_beg = b.frags[j].rpos;
             uint32_t ai_end = a.frags[i].rpos + a.frags[i].len - 1;
             if (ai_end < bj_beg) {
+                int distance = bj_beg - ai_end;
+                if (distance < best_distance) {
+                    best_distance = distance;
+                    best_i = i;
+                    best_j = j;
+                }
                 ++i;
                 continue;
             }
             uint32_t ai_beg = a.frags[i].rpos;
             uint32_t bj_end = b.frags[j].rpos + b.frags[j].len - 1;
             if (bj_end < ai_beg) {
+                int distance = ai_beg - bj_end;
+                if (distance < best_distance) {
+                    best_distance = distance;
+                    best_i = i;
+                    best_j = j;
+                }
                 ++j;
                 continue;
             }
+            best_i = i;
+            best_j = j;
+            break;
+        }
 
-            uint32_t common_bp = maxM(ai_beg, bj_beg);
-            int32_t a_ov_qpos = a.frags[i].qpos + (common_bp - a.frags[i].rpos);
-            int32_t b_ov_qpos = b.frags[j].qpos + (common_bp - b.frags[j].rpos);
+
+        uint32_t common_bp = maxM(a.frags[best_i].rpos, b.frags[best_j].rpos);
+        int32_t a_ov_qpos = a.frags[best_i].qpos + (common_bp - a.frags[best_i].rpos);
+        int32_t b_ov_qpos = b.frags[best_j].qpos + (common_bp - b.frags[best_j].rpos);
 
 // 			fprintf(stderr, "OV -> Decision made\n");
+        if (a_ov_qpos < read_length and b_ov_qpos < read_length)
             return a_ov_qpos >= b_ov_qpos;
-        }
+
 // 		fprintf(stderr, "OV -> Ambiguous\n");
 
         return a_beg < b_beg;

diff --git a/src/utils.h b/src/utils.h
@@ -56,7 +56,7 @@ string get_consensus(const vector <string> &vseq);
 
 void reverse_str(char *s, int n, char *revs);
 
-bool is_left_chain(chain_t a, chain_t b);
+bool is_left_chain(chain_t a, chain_t b, int read_length);
 
 void remove_side_introns(MatchedMate &mm, int rlen);