sam_rsem_cvt.h

   1 #ifndef SAM_RSEM_CVT_H_
   2 #define SAM_RSEM_CVT_H_
   3
   4 #include<vector>
   5
   6 #include "stdint.h"
   7 #include "sam/bam.h"
   8
   9 #include "Transcript.h"
  10 #include "Transcripts.h"
  11
  12 uint8_t getMAPQ(double val) {
  13         double err = 1.0 - val;
  14         if (err <= 1e-10) return 100;
  15         return (uint8_t)(-10 * log10(err) + .5); // round it
  16 }
  17
  18 //convert transcript coordinate to chromosome coordinate and generate CIGAR string
  19 void tr2chr(const Transcript& transcript, int sp, int ep, int& pos, int& n_cigar, std::vector<uint32_t>& data) {
  20         int length = transcript.getLength();
  21         char strand = transcript.getStrand();
  22         const std::vector<Interval>& structure = transcript.getStructure();
  23
  24         int s, i;
  25         int oldlen, curlen;
  26
  27         uint32_t operation;
  28
  29         n_cigar = 0;
  30         s = structure.size();
  31
  32         if (strand == '-') {
  33                 int tmp = sp;
  34                 sp = length - ep + 1;
  35                 ep = length - tmp + 1;
  36         }
  37
  38         if (ep < 1 || sp > length) { // a read which align to polyA tails totally!
  39           pos = (sp > length ? structure[s - 1].end : structure[0].start - 1); // 0 based
  40
  41           n_cigar = 1;
  42           operation = (ep - sp + 1) << BAM_CIGAR_SHIFT | BAM_CINS; //BAM_CSOFT_CLIP;
  43           data.push_back(operation);
  44
  45           return;
  46         }
  47
  48         if (sp < 1) {
  49                 n_cigar++;
  50                 operation = (1 - sp) << BAM_CIGAR_SHIFT | BAM_CINS; //BAM_CSOFT_CLIP;
  51                 data.push_back(operation);
  52                 sp = 1;
  53         }
  54
  55         oldlen = curlen = 0;
  56
  57         for (i = 0; i < s; i++) {
  58                 oldlen = curlen;
  59                 curlen += structure[i].end - structure[i].start + 1;
  60                 if (curlen >= sp) break;
  61         }
  62         assert(i < s);
  63         pos = structure[i].start + (sp - oldlen - 1) - 1; // 0 based
  64
  65         while (curlen < ep && i < s) {
  66                 n_cigar++;
  67                 operation = (curlen - sp + 1) << BAM_CIGAR_SHIFT | BAM_CMATCH;
  68                 data.push_back(operation);
  69                 ++i;
  70                 if (i >= s) continue;
  71                 n_cigar++;
  72                 operation = (structure[i].start - structure[i - 1].end - 1) << BAM_CIGAR_SHIFT | BAM_CREF_SKIP;
  73                 data.push_back(operation);
  74
  75                 oldlen = curlen;
  76                 sp = oldlen + 1;
  77                 curlen += structure[i].end - structure[i].start + 1;
  78         }
  79
  80         if (i >= s) {
  81                 n_cigar++;
  82                 operation = (ep - length) << BAM_CIGAR_SHIFT | BAM_CINS; //BAM_CSOFT_CLIP;
  83                 data.push_back(operation);
  84         }
  85         else {
  86                 n_cigar++;
  87                 operation = (ep - sp + 1) << BAM_CIGAR_SHIFT | BAM_CMATCH;
  88                 data.push_back(operation);
  89         }
  90 }
  91
  92 #endif /* SAM_RSEM_CVT_H_ */