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

[rsem.git] / BamConverter.h

diff --git a/BamConverter.h b/BamConverter.h
index 4ea7d3e673e3942e4de7a678acb5be15faf03c11..d2b8b7324df6d61b59e3f03d92dbe0a4941a3bff 100644 (file)
--- a/BamConverter.h
+++ b/BamConverter.h
@@ -39,7 +39,7 @@ private:
        void writeCollapsedLines();
        void flipSeq(uint8_t*, int);
        void flipQual(uint8_t*, int);
-       void addXSTag(bam1_t*, const Transcript&);
+       void modifyTags(bam1_t*, const Transcript&); // modify MD tag and XS tag if needed
 };
 
 BamConverter::BamConverter(const char* inpF, const char* outF, const char* chr_list, Transcripts& transcripts)
@@ -119,14 +119,13 @@ void BamConverter::process() {
                // at least one segment is not properly mapped
                bool notgood = (b->core.flag & 0x0004) || (isPaired && (b2->core.flag & 0x0004));
                
-
-               if (isPaired && notgood) assert((b->core.flag & 0x0004) && (b2->core.flag & 0x0004));
-
+               if (isPaired && notgood) general_assert((b->core.flag & 0x0004) && (b2->core.flag & 0x0004), cstrtos(bam1_qname(b)) + "'s two mates are not all marked as unalignable!");
 
                if (!notgood) {
                        // for collapsing
                        if (isPaired) {
                                assert(b->core.tid == b2->core.tid);
+                               general_assert(b->core.tid == b2->core.tid, cstrtos(bam1_qname(b)) + "'s two mates are aligned to two different transcripts!");
                                if ((b->core.flag & 0x0080) && (b2->core.flag & 0x0040)) {
                                        bam1_t *tmp = b; b = b2; b2 = tmp;
                                }
@@ -211,7 +210,7 @@ void BamConverter::convert(bam1_t* b, const Transcript& transcript) {
        b->core.n_cigar = core_n_cigar;
        b->core.bin = bam_reg2bin(b->core.pos, bam_calend(&(b->core), bam1_cigar(b)));
 
-       addXSTag(b, transcript); // check if need to add XS tag, if need, add it
+       modifyTags(b, transcript); // check if need to add XS tag, if need, add it
 }
 
 inline void BamConverter::writeCollapsedLines() {
@@ -227,7 +226,8 @@ inline void BamConverter::writeCollapsedLines() {
                                memcpy(bam_aux_get(tmp_b, "ZW") + 1, (uint8_t*)&(prb), bam_aux_type2size('f'));
                                tmp_b->core.qual = getMAPQ(prb);
                        }
-                       else tmp_b->core.qual = getMAPQ(1.0);
+                       // otherwise, just use the MAPQ score of the orignal alignment
+
                        samwrite(out, tmp_b);
                        if (isPaired) {
                                if (p != NULL) memcpy(bam_aux_get(tmp_b2, "ZW") + 1, (uint8_t*)&(prb), bam_aux_type2size('f'));
@@ -271,14 +271,58 @@ inline void BamConverter::flipQual(uint8_t* q, int readlen) {
        }
 }
 
-inline void BamConverter::addXSTag(bam1_t* b, const Transcript& transcript) {
-       uint32_t* p = bam1_cigar(b);
+inline void BamConverter::modifyTags(bam1_t* b, const Transcript& transcript) {
+       char strand = transcript.getStrand();
+       uint8_t *s = NULL;
+
+       if (strand == '-') {
+         s = bam_aux_get(b, "MD");
+         if ((s != NULL) && (*(s) == 'Z') && (bam_aux2Z(s) != NULL)) {
+           char *mis = bam_aux2Z(s);
+           int len = strlen(mis);
+           char *tmp = new char[len];
+           int cur_type = -1, fr = -1, type, base;
+           for (int i = 0; i < len; i++) {
+             type = (mis[i] >= '0' && mis[i] <= '9');
+             if (cur_type != type) {
+               switch(cur_type) {
+               case 0:
+                 base = len - 1;
+                 if (mis[fr] == '^') { tmp[len - i] = mis[fr]; ++fr; ++base; }
+                 for (int j = fr; j < i; j++) tmp[base - j] = ((mis[j] == 'A' || mis[j] == 'C' || mis[j] == 'G' || mis[j] == 'T') ? getOpp(mis[j]) : mis[j]);
+                 break;
+               case 1: 
+                 base = len - i - fr;
+                 for (int j = fr; j < i; j++) tmp[base + j] = mis[j]; 
+                 break; 
+               }
+               cur_type = type;
+               fr = i;
+             }
+           }
+           switch(cur_type) {
+           case 0:
+             base = len - 1;
+             if (mis[fr] == '^') { tmp[0] = mis[fr]; ++fr; ++base; }
+             for (int j = fr; j < len; j++) tmp[base - j] = ((mis[j] == 'A' || mis[j] == 'C' || mis[j] == 'G' || mis[j] == 'T') ? getOpp(mis[j]) : mis[j]);
+             break;
+           case 1: 
+             for (int j = fr; j < len; j++) tmp[j - fr] = mis[j]; 
+             break; 
+           }
+           strncpy(mis, tmp, len);
+           delete[] tmp;
+         }
+       }
+
+       // append XS:A field if necessary
+       s = bam_aux_get(b, "XS");
+       if (s != NULL) bam_aux_del(b, s);
        bool hasN = false;
+       uint32_t* p = bam1_cigar(b);
        for (int i = 0; i < (int)b->core.n_cigar; i++)
                if ((*(p + i) & BAM_CIGAR_MASK) == BAM_CREF_SKIP) { hasN = true; break; }
-       if (!hasN) return;
-       char strand = transcript.getStrand();
-       bam_aux_append(b, "XS", 'A', 1, (uint8_t*)&strand);
+       if (hasN) bam_aux_append(b, "XS", 'A', 1, (uint8_t*)&strand);
 }
 
 #endif /* BAMCONVERTER_H_ */