hit_set.h

   1 /*
   2  * hit_set.h
   3  *
   4  *  Created on: Jul 31, 2009
   5  *      Author: Ben Langmead
   6  */
   7
   8 #ifndef HIT_SET_H_
   9 #define HIT_SET_H_
  10
  11 #include <seqan/sequence.h>
  12 #include <vector>
  13 #include <algorithm>
  14 #include "assert_helpers.h"
  15 #include "filebuf.h"
  16 #include "edit.h"
  17 #include "alphabet.h"
  18 #include "annot.h"
  19 #include "refmap.h"
  20
  21 /**
  22  * Encapsulates a hit contained within a HitSet that can be
  23  * (de)serialized to/from FileBufs.  Used for chaining.
  24  */
  25 struct HitSetEnt {
  26         typedef std::pair<uint32_t,uint32_t> U32Pair;
  27
  28         HitSetEnt() { }
  29
  30         /**
  31          * Write binary representation of HitSetEnt to an OutFileBuf.
  32          */
  33         void serialize(OutFileBuf& fb) const {
  34                 fb.writeChars((const char*)&h.first, 4);
  35                 fb.writeChars((const char*)&h.second, 4);
  36                 assert(fw == 0 || fw == 1);
  37                 fb.write(fw);
  38                 assert_geq(stratum, 0);
  39                 assert_lt(stratum, 4);
  40                 fb.write(stratum);
  41                 assert_eq(stratum, (cost >> 14));
  42                 fb.writeChars((const char*)&cost, 2);
  43                 fb.writeChars((const char*)&oms, 4);
  44                 uint32_t sz = edits.size();
  45                 fb.writeChars((const char*)&sz, 4);
  46                 std::vector<Edit>::const_iterator it;
  47                 for(it = edits.begin(); it != edits.end(); it++) {
  48                         it->serialize(fb);
  49                 }
  50                 sz = cedits.size();
  51                 fb.writeChars((const char*)&sz, 4);
  52                 for(it = cedits.begin(); it != cedits.end(); it++) {
  53                         it->serialize(fb);
  54                 }
  55         }
  56
  57         /**
  58          * Repopulate a HitSetEnt from its binary representation in FileBuf.
  59          */
  60         void deserialize(FileBuf& fb) {
  61                 fb.get((char*)&h.first, 4);
  62                 fb.get((char*)&h.second, 4);
  63                 fw = fb.get();
  64                 assert(fw == 0 || fw == 1);
  65                 stratum = fb.get();
  66                 assert_geq(stratum, 0);
  67                 assert_lt(stratum, 4);
  68                 fb.get((char*)&cost, 2);
  69                 assert_eq(stratum, (cost >> 14));
  70                 fb.get((char*)&oms, 4);
  71                 uint32_t sz = 0;
  72                 fb.get((char*)&sz, 4);
  73                 assert_lt(sz, 1024);
  74                 edits.resize(sz);
  75                 for(uint32_t i = 0; i < sz; i++) {
  76                         edits[i].deserialize(fb);
  77                 }
  78                 fb.get((char*)&sz, 4);
  79                 assert_lt(sz, 1024);
  80                 cedits.resize(sz);
  81                 for(uint32_t i = 0; i < sz; i++) {
  82                         cedits[i].deserialize(fb);
  83                 }
  84         }
  85
  86         /**
  87          * Less than operator.  Break HitSetEnt ties by:
  88          *  - Stratum, then
  89          *  - Cost, then
  90          *  - Position, then
  91          *  - Orientation
  92          */
  93         int operator< (const HitSetEnt &rhs) const {
  94                 if(stratum < rhs.stratum) return 1;
  95                 if(stratum > rhs.stratum) return 0;
  96                 if(cost < rhs.cost) return 1;
  97                 if(cost > rhs.cost) return 0;
  98                 if(h < rhs.h) return 1;
  99                 if(h > rhs.h) return 0;
 100                 return (fw < rhs.fw)? 1 : 0;
 101         }
 102
 103         /**
 104          * Greater than operator.
 105          */
 106         int operator> (const HitSetEnt &rhs) const {
 107                 if(stratum < rhs.stratum) return 0;
 108                 if(stratum > rhs.stratum) return 1;
 109                 if(cost < rhs.cost) return 0;
 110                 if(cost > rhs.cost) return 1;
 111                 if(h < rhs.h) return 0;
 112                 if(h > rhs.h) return 1;
 113                 return (fw <= rhs.fw)? 0 : 1;
 114         }
 115
 116         /**
 117          * Equality comparison operator.
 118          */
 119         int operator== (const HitSetEnt &rhs) const {
 120                 return(stratum == rhs.stratum &&
 121                        cost == rhs.cost &&
 122                        fw == rhs.fw &&
 123                        h == rhs.h);
 124         }
 125
 126         /**
 127          * Indexing returns edits.
 128          */
 129         Edit& operator[](unsigned x) {
 130                 return edits[x];
 131         }
 132
 133         /**
 134          * Indexing returns edits.
 135          */
 136         const Edit& operator[](unsigned x) const {
 137                 return edits[x];
 138         }
 139
 140         /**
 141          * Another way to get at an edit.
 142          */
 143         Edit& editAt(unsigned i) {
 144                 return edits[i];
 145         }
 146
 147         /**
 148          * Another way to get at a const edit.
 149          */
 150         const Edit& editAt(unsigned i) const {
 151                 return edits[i];
 152         }
 153
 154         /**
 155          * Get the ith color edit.
 156          */
 157         Edit& colorEditAt(unsigned i) {
 158                 return cedits[i];
 159         }
 160
 161         /**
 162          * Another way to get at an edit.
 163          */
 164         const Edit& colorEditAt(unsigned i) const {
 165                 return cedits[i];
 166         }
 167
 168         /**
 169          * Return the front entry.
 170          */
 171         Edit& front() {
 172                 return edits.front();
 173         }
 174
 175         /**
 176          * Return the back entry.
 177          */
 178         Edit& back() {
 179                 return edits.back();
 180         }
 181
 182         /**
 183          * Expand the entry list by one.
 184          */
 185         void expand() {
 186                 edits.resize(edits.size() + 1);
 187         }
 188
 189         /**
 190          * Sort edits by position
 191          */
 192         void sort() {
 193                 if(edits.size() > 1) std::sort(edits.begin(), edits.end());
 194         }
 195
 196         /**
 197          * Return number of edits.
 198          */
 199         size_t size() const {
 200                 return edits.size();
 201         }
 202
 203         bool empty() const {
 204                 return edits.empty();
 205         }
 206
 207         /**
 208          * Write HitSetEnt to an output stream.
 209          */
 210         friend std::ostream& operator << (std::ostream& os, const HitSetEnt& hse);
 211
 212         U32Pair h; // reference coordinates
 213         uint8_t fw; // orientation
 214         int8_t stratum; // stratum
 215         uint16_t cost; // cost, including stratum
 216         uint32_t oms; // # others
 217         std::vector<Edit> edits; // edits to get from reference to subject
 218         std::vector<Edit> cedits; // color edits to get from reference to subject
 219 };
 220
 221 /**
 222  * Encapsulates a set of hits that can be (de)serialized to/from
 223  * FileBufs.  Used for chaining.
 224  */
 225 struct HitSet {
 226
 227         typedef std::vector<HitSetEnt> EntVec;
 228         typedef EntVec::const_iterator Iter;
 229         typedef std::pair<uint32_t,uint32_t> U32Pair;
 230
 231         HitSet() {
 232                 maxedStratum = -1;
 233         }
 234
 235         HitSet(FileBuf& fb) {
 236                 deserialize(fb);
 237         }
 238
 239         /**
 240          * Write binary representation of HitSet to an OutFileBuf.
 241          */
 242         void serialize(OutFileBuf& fb) const {
 243                 fb.write(color ? 1 : 0);
 244                 uint32_t i = seqan::length(name);
 245                 assert_gt(i, 0);
 246                 fb.writeChars((const char*)&i, 4);
 247                 fb.writeChars(seqan::begin(name), i);
 248                 i = seqan::length(seq);
 249                 assert_gt(i, 0);
 250                 assert_lt(i, 1024);
 251                 fb.writeChars((const char*)&i, 4);
 252                 for(size_t j = 0; j < i; j++) {
 253                         fb.write("ACGTN"[(int)seq[j]]);
 254                 }
 255                 fb.writeChars(seqan::begin(qual), i);
 256                 i = ents.size();
 257                 fb.writeChars((const char*)&i, 4);
 258                 std::vector<HitSetEnt>::const_iterator it;
 259                 for(it = ents.begin(); it != ents.end(); it++) {
 260                         it->serialize(fb);
 261                 }
 262                 fb.write(maxedStratum);
 263         }
 264
 265         /**
 266          * Repopulate a HitSet from its binary representation in FileBuf.
 267          */
 268         void deserialize(FileBuf& fb) {
 269                 color = (fb.get() != 0 ? true : false);
 270                 uint32_t sz = 0;
 271                 if(fb.get((char*)&sz, 4) != 4) {
 272                         seqan::clear(name);
 273                         seqan::clear(seq);
 274                         return;
 275                 }
 276                 assert_gt(sz, 0);
 277                 assert_lt(sz, 1024);
 278                 seqan::resize(name, sz);
 279                 fb.get(seqan::begin(name), sz);
 280                 fb.get((char*)&sz, 4);
 281                 assert_gt(sz, 0);
 282                 assert_lt(sz, 1024);
 283                 seqan::resize(seq, sz);
 284                 for(size_t j = 0; j < sz; j++) {
 285                         seq[j] = charToDna5[fb.get()];
 286                 }
 287                 seqan::resize(qual, sz);
 288                 fb.get(seqan::begin(qual), sz);
 289                 fb.get((char*)&sz, 4);
 290                 if(sz > 0) {
 291                         ents.resize(sz);
 292                         for(size_t i = 0; i < sz; i++) {
 293                                 ents[i].deserialize(fb);
 294                         }
 295                 } else {
 296                         ents.clear();
 297                 }
 298                 maxedStratum = fb.get();
 299         }
 300
 301         /**
 302          * Return true iff this HitSet is initialized with a read (but not
 303          * necessarily any alignments).
 304          */
 305         bool initialized() const {
 306                 return !seqan::empty(seq);
 307         }
 308
 309         /**
 310          * Return true iff this HitSet has no hits.
 311          */
 312         bool empty() const {
 313                 return ents.empty();
 314         }
 315
 316         /**
 317          * Return number of entries in this HitSet.
 318          */
 319         size_t size() const {
 320                 return ents.size();
 321         }
 322
 323         /**
 324          * Remove all entries from this HitSet.
 325          */
 326         void clear() {
 327                 maxedStratum = -1;
 328                 ents.clear();
 329         }
 330
 331         Iter begin() const {
 332                 return ents.begin();
 333         }
 334
 335         Iter end() const {
 336                 return ents.end();
 337         }
 338
 339         /**
 340          * Return the front entry.
 341          */
 342         HitSetEnt& front() {
 343                 return ents.front();
 344         }
 345
 346         /**
 347          * Return the back entry.
 348          */
 349         HitSetEnt& back() {
 350                 return ents.back();
 351         }
 352
 353         /**
 354          * Expand the entry list by one.
 355          */
 356         void expand() {
 357                 ents.resize(ents.size() + 1);
 358                 assert(ents.back().empty());
 359         }
 360
 361         /**
 362          * Resize entry list
 363          */
 364         void resize(size_t sz) {
 365                 ents.resize(sz);
 366         }
 367
 368         /**
 369          * Sort hits by stratum/penalty.
 370          */
 371         void sort() {
 372                 if(ents.size() > 1) std::sort(ents.begin(), ents.end());
 373         }
 374
 375         /**
 376          * Return true iff Ents are sorted
 377          */
 378         bool sorted() const {
 379                 if(ents.empty()) return true;
 380                 for(size_t i = 0; i < ents.size()-1; i++) {
 381                         if(!(ents[i] < ents[i+1])) return false;
 382                 }
 383                 return true;
 384         }
 385
 386         /**
 387          * Remove the ith hit from the HitSet, shifting all subsequent hits
 388          * up by one.
 389          */
 390         void remove(size_t i) {
 391                 ents.erase(ents.begin() + i);
 392         }
 393
 394         /**
 395          * Return true if strata are uniform across hits; assert if they're
 396          * not.
 397          */
 398         bool uniformStrata() const {
 399                 for(size_t i = 0; i < ents.size()-1; i++) {
 400                         assert(ents[i].stratum == ents[i+1].stratum);
 401                 }
 402                 return true;
 403         }
 404
 405         /**
 406          * Indexing returns entries.
 407          */
 408         HitSetEnt& operator[](unsigned x) {
 409                 return ents[x];
 410         }
 411
 412         /**
 413          * Indexing returns entries.
 414          */
 415         const HitSetEnt& operator[](unsigned x) const {
 416                 return ents[x];
 417         }
 418
 419         /**
 420          * Apply a reference mappings to all the contained hits.
 421          */
 422         void applyReferenceMap(const ReferenceMap& map) {
 423                 std::vector<HitSetEnt>::iterator it;
 424                 for(it = ents.begin(); it != ents.end(); it++) map.map(it->h);
 425         }
 426
 427         /**
 428          * Clear out all the strings and all the entries.
 429          */
 430         void clearAll() {
 431                 seqan::clear(name);
 432                 seqan::clear(seq);
 433                 seqan::clear(qual);
 434                 ents.clear();
 435                 color = false;
 436         }
 437
 438         /**
 439          *
 440          */
 441         bool tryReplacing(U32Pair h,
 442                           bool fw,
 443                           uint16_t cost,
 444                           size_t& pos)
 445         {
 446                 for(size_t i = 0; i < ents.size(); i++) {
 447                         if(ents[i].h == h && ents[i].fw == fw) {
 448                                 if(cost < ents[i].cost) {
 449                                         // New hit at same position is better!  Replace it.
 450                                         ents[i].h = h;
 451                                         ents[i].fw = fw;
 452                                         ents[i].stratum = cost >> 14;
 453                                         ents[i].cost = cost;
 454                                         pos = i;
 455                                         return true;
 456                                 } else {
 457                                         pos = 0xffffffff;
 458                                         return true;
 459                                 }
 460                         }
 461                 }
 462                 return false;
 463         }
 464
 465         /**
 466          * Report up to 'khits' hits from this HitSet.
 467          */
 468         void reportUpTo(std::ostream& os, int khits = INT_MAX);
 469
 470         /**
 471          * Print this HitSet.
 472          */
 473         friend std::ostream& operator << (std::ostream& os, const HitSet& hs);
 474
 475         seqan::String<char> name;
 476         seqan::String<seqan::Dna5> seq;
 477         seqan::String<char> qual;
 478         int8_t maxedStratum;
 479         std::vector<HitSetEnt> ents;
 480         bool color; // whether read was orginally in colorspace
 481 };
 482
 483 #endif /* HIT_SET_H_ */