- std::map<std::string, int> refmap;
- std::map<std::string, int>::iterator iter;
-
- struct SingleEndT {
- bam1_t *b;
-
- SingleEndT(bam1_t *b = NULL) {
- this->b = b;
- }
-
- bool operator< (const SingleEndT& o) const {
- int strand1, strand2;
- uint32_t *p1, *p2;
-
- if (b->core.tid != o.b->core.tid) return b->core.tid < o.b->core.tid;
- if (b->core.pos != o.b->core.pos) return b->core.pos < o.b->core.pos;
- strand1 = b->core.flag & 0x0010; strand2 = o.b->core.flag & 0x0010;
- if (strand1 != strand2) return strand1 < strand2;
- if (b->core.n_cigar != o.b->core.n_cigar) return b->core.n_cigar < o.b->core.n_cigar;
- p1 = bam1_cigar(b); p2 = bam1_cigar(o.b);
- for (int i = 0; i < (int)b->core.n_cigar; i++) {
- if (*p1 != *p2) return *p1 < *p2;
- ++p1; ++p2;
- }
- return false;
- }
- };
-
- //b is mate 1, b2 is mate 2
- struct PairedEndT {
- bam1_t *b, *b2;
-
- PairedEndT() { b = NULL; b2 = NULL;}
-
- PairedEndT(bam1_t *b, bam1_t *b2) {
- this->b = b;
- this->b2 = b2;
- }
-
- bool operator< (const PairedEndT& o) const {
- int strand1, strand2;
- uint32_t *p1, *p2;
-
- //compare b
- if (b->core.tid != o.b->core.tid) return b->core.tid < o.b->core.tid;
- if (b->core.pos != o.b->core.pos) return b->core.pos < o.b->core.pos;
- strand1 = b->core.flag & 0x0010; strand2 = o.b->core.flag & 0x0010;
- if (strand1 != strand2) return strand1 < strand2;
- if (b->core.n_cigar != o.b->core.n_cigar) return b->core.n_cigar < o.b->core.n_cigar;
- p1 = bam1_cigar(b); p2 = bam1_cigar(o.b);
- for (int i = 0; i < (int)b->core.n_cigar; i++) {
- if (*p1 != *p2) return *p1 < *p2;
- ++p1; ++p2;
- }
-
- //compare b2
- if (b2->core.tid != o.b2->core.tid) return b2->core.tid < o.b2->core.tid;
- if (b2->core.pos != o.b2->core.pos) return b2->core.pos < o.b2->core.pos;
- strand1 = b2->core.flag & 0x0010; strand2 = o.b2->core.flag & 0x0010;
- if (strand1 != strand2) return strand1 < strand2;
- if (b2->core.n_cigar != o.b2->core.n_cigar) return b2->core.n_cigar < o.b2->core.n_cigar;
- p1 = bam1_cigar(b2); p2 = bam1_cigar(o.b2);
- for (int i = 0; i < (int)b2->core.n_cigar; i++) {
- if (*p1 != *p2) return *p1 < *p2;
- ++p1; ++p2;
- }
-
- return false;
- }
- };
-
- uint8_t getMAPQ(double val) {
- double err = 1.0 - val;
- if (err <= 1e-10) return 100;
- return (uint8_t)(-10 * log10(err) + .5); // round it
- }
-
- void push_qname(const uint8_t* qname, int l_qname, std::vector<uint8_t>& data) {
- for (int i = 0; i < l_qname; i++) data.push_back(*(qname + i));
- }
-
- void push_seq(const uint8_t* seq, int readlen, char strand, std::vector<uint8_t>& data) {
- int seq_len = (readlen + 1) / 2;
-
- switch (strand) {
- case '+': for (int i = 0; i < seq_len; i++) data.push_back(*(seq + i)); break;
- case '-':
- uint8_t code, base;
- code = 0; base = 0;
- for (int i = 0; i < readlen; i++) {
- switch (bam1_seqi(seq, readlen - i - 1)) {
- case 1: base = 8; break;
- case 2: base = 4; break;
- case 4: base = 2; break;
- case 8: base = 1; break;
- case 15: base = 15; break;
- default: assert(false);
- }
- code |= base << (4 * (1 - i % 2));
- if (i % 2 == 1) { data.push_back(code); code = 0; }
- }
-
- if (readlen % 2 == 1) { data.push_back(code); }
- break;
- default: assert(false);
- }
- }
-
- void push_qual(const uint8_t* qual, int readlen, char strand, std::vector<uint8_t>& data) {
- switch (strand) {
- case '+': for (int i = 0; i < readlen; i++) data.push_back(*(qual + i)); break;
- case '-': for (int i = readlen - 1; i >= 0; i--) data.push_back(*(qual + i)); break;
- default: assert(false);
- }
- }
-
- //convert transcript coordinate to chromosome coordinate and generate CIGAR string
- void tr2chr(const Transcript&, int, int, int&, int&, std::vector<uint8_t>&);