Merge pull request #25 from smithlabcode/more_independence_between_ends

More independence between ends

src/abismal.cpp

@@ -1596,15 +1596,19 @@ map_fragments(const uint32_t max_candidates,
   res1.reset(read1.size());
   res2.reset(read2.size());
 
-  if (!read1.empty() && !read2.empty()) {
+  if (read1.empty() && read2.empty()) return false;
+
+  if (!read1.empty()) {
     prep_read<cmp>(read1, pread1);
     pack_read(pread1, packed_pread);
     process_seeds<strand_code1>(max_candidates,
         counter_st, counter_three_st,
         index_st, index_three_st,
         genome_st, pread1, packed_pread, res1
     );
+  }
 
+  if (!read2.empty()) {
     const string read_rc(revcomp(read2));
     prep_read<cmp>(read_rc, pread2);
     pack_read(pread2, packed_pread);
@@ -1613,11 +1617,10 @@ map_fragments(const uint32_t max_candidates,
         index_st, index_three_st,
       genome_st, pread2, packed_pread, res2
     );
+  }
 
     return select_maps<swap_ends>(pread1, pread2, cigar1, cigar2,
       res1, res2, mem_scr1, res_se1, res_se2, aln, best);
-  }
-  return false;
 }
 
 template <const conversion_type conv>
@@ -1721,8 +1724,8 @@ map_paired_ended(const bool VERBOSE,
       bests_se1[i].reset(readlen1);
       bests_se2[i].reset(readlen2);
 
-      if (
-        !map_fragments<conv, false,
+      const bool strand_pm_success =
+        map_fragments<conv, false,
                      get_strand_code('+',conv),
                      get_strand_code('-', flip_conv(conv))>(
           max_candidates, reads1[i], reads2[i],
@@ -1733,9 +1736,10 @@ map_paired_ended(const bool VERBOSE,
           genome_st, pread1, pread2_rc, packed_pread,
           cigar1[i], cigar2[i], aln, res1, res2, mem_scr1,
           res_se1, res_se2, bests[i]
-        ) ||
+        );
 
-        !map_fragments<!conv, true,
+      const bool strand_mp_success =
+        map_fragments<!conv, true,
                      get_strand_code('+', flip_conv(conv)),
                      get_strand_code('-', conv)>(
           max_candidates, reads2[i], reads1[i],
@@ -1746,7 +1750,8 @@ map_paired_ended(const bool VERBOSE,
           genome_st, pread2, pread1_rc, packed_pread,
           cigar2[i], cigar1[i], aln, res2, res1, mem_scr1,
           res_se2, res_se1, bests[i]
-        )) {
+        );
+      if (!strand_pm_success && !strand_mp_success) {
         bests[i].reset();
         res_se1.reset();
         res_se2.reset();
@@ -1881,54 +1886,55 @@ map_paired_ended_rand(const bool VERBOSE, const bool allow_ambig,
       bests_se1[i].reset(readlen1);
       bests_se2[i].reset(readlen2);
 
-      if (
       // GS: (1) T/A-rich +/- strand
-      !map_fragments<t_rich, false,
+      const bool richness_ta_strand_pm_success =
+       map_fragments<t_rich, false,
                    get_strand_code('+', t_rich),
                    get_strand_code('-', a_rich)>(
         max_candidates, reads1[i], reads2[i],
         counter_st, counter_t_st,
         index_st, index_t_st,
         genome_st, pread1_t, pread2_t_rc, packed_pread,
         cigar1[i], cigar2[i], aln, res1, res2, mem_scr1,
-        res_se1, res_se2, bests[i]
-        ) ||
-
+        res_se1, res_se2, bests[i]);
       // GS: (2) T/A-rich, -/+ strand
-      !map_fragments<a_rich, true,
+      const bool richness_ta_strand_mp_success =
+       map_fragments<a_rich, true,
                    get_strand_code('+', a_rich),
                    get_strand_code('-', t_rich)>(
         max_candidates, reads2[i], reads1[i],
         counter_st, counter_a_st,
         index_st, index_a_st,
         genome_st, pread2_a, pread1_a_rc, packed_pread,
         cigar2[i], cigar1[i], aln, res2, res1, mem_scr1,
-        res_se2, res_se1, bests[i]
-      ) ||
-
+        res_se2, res_se1, bests[i]);
       // GS: (3) A/T-rich +/- strand
-      !map_fragments<a_rich, false,
+      const bool richness_at_strand_pm_success =
+       map_fragments<a_rich, false,
                    get_strand_code('+', a_rich),
                    get_strand_code('-', t_rich)>(
         max_candidates, reads1[i], reads2[i],
         counter_st, counter_a_st,
         index_st, index_a_st,
         genome_st, pread1_a, pread2_a_rc, packed_pread,
         cigar1[i], cigar2[i], aln, res1, res2, mem_scr1,
-        res_se1, res_se2, bests[i]
-      ) ||
-
+        res_se1, res_se2, bests[i]);
       // GS: (4) A/T-rich, -/+ strand
-      !map_fragments<t_rich, true,
+      const bool richness_at_strand_mp_success =
+       map_fragments<t_rich, true,
                    get_strand_code('+', t_rich),
                    get_strand_code('-', a_rich)>(
         max_candidates, reads2[i], reads1[i],
         counter_st, counter_t_st,
         index_st, index_t_st,
         genome_st, pread2_t, pread1_t_rc, packed_pread,
         cigar2[i], cigar1[i], aln, res2, res1, mem_scr1,
-        res_se2, res_se1, bests[i]
-      )) {
+        res_se2, res_se1, bests[i]);
+
+      if (!richness_ta_strand_pm_success &&
+          !richness_ta_strand_mp_success &&
+          !richness_at_strand_pm_success &&
+          !richness_at_strand_mp_success) {
         bests[i].reset();
         res_se1.reset();
         res_se2.reset();