bc_aux.h

   1 #ifndef BC_AUX_H_
   2 #define BC_AUX_H_
   3
   4 #include<map>
   5
   6 #include <stdint.h>
   7 #include "sam/bam.h"
   8
   9 struct SingleEndT {
  10         bam1_t *b;
  11
  12         SingleEndT(bam1_t *b) {
  13                 this->b = b;
  14         }
  15
  16         int getSign(bool value) const { return value ? -1 : 1; }
  17
  18         int compare(const SingleEndT& o) const {
  19                 int strand1, strand2;
  20                 uint32_t *p1, *p2;
  21
  22                 if (b->core.tid != o.b->core.tid) return getSign(b->core.tid < o.b->core.tid);
  23                 if (b->core.pos != o.b->core.pos) return getSign(b->core.pos < o.b->core.pos);
  24                 strand1 = b->core.flag & 0x0010; strand2 = o.b->core.flag & 0x0010;
  25                 if (strand1 != strand2) return getSign(strand1 < strand2);
  26                 if (b->core.n_cigar != o.b->core.n_cigar) return getSign(b->core.n_cigar < o.b->core.n_cigar);
  27                 p1 = bam1_cigar(b); p2 = bam1_cigar(o.b);
  28                 for (int i = 0; i < (int)b->core.n_cigar; i++) {
  29                         if (*p1 != *p2) return getSign(*p1 < *p2);
  30                         ++p1; ++p2;
  31                 }
  32
  33                 return 0;
  34         }
  35
  36         bool operator< (const SingleEndT& o) const {
  37                 return compare(o) < 0;
  38         }
  39 };
  40
  41 struct PairedEndT {
  42         SingleEndT mate1, mate2;
  43
  44         PairedEndT(const SingleEndT& mate1, const SingleEndT& mate2) : mate1(mate1), mate2(mate2) {
  45         }
  46
  47         bool operator< (const PairedEndT& o) const {
  48                 int value = mate1.compare(o.mate1);
  49                 return value < 0 || (value == 0 && mate2 < o.mate2);
  50         }
  51 };
  52
  53 class CollapseMap {
  54 public:
  55         CollapseMap() { isPaired = false; smap.clear(); pmap.clear(); }
  56
  57         void init(bool isPaired) {
  58                 this->isPaired = isPaired;
  59                 isPaired ? pmap.clear() : smap.clear();
  60         }
  61
  62         void insert(bam1_t *b, bam1_t *b2, float prb) {
  63                 if (!isPaired) {
  64                         smapIter = smap.find(SingleEndT(b));
  65                         if (smapIter == smap.end()) { smap[SingleEndT(bam_dup1(b))] = prb; }
  66                         else smapIter->second += prb;
  67                 }
  68                 else {
  69                         pmapIter = pmap.find(PairedEndT(SingleEndT(b), SingleEndT(b2)));
  70                         if (pmapIter == pmap.end()) { pmap[PairedEndT(SingleEndT(bam_dup1(b)), SingleEndT(bam_dup1(b2)))] = prb; }
  71                         else pmapIter->second += prb;
  72                 }
  73         }
  74
  75         //once this function is called, "insert" cannot be called anymore
  76         bool empty(bool& par) {
  77                 bool value;
  78
  79                 par = isPaired;
  80                 if (!isPaired) { value = smap.empty(); smapIter = smap.begin(); }
  81                 else { value = pmap.empty(); pmapIter = pmap.begin(); }
  82
  83                 return value;
  84         }
  85
  86         bool next(bam1_t*& b, bam1_t*& b2, float& prb) {
  87                 bool value;
  88
  89                 if (!isPaired) {
  90                         value = smapIter != smap.end();
  91                         if (value) {
  92                                 b = smapIter->first.b;
  93                                 prb = smapIter->second;
  94                                 smapIter++;
  95                         }
  96                 }
  97                 else {
  98                         value = pmapIter != pmap.end();
  99                         if (value) {
 100                                 b = pmapIter->first.mate1.b;
 101                                 b2 = pmapIter->first.mate2.b;
 102                                 prb = pmapIter->second;
 103                                 pmapIter++;
 104                         }
 105                 }
 106
 107                 return value;
 108         }
 109
 110 private:
 111         bool isPaired;
 112
 113         std::map<SingleEndT, float> smap;
 114         std::map<SingleEndT, float>::iterator smapIter;
 115
 116         std::map<PairedEndT, float> pmap;
 117         std::map<PairedEndT, float>::iterator pmapIter;
 118 };
 119
 120 #endif /* BC_AUX_H_ */