BamConverter.h

   1 #ifndef BAMCONVERTER_H_
   2 #define BAMCONVERTER_H_
   3
   4 #include<cstdio>
   5 #include<cstring>
   6 #include<cassert>
   7 #include<string>
   8 #include<map>
   9
  10 #include <stdint.h>
  11 #include "sam/bam.h"
  12 #include "sam/sam.h"
  13 #include "sam_rsem_aux.h"
  14 #include "sam_rsem_cvt.h"
  15
  16 #include "utils.h"
  17 #include "my_assert.h"
  18 #include "bc_aux.h"
  19 #include "Transcript.h"
  20 #include "Transcripts.h"
  21
  22 class BamConverter {
  23 public:
  24         BamConverter(const char*, const char*, const char*, Transcripts&);
  25         ~BamConverter();
  26
  27         void process();
  28 private:
  29         samfile_t *in, *out;
  30         Transcripts& transcripts;
  31
  32         std::map<std::string, int> refmap;
  33         std::map<std::string, int>::iterator iter;
  34
  35         CollapseMap collapseMap;
  36
  37         void convert(bam1_t*, const Transcript&);
  38
  39         void writeCollapsedLines();
  40         void flipSeq(uint8_t*, int);
  41         void flipQual(uint8_t*, int);
  42         void modifyTags(bam1_t*, const Transcript&); // modify MD tag and XS tag if needed
  43 };
  44
  45 BamConverter::BamConverter(const char* inpF, const char* outF, const char* chr_list, Transcripts& transcripts)
  46         : transcripts(transcripts)
  47 {
  48         general_assert(transcripts.getType() == 0, "Genome information is not provided! RSEM cannot convert the transcript bam file!");
  49
  50         in = samopen(inpF, "rb", NULL);
  51         assert(in != 0);
  52
  53         transcripts.buildMappings(in->header->n_targets, in->header->target_name);
  54
  55         bam_header_t *out_header = sam_header_read2(chr_list);
  56
  57         refmap.clear();
  58         for (int i = 0; i < out_header->n_targets; i++) {
  59                 refmap[out_header->target_name[i]] = i;
  60         }
  61
  62         if (in->header->l_text > 0) {
  63                 char comment[] = "@CO\tThis BAM file is processed by rsem-tbam2gam to convert from transcript coordinates into genomic coordinates.\n";
  64                 int comment_len = strlen(comment);
  65
  66                 //Filter @SQ fields if the BAM file is user provided
  67                 char *text = in->header->text;
  68                 int l_text = in->header->l_text;
  69                 char *new_text = new char[l_text + comment_len];
  70                 int pos = 0, s = 0;
  71                 while (pos < l_text) {
  72                         if ((pos + 2 < l_text) && (text[pos] == '@') && (text[pos + 1] == 'S') && (text[pos + 2] == 'Q')) {
  73                                 pos += 3;
  74                                 while (pos < l_text && text[pos] != '\n') ++pos;
  75                         }
  76                         else new_text[s++] = text[pos];
  77                         ++pos;
  78                 }
  79                 strncpy(new_text + s, comment, comment_len);
  80                 s += comment_len;
  81
  82                 append_header_text(out_header, new_text, s);
  83                 delete[] new_text;
  84         }
  85
  86         out = samopen(outF, "wb", out_header);
  87         assert(out != 0);
  88
  89         bam_header_destroy(out_header);
  90 }
  91
  92 BamConverter::~BamConverter() {
  93         samclose(in);
  94         samclose(out);
  95 }
  96
  97 void BamConverter::process() {
  98         bam1_t *b, *b2;
  99         std::string cqname;
 100         bool isPaired = false;
 101
 102         HIT_INT_TYPE cnt = 0;
 103
 104         cqname = "";
 105         b = bam_init1(); b2 = bam_init1();
 106
 107         while (samread(in, b) >= 0) {
 108                 ++cnt;
 109                 isPaired = (b->core.flag & 0x0001) > 0;
 110                 if (isPaired) {
 111                         assert(samread(in, b2) >= 0 && (b2->core.flag & 0x0001));
 112                         assert((b->core.flag & 0x0001) && (b2->core.flag & 0x0001));
 113                         assert(((b->core.flag & 0x0040) && (b2->core.flag & 0x0080)) || ((b->core.flag & 0x0080) && (b2->core.flag & 0x0040)));
 114                         ++cnt;
 115                 }
 116
 117                 if (cnt % 1000000 == 0) { printf("."); fflush(stdout); }
 118
 119                 // at least one segment is not properly mapped
 120                 bool notgood = (b->core.flag & 0x0004) || (isPaired && (b2->core.flag & 0x0004));
 121
 122                 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!");
 123
 124                 if (!notgood) {
 125                         // for collapsing
 126                         if (isPaired) {
 127                                 assert(b->core.tid == b2->core.tid);
 128                                 general_assert(b->core.tid == b2->core.tid, cstrtos(bam1_qname(b)) + "'s two mates are aligned to two different transcripts!");
 129                                 if ((b->core.flag & 0x0080) && (b2->core.flag & 0x0040)) {
 130                                         bam1_t *tmp = b; b = b2; b2 = tmp;
 131                                 }
 132                         }
 133
 134                         const Transcript& transcript = transcripts.getTranscriptViaEid(b->core.tid + 1);
 135
 136                         convert(b, transcript);
 137                         if (isPaired) {
 138                                 convert(b2, transcript);
 139                                 b->core.mpos = b2->core.pos;
 140                                 b2->core.mpos = b->core.pos;
 141                         }
 142
 143                         if (cqname != bam1_qname(b)) {
 144                                 writeCollapsedLines();
 145                                 cqname = bam1_qname(b);
 146                                 collapseMap.init(isPaired);
 147                         }
 148
 149                         uint8_t *p = bam_aux_get(b, "ZW");
 150                         float prb = (p != NULL? bam_aux2f(p) : 1.0);
 151                         collapseMap.insert(b, b2, prb);
 152                 }
 153                 else {
 154                         assert(cqname != bam1_qname(b));
 155
 156                         writeCollapsedLines();
 157                         cqname = bam1_qname(b);
 158                         collapseMap.init(isPaired);
 159
 160                         samwrite(out, b);
 161                         if (isPaired) samwrite(out, b2);
 162                 }
 163         }
 164
 165         writeCollapsedLines();
 166
 167         bam_destroy1(b);
 168         bam_destroy1(b2);
 169
 170         if (cnt >= 1000000) printf("\n");
 171 }
 172
 173 void BamConverter::convert(bam1_t* b, const Transcript& transcript) {
 174         int pos = b->core.pos;
 175         int readlen = b->core.l_qseq;
 176
 177         general_assert(readlen > 0, "One alignment line has SEQ field as *. RSEM does not support this currently!");
 178
 179         iter = refmap.find(transcript.getSeqName());
 180         assert(iter != refmap.end());
 181         b->core.tid = iter->second;
 182         if (b->core.flag & 0x0001) { b->core.mtid = b->core.tid; }
 183         b->core.qual = 255; // set to not available temporarily
 184
 185         if (transcript.getStrand() == '-') {
 186                 b->core.flag ^= 0x0010;
 187                 if (b->core.flag & 0x0001) {
 188                         b->core.flag ^= 0x0020;
 189                         b->core.isize = -b->core.isize;
 190                 }
 191                 flipSeq(bam1_seq(b), readlen);
 192                 flipQual(bam1_qual(b), readlen);
 193         }
 194
 195         std::vector<uint32_t> data;
 196         data.clear();
 197
 198         int core_pos, core_n_cigar;
 199         tr2chr(transcript, pos + 1, pos + readlen, core_pos, core_n_cigar, data);
 200         assert(core_pos >= 0);
 201
 202         int rest_len = b->data_len - b->core.l_qname - b->core.n_cigar * 4;
 203         b->data_len = b->core.l_qname + core_n_cigar * 4 + rest_len;
 204         expand_data_size(b);
 205         uint8_t* pt = b->data + b->core.l_qname;
 206         memmove(pt + core_n_cigar * 4, pt + b->core.n_cigar * 4, rest_len);
 207         for (int i = 0; i < core_n_cigar; i++) { memmove(pt, &data[i], 4); pt += 4; }
 208
 209         b->core.pos = core_pos;
 210         b->core.n_cigar = core_n_cigar;
 211         b->core.bin = bam_reg2bin(b->core.pos, bam_calend(&(b->core), bam1_cigar(b)));
 212
 213         modifyTags(b, transcript); // check if need to add XS tag, if need, add it
 214 }
 215
 216 inline void BamConverter::writeCollapsedLines() {
 217         bam1_t *tmp_b = NULL,*tmp_b2 = NULL;
 218         float prb;
 219         bool isPaired;
 220         uint8_t *p;
 221
 222         if (!collapseMap.empty(isPaired)) {
 223                 while (collapseMap.next(tmp_b, tmp_b2, prb)) {
 224                         p = bam_aux_get(tmp_b, "ZW");
 225                         if (p != NULL) {
 226                                 memcpy(bam_aux_get(tmp_b, "ZW") + 1, (uint8_t*)&(prb), bam_aux_type2size('f'));
 227                                 tmp_b->core.qual = getMAPQ(prb);
 228                         }
 229                         // otherwise, just use the MAPQ score of the orignal alignment
 230
 231                         samwrite(out, tmp_b);
 232                         if (isPaired) {
 233                                 if (p != NULL) memcpy(bam_aux_get(tmp_b2, "ZW") + 1, (uint8_t*)&(prb), bam_aux_type2size('f'));
 234                                 tmp_b2->core.qual = tmp_b->core.qual;
 235                                 samwrite(out, tmp_b2);
 236                         }
 237                         bam_destroy1(tmp_b);
 238                         if (isPaired) bam_destroy1(tmp_b2);
 239                 }
 240         }
 241 }
 242
 243 inline void BamConverter::flipSeq(uint8_t* s, int readlen) {
 244         uint8_t code, base;
 245         std::vector<uint8_t> seq;
 246
 247         code = 0; base = 0;
 248         seq.clear();
 249         for (int i = 0; i < readlen; i++) {
 250                 switch (bam1_seqi(s, readlen - i - 1)) {
 251                 case 1: base = 8; break;
 252                 case 2: base = 4; break;
 253                 case 4: base = 2; break;
 254                 case 8: base = 1; break;
 255                 case 15: base = 15; break;
 256                 default: assert(false);
 257                 }
 258                 code |=  base << (4 * (1 - i % 2));
 259                 if (i % 2 == 1) { seq.push_back(code); code = 0; }
 260         }
 261         if (readlen % 2 == 1) { seq.push_back(code); }
 262
 263         for (int i = 0; i < (int)seq.size(); i++) s[i] = seq[i];
 264 }
 265
 266 inline void BamConverter::flipQual(uint8_t* q, int readlen) {
 267         int32_t mid = readlen / 2;
 268         uint8_t tmp;
 269         for (int i = 0; i < mid; i++) {
 270                 tmp = q[i]; q[i] = q[readlen - i - 1]; q[readlen -i -1] = tmp;
 271         }
 272 }
 273
 274 inline void BamConverter::modifyTags(bam1_t* b, const Transcript& transcript) {
 275         char strand = transcript.getStrand();
 276         uint8_t *s = NULL;
 277
 278         if (strand == '-') {
 279           s = bam_aux_get(b, "MD");
 280           if ((s != NULL) && (*(s) == 'Z') && (bam_aux2Z(s) != NULL)) {
 281             char *mis = bam_aux2Z(s);
 282             int len = strlen(mis);
 283             char *tmp = new char[len];
 284             int cur_type = -1, fr = -1, type, base;
 285             for (int i = 0; i < len; i++) {
 286               type = (mis[i] >= '0' && mis[i] <= '9');
 287               if (cur_type != type) {
 288                 switch(cur_type) {
 289                 case 0:
 290                   base = len - 1;
 291                   if (mis[fr] == '^') { tmp[len - i] = mis[fr]; ++fr; ++base; }
 292                   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]);
 293                   break;
 294                 case 1:
 295                   base = len - i - fr;
 296                   for (int j = fr; j < i; j++) tmp[base + j] = mis[j];
 297                   break;
 298                 }
 299                 cur_type = type;
 300                 fr = i;
 301               }
 302             }
 303             switch(cur_type) {
 304             case 0:
 305               base = len - 1;
 306               if (mis[fr] == '^') { tmp[0] = mis[fr]; ++fr; ++base; }
 307               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]);
 308               break;
 309             case 1:
 310               for (int j = fr; j < len; j++) tmp[j - fr] = mis[j];
 311               break;
 312             }
 313             strncpy(mis, tmp, len);
 314             delete[] tmp;
 315           }
 316         }
 317
 318         // append XS:A field if necessary
 319         s = bam_aux_get(b, "XS");
 320         if (s != NULL) bam_aux_del(b, s);
 321         bool hasN = false;
 322         uint32_t* p = bam1_cigar(b);
 323         for (int i = 0; i < (int)b->core.n_cigar; i++)
 324                 if ((*(p + i) & BAM_CIGAR_MASK) == BAM_CREF_SKIP) { hasN = true; break; }
 325         if (hasN) bam_aux_append(b, "XS", 'A', 1, (uint8_t*)&strand);
 326 }
 327
 328 #endif /* BAMCONVERTER_H_ */