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 "bc_aux.h"
  18 #include "Transcript.h"
  19 #include "Transcripts.h"
  20
  21 class BamConverter {
  22 public:
  23         BamConverter(const char*, const char*, const char*, Transcripts&);
  24         ~BamConverter();
  25
  26         void process();
  27 private:
  28         samfile_t *in, *out;
  29         Transcripts& transcripts;
  30
  31         std::map<std::string, int> refmap;
  32         std::map<std::string, int>::iterator iter;
  33
  34         CollapseMap collapseMap;
  35
  36         void convert(bam1_t*, const Transcript&);
  37
  38         void writeCollapsedLines();
  39         void flipSeq(uint8_t*, int);
  40         void flipQual(uint8_t*, int);
  41         void addXSTag(bam1_t*, const Transcript&);
  42 };
  43
  44 BamConverter::BamConverter(const char* inpF, const char* outF, const char* chr_list, Transcripts& transcripts)
  45         : transcripts(transcripts)
  46 {
  47         if (transcripts.getType() != 0)
  48                 exitWithError("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         bam_header_t *out_header = sam_header_read2(chr_list);
  54         refmap.clear();
  55         for (int i = 0; i < out_header->n_targets; i++) {
  56                 refmap[out_header->target_name[i]] = i;
  57         }
  58
  59         append_header_text(out_header, in->header->text, in->header->l_text);
  60
  61         out = samopen(outF, "wb", out_header);
  62         assert(out != 0);
  63
  64         bam_header_destroy(out_header);
  65 }
  66
  67 BamConverter::~BamConverter() {
  68         samclose(in);
  69         samclose(out);
  70 }
  71
  72 void BamConverter::process() {
  73         bam1_t *b, *b2;
  74         std::string cqname;
  75         bool isPaired = false;
  76
  77         int cnt = 0;
  78
  79         cqname = "";
  80         b = bam_init1(); b2 = bam_init1();
  81
  82         while (samread(in, b) >= 0) {
  83                 ++cnt;
  84                 isPaired = (b->core.flag & 0x0001) > 0;
  85                 if (isPaired) {
  86                         assert(samread(in, b2) >= 0 && (b2->core.flag & 0x0001) && b->core.tid == b2->core.tid);
  87                         assert((b->core.flag & 0x0040) && (b2->core.flag & 0x0080)); // for collapsing
  88                         ++cnt;
  89                 }
  90
  91                 if (cnt % 1000000 == 0) { printf("."); fflush(stdout); }
  92
  93                 // at least one segment is not properly mapped
  94                 if ((b->core.flag & 0x0004) || isPaired && (b2->core.flag & 0x0004)) continue;
  95
  96                 const Transcript& transcript = transcripts.getTranscriptAt(b->core.tid + 1);
  97
  98                 convert(b, transcript);
  99                 if (isPaired) {
 100                         convert(b2, transcript);
 101                         b->core.mpos = b2->core.pos;
 102                         b2->core.mpos = b->core.pos;
 103                 }
 104
 105                 if (cqname != bam1_qname(b)) {
 106                         writeCollapsedLines();
 107                         cqname = bam1_qname(b);
 108                         collapseMap.init(isPaired);
 109                 }
 110
 111                 collapseMap.insert(b, b2, bam_aux2f(bam_aux_get(b, "ZW")));
 112         }
 113
 114         writeCollapsedLines();
 115
 116         bam_destroy1(b);
 117         bam_destroy1(b2);
 118
 119         if (cnt >= 1000000) printf("\n");
 120 }
 121
 122 void BamConverter::convert(bam1_t* b, const Transcript& transcript) {
 123         int pos = b->core.pos;
 124         int readlen = b->core.l_qseq;
 125
 126         if (readlen == 0) exitWithError("One alignment line has SEQ field as *. RSEM does not support this currently!");
 127
 128         iter = refmap.find(transcript.getSeqName());
 129         assert(iter != refmap.end());
 130         b->core.tid = iter->second;
 131         if (b->core.flag & 0x0001) { b->core.mtid = b->core.tid; }
 132         b->core.qual = 255; // set to not available temporarily
 133
 134         if (transcript.getStrand() == '-') {
 135                 b->core.flag ^= 0x0010;
 136                 if (b->core.flag & 0x0001) {
 137                         b->core.flag ^= 0x0020;
 138                         b->core.isize = -b->core.isize;
 139                 }
 140                 flipSeq(bam1_seq(b), readlen);
 141                 flipQual(bam1_qual(b), readlen);
 142         }
 143
 144         std::vector<uint32_t> data;
 145         data.clear();
 146
 147         int core_pos, core_n_cigar;
 148         tr2chr(transcript, pos + 1, pos + readlen, core_pos, core_n_cigar, data);
 149         assert(core_pos >= 0);
 150
 151         int rest_len = b->data_len - b->core.l_qname - b->core.n_cigar * 4;
 152         b->data_len = b->core.l_qname + core_n_cigar * 4 + rest_len;
 153         expand_data_size(b);
 154         uint8_t* pt = b->data + b->core.l_qname;
 155         memmove(pt + core_n_cigar * 4, pt + b->core.n_cigar * 4, rest_len);
 156         for (int i = 0; i < core_n_cigar; i++) { memmove(pt, &data[i], 4); pt += 4; }
 157
 158         b->core.pos = core_pos;
 159         b->core.n_cigar = core_n_cigar;
 160         b->core.bin = bam_reg2bin(b->core.pos, bam_calend(&(b->core), bam1_cigar(b)));
 161
 162         addXSTag(b, transcript); // check if need to add XS tag, if need, add it
 163 }
 164
 165 inline void BamConverter::writeCollapsedLines() {
 166         bam1_t *tmp_b = NULL,*tmp_b2 = NULL;
 167         float prb;
 168         bool isPaired;
 169
 170         if (!collapseMap.empty(isPaired)) {
 171                 while (collapseMap.next(tmp_b, tmp_b2, prb)) {
 172                         memcpy(bam_aux_get(tmp_b, "ZW") + 1, (uint8_t*)&(prb), bam_aux_type2size('f'));
 173                         tmp_b->core.qual = getMAPQ(prb);
 174                         if (tmp_b->core.qual > 0) {
 175                                 samwrite(out, tmp_b);
 176                                 if (isPaired) {
 177                                         memcpy(bam_aux_get(tmp_b2, "ZW") + 1, (uint8_t*)&(prb), bam_aux_type2size('f'));
 178                                         tmp_b2->core.qual = tmp_b->core.qual;
 179                                         samwrite(out, tmp_b2);
 180                                 }
 181                         }
 182                         bam_destroy1(tmp_b);
 183                         if (isPaired) bam_destroy1(tmp_b2);
 184
 185                 }
 186         }
 187 }
 188
 189 inline void BamConverter::flipSeq(uint8_t* s, int readlen) {
 190         uint8_t code, base;
 191         std::vector<uint8_t> seq;
 192
 193         code = 0; base = 0;
 194         seq.clear();
 195         for (int i = 0; i < readlen; i++) {
 196                 switch (bam1_seqi(s, readlen - i - 1)) {
 197                 case 1: base = 8; break;
 198                 case 2: base = 4; break;
 199                 case 4: base = 2; break;
 200                 case 8: base = 1; break;
 201                 case 15: base = 15; break;
 202                 default: assert(false);
 203                 }
 204                 code |=  base << (4 * (1 - i % 2));
 205                 if (i % 2 == 1) { seq.push_back(code); code = 0; }
 206         }
 207         if (readlen % 2 == 1) { seq.push_back(code); }
 208
 209         for (int i = 0; i < (int)seq.size(); i++) s[i] = seq[i];
 210 }
 211
 212 inline void BamConverter::flipQual(uint8_t* q, int readlen) {
 213         int32_t mid = readlen / 2;
 214         uint8_t tmp;
 215         for (int i = 0; i < mid; i++) {
 216                 tmp = q[i]; q[i] = q[readlen - i - 1]; q[readlen -i -1] = tmp;
 217         }
 218 }
 219
 220 inline void BamConverter::addXSTag(bam1_t* b, const Transcript& transcript) {
 221         uint32_t* p = bam1_cigar(b);
 222         bool hasN = false;
 223         for (int i = 0; i < (int)b->core.n_cigar; i++)
 224                 if ((*(p + i) & BAM_CIGAR_MASK) == BAM_CREF_SKIP) { hasN = true; break; }
 225         if (!hasN) return;
 226         char strand = transcript.getStrand();
 227         bam_aux_append(b, "XS", 'A', 1, (uint8_t*)&strand);
 228 }
 229
 230 #endif /* BAMCONVERTER_H_ */