Added error detection for cases such as a read's two mates having different names...

[rsem.git] / scanForPairedEndReads.cpp

diff --git a/scanForPairedEndReads.cpp b/scanForPairedEndReads.cpp
index 4d1fbfb79ea16d52a0e1b903bd0e40cfc11f66cd..3669f5faad46972a6b2949a2226a6a159f284b67 100644 (file)
--- a/scanForPairedEndReads.cpp
+++ b/scanForPairedEndReads.cpp
@@ -10,6 +10,7 @@
 #include "sam/bam.h"
 #include "sam/sam.h"
 
+#include "utils.h"
 #include "my_assert.h"
 
 using namespace std;
@@ -28,15 +29,23 @@ inline void add_to_appropriate_arr(bam1_t *b) {
        else arr_partial_unknown.push_back(bam_dup1(b));
 }
 
+char get_pattern_code(uint32_t flag) {
+  if (flag & 0x0040) return (flag & 0x0010 ? 1 : 0);
+  else return (flag & 0x0010 ? 0 : 1);
+}
+
 bool less_than(bam1_t *a, bam1_t *b) {
        int32_t ap1 = min(a->core.pos, a->core.mpos);
        int32_t ap2 = max(a->core.pos, a->core.mpos);
        int32_t bp1 = min(b->core.pos, b->core.mpos);
        int32_t bp2 = max(b->core.pos, b->core.mpos);
+       char apat = get_pattern_code(a->core.flag); // apt: a's pattern of strand and mate information
+       char bpat = get_pattern_code(b->core.flag);
 
        if (a->core.tid != b->core.tid) return a->core.tid < b->core.tid;
        if (ap1 != bp1) return ap1 < bp1;
-       return ap2 < bp2;
+       if (ap2 != bp2) return ap2 < bp2;
+       return apat < bpat;
 }
 
 int main(int argc, char* argv[]) {
@@ -55,7 +64,7 @@ int main(int argc, char* argv[]) {
        string qname;
        bool go_on = (samread(in, b) >= 0);
        bool isPaired;
-       long cnt = 0;
+       HIT_INT_TYPE cnt = 0;
 
        printf("."); fflush(stdout);
 
@@ -66,12 +75,12 @@ int main(int argc, char* argv[]) {
                if (isPaired) {
                        add_to_appropriate_arr(b);
                        while ((go_on = (samread(in, b) >= 0)) && (qname == bam1_qname(b))) {
-                               general_assert(b->core.flag & 0x0001, "Read " + qname + " is detected as both single-end and paired-end read!", true);
+                               general_assert_1(b->core.flag & 0x0001, "Read " + qname + " is detected as both single-end and paired-end read!");
                                add_to_appropriate_arr(b);
                        }
 
-                       general_assert(arr_both.size() % 2 == 0, "Number of first and second mates in read " + qname + "'s full alignments (both mates are aligned) are not matched!", true);
-                       general_assert((arr_partial_1.size() + arr_partial_2.size() + arr_partial_unknown.size()) % 2 == 0, "Number of first and second mates in read " + qname + "'s partial alignments (at most one mate is aligned) are not matched!", true);
+                       general_assert_1(arr_both.size() % 2 == 0, "Number of first and second mates in read " + qname + "'s full alignments (both mates are aligned) are not matched!");
+                       general_assert_1((arr_partial_1.size() + arr_partial_2.size() + arr_partial_unknown.size()) % 2 == 0, "Number of first and second mates in read " + qname + "'s partial alignments (at most one mate is aligned) are not matched!");
 
                        if (!arr_both.empty()) {
                                sort(arr_both.begin(), arr_both.end(), less_than);