1 #ifndef SAM_RSEM_CVT_H_
2 #define SAM_RSEM_CVT_H_
9 #include "Transcript.h"
10 #include "Transcripts.h"
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
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();
35 ep = length - tmp + 1;
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
42 operation = (ep - sp + 1) << BAM_CIGAR_SHIFT | BAM_CINS; //BAM_CSOFT_CLIP;
43 data.push_back(operation);
50 operation = (1 - sp) << BAM_CIGAR_SHIFT | BAM_CINS; //BAM_CSOFT_CLIP;
51 data.push_back(operation);
57 for (i = 0; i < s; i++) {
59 curlen += structure[i].end - structure[i].start + 1;
60 if (curlen >= sp) break;
63 pos = structure[i].start + (sp - oldlen - 1) - 1; // 0 based
65 while (curlen < ep && i < s) {
67 operation = (curlen - sp + 1) << BAM_CIGAR_SHIFT | BAM_CMATCH;
68 data.push_back(operation);
72 operation = (structure[i].start - structure[i - 1].end - 1) << BAM_CIGAR_SHIFT | BAM_CREF_SKIP;
73 data.push_back(operation);
77 curlen += structure[i].end - structure[i].start + 1;
82 operation = (ep - length) << BAM_CIGAR_SHIFT | BAM_CINS; //BAM_CSOFT_CLIP;
83 data.push_back(operation);
87 operation = (ep - sp + 1) << BAM_CIGAR_SHIFT | BAM_CMATCH;
88 data.push_back(operation);
92 #endif /* SAM_RSEM_CVT_H_ */