Fixed a bug to make RSEM generated BAM file valid input for RSEM

[rsem.git] / BamWriter.h

diff --git a/BamWriter.h b/BamWriter.h
index 782f5bdeb6705748a8919265e854b54d1f418133..8c568f348f361fd009c964cd2b9ad021e50fc257 100644 (file)
--- a/BamWriter.h
+++ b/BamWriter.h
@@ -125,25 +125,27 @@ void BamWriter::work(HitWrapper<PairedEndHit> wrapper) {
                bool notgood = (b->core.flag & 0x0004) || (b2->core.flag & 0x0004);
 
                if (!notgood) {
-                 //swap if b is mate 2
-                 if (b->core.flag & 0x0080) {
-                   assert(b2->core.flag & 0x0040);
-                   bam1_t *tmp = b;
-                   b = b2; b2 = tmp;
-                 }
+                       //swap if b is mate 2
+                       if (b->core.flag & 0x0080) {
+                               assert(b2->core.flag & 0x0040);
+                               bam1_t *tmp = b;
+                               b = b2; b2 = tmp;
+                       }
 
-                 hit = wrapper.getNextHit();
-                 assert(hit != NULL);
+                       hit = wrapper.getNextHit();
+                       assert(hit != NULL);
 
-                 assert(transcripts.getInternalSid(b->core.tid + 1) == hit->getSid());
-                 assert(transcripts.getInternalSid(b2->core.tid + 1) == hit->getSid());
+                       assert(transcripts.getInternalSid(b->core.tid + 1) == hit->getSid());
+                       assert(transcripts.getInternalSid(b2->core.tid + 1) == hit->getSid());
 
-                 convert(b, hit->getConPrb());
-                 convert(b2, hit->getConPrb());
+                       convert(b, hit->getConPrb());
+                       convert(b2, hit->getConPrb());
 
-                 b->core.mpos = b2->core.pos;
-                 b2->core.mpos = b->core.pos;
+                       b->core.mpos = b2->core.pos;
+                       b2->core.mpos = b->core.pos;
                }
+
+               /*
                else {
                  // if only one mate can be aligned, mask it as unaligned and put an additional tag Z0:A:!
                  char exclamation = '!';
@@ -156,7 +158,7 @@ void BamWriter::work(HitWrapper<PairedEndHit> wrapper) {
                    bam_aux_append(b2, "Z0", 'A', bam_aux_type2size('A'), (uint8_t*)&exclamation);
                  }
                }
-
+               */
 
                samwrite(out, b);
                samwrite(out, b2);
@@ -171,33 +173,7 @@ void BamWriter::work(HitWrapper<PairedEndHit> wrapper) {
 }
 
 void BamWriter::convert(bam1_t *b, double prb) {
-       const Transcript& transcript = transcripts.getTranscriptViaEid(b->core.tid + 1);
-
-       int pos = b->core.pos;
-       int readlen = b->core.l_qseq;
-
-       std::vector<uint32_t> data;
-       data.clear();
-
-       int core_pos, core_n_cigar;
-       std::vector<Interval> vec;
-       vec.assign(1, Interval(1, transcript.getLength()));
-       // make an artificial chromosome coordinates for the transcript to get new CIGAR strings
-       tr2chr(Transcript("", "", "", '+', vec, ""), pos + 1, pos + readlen, core_pos, core_n_cigar, data);
-       assert(core_pos >= 0);
-
-       int rest_len = b->data_len - b->core.l_qname - b->core.n_cigar * 4;
-       b->data_len = b->core.l_qname + core_n_cigar * 4 + rest_len;
-       expand_data_size(b);
-       uint8_t* pt = b->data + b->core.l_qname;
-       memmove(pt + core_n_cigar * 4, pt + b->core.n_cigar * 4, rest_len);
-       for (int i = 0; i < core_n_cigar; i++) { memmove(pt, &data[i], 4); pt += 4; }
-
-       b->core.pos = core_pos;
-       b->core.n_cigar = core_n_cigar;
        b->core.qual = getMAPQ(prb);
-       b->core.bin = bam_reg2bin(b->core.pos, bam_calend(&(b->core), bam1_cigar(b)));
-
        float val = (float)prb;
        bam_aux_append(b, "ZW", 'f', bam_aux_type2size('f'), (uint8_t*)&val);
 }